source: LMDZ6/branches/Amaury_dev/libf/phylmd/phys_output_mod.F90 @ 5441

Last change on this file since 5441 was 5160, checked in by abarral, 5 months ago

Put .h into modules

  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 29.3 KB
Line 
1! $Id: phys_output_mod.F90 5160 2024-08-03 12:56:58Z fhourdin $
2
3MODULE phys_output_mod
4  USE indice_sol_mod
5  USE phys_output_var_mod
6  USE phys_output_write_mod, ONLY: phys_output_write
7  REAL, DIMENSION(nfiles), SAVE :: ecrit_files
8
9  ! Abderrahmane 12 2007
10  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
11  !!! Ecreture des Sorties du modele dans les fichiers Netcdf :
12  ! histmth.nc : moyennes mensuelles
13  ! histday.nc : moyennes journalieres
14  ! histhf.nc  : moyennes toutes les 3 heures
15  ! histins.nc : valeurs instantanees
16  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
17
18CONTAINS
19
20  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
21  !!!!!!!!! Ouverture des fichier et definition des variable de sortie !!!!!!!!
22  !! histbeg, histvert et histdef
23  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
24
25  SUBROUTINE phys_output_open(rlon, rlat, pim, tabij, ipt, jpt, plon, plat, &
26          jjmp1, nlevSTD, clevSTD, rlevSTD, dtime, ok_veget, &
27          type_ocean, iflag_pbl, iflag_pbl_split, ok_mensuel, ok_journe, &
28          ok_hf, ok_instan, ok_LES, ok_ade, ok_aie, read_climoz, &
29          phys_out_filestations, &
30          aerosol_couple, flag_aerosol_strat, &
31          pdtphys, paprs, pphis, pplay, lmax_th, ptconv, ptconvth, ivap, &
32          d_u, d_t, qx, d_qx, zmasse, ok_sync)
33
34    USE iophy
35    USE dimphy
36    USE infotrac_phy, ONLY: nqtot, tracers, niso, ntraciso => ntiso
37    USE lmdz_strings, ONLY: maxlen
38    USE ioipsl
39    USE phys_cal_mod, ONLY: hour, calend
40    USE lmdz_phys_para
41    !Martin
42    USE surface_data, ONLY: landice_opt
43    USE phys_output_ctrlout_mod
44    USE lmdz_grid_phy, ONLY: klon_glo, nbp_lon, nbp_lat
45    USE lmdz_print_control, ONLY: prt_level, lunout
46    USE lmdz_vertical_layers, ONLY: ap, bp, preff, presnivs, aps, bps, pseudoalt, presinter
47    USE time_phylmdz_mod, ONLY: day_ini, itau_phy, start_time, annee_ref, day_ref
48    ! ug Pour les sorties XIOS
49    USE lmdz_wxios
50    USE infotrac_phy, ONLY: nbtr_bin
51#ifdef ISO
52    USE isotopes_mod, ONLY: isoName,iso_HTO
53#ifdef ISOTRAC
54    USE isotrac_mod, ONLY: index_zone,index_iso,strtrac
55#endif
56#endif
57    USE lmdz_cppkeys_wrapper, ONLY: CPPKEY_STRATAER
58    USE lmdz_clesphys
59    USE lmdz_yomcst
60
61    IMPLICIT NONE
62
63    ! ug Nouveaux arguments n\'ecessaires au histwrite_mod:
64    INTEGER, INTENT(IN) :: ivap
65    INTEGER, DIMENSION(klon), INTENT(IN) :: lmax_th
66    LOGICAL, INTENT(IN) :: ok_sync
67    LOGICAL, DIMENSION(klon, klev), INTENT(IN) :: ptconv, ptconvth
68    REAL, INTENT(IN) :: pdtphys
69    REAL, DIMENSION(klon), INTENT(IN) :: pphis
70    REAL, DIMENSION(klon, klev), INTENT(IN) :: pplay, d_u, d_t
71    REAL, DIMENSION(klon, klev + 1), INTENT(IN) :: paprs
72    REAL, DIMENSION(klon, klev, nqtot), INTENT(IN) :: qx, d_qx
73    REAL, DIMENSION(klon, klev), INTENT(IN) :: zmasse
74
75    REAL, DIMENSION(klon), INTENT(IN) :: rlon
76    REAL, DIMENSION(klon), INTENT(IN) :: rlat
77    INTEGER, INTENT(IN) :: pim
78    INTEGER, DIMENSION(pim) :: tabij
79    INTEGER, DIMENSION(pim), INTENT(IN) :: ipt, jpt
80    REAL, DIMENSION(pim), INTENT(IN) :: plat, plon
81    REAL, DIMENSION(pim, 2) :: plat_bounds, plon_bounds
82
83    INTEGER :: jjmp1
84    INTEGER :: nlevSTD, radpas
85    LOGICAL :: ok_mensuel, ok_journe, ok_hf, ok_instan
86    LOGICAL :: ok_LES, ok_ade, ok_aie
87    INTEGER :: flag_aerosol_strat
88    LOGICAL :: aerosol_couple
89    INTEGER, INTENT(IN) :: read_climoz ! read ozone climatology
90    !     Allowed values are 0, 1 and 2
91    !     0: do not read an ozone climatology
92    !     1: read a single ozone climatology that will be used day and night
93    !     2: read two ozone climatologies, the average day and night
94    !     climatology and the daylight climatology
95
96    REAL :: dtime
97    INTEGER :: idayref
98    REAL :: zjulian_start, zjulian
99    CHARACTER(LEN = 4), DIMENSION(nlevSTD) :: clevSTD
100    REAL, DIMENSION(nlevSTD) :: rlevSTD
101    INTEGER :: nsrf, k, iq, iff, i, j, ilev, itr, itrb, ixt, iiso, izone
102    INTEGER :: naero
103    LOGICAL :: ok_veget
104    INTEGER :: iflag_pbl
105    INTEGER :: iflag_pbl_split
106    CHARACTER(LEN = 4) :: bb2
107    CHARACTER(LEN = 2) :: bb3
108    CHARACTER(LEN = 6) :: type_ocean
109    INTEGER, DIMENSION(nbp_lon * jjmp1) :: ndex2d
110    INTEGER, DIMENSION(nbp_lon * jjmp1 * klev) :: ndex3d
111    INTEGER :: imin_ins, imax_ins
112    INTEGER :: jmin_ins, jmax_ins
113    INTEGER, DIMENSION(nfiles) :: phys_out_levmin, phys_out_levmax
114    INTEGER, DIMENSION(nfiles) :: phys_out_filelevels
115    CHARACTER(LEN = 20), DIMENSION(nfiles) :: chtimestep = (/ 'Default', 'Default', 'Default', 'Default', 'Default', &
116            'Default', 'Default', 'Default', 'Default', 'Default' /)
117    LOGICAL, DIMENSION(nfiles) :: phys_out_filekeys
118    LOGICAL, DIMENSION(nfiles) :: phys_out_filestations
119
120#ifdef ISO
121    CHARACTER(LEN=maxlen) :: outiso
122    CHARACTER(LEN=20) :: unit
123#endif
124    CHARACTER(LEN = maxlen) :: tnam, lnam, dn
125    INTEGER :: flag(nfiles)
126
127    !!!!!!!!!! stockage dans une region limitee pour chaque fichier !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
128    !                 entre [phys_out_lonmin,phys_out_lonmax] et [phys_out_latmin,phys_out_latmax]
129    LOGICAL, DIMENSION(nfiles), SAVE :: &
130            phys_out_regfkey = [.FALSE., .FALSE., .FALSE., .FALSE., .FALSE., .FALSE., .FALSE., .FALSE., .FALSE., .FALSE.]
131    REAL, DIMENSION(nfiles), SAVE :: &
132            phys_out_lonmin = [  -180., -180., -180., -180., -180., -180., -180., -180., -180., -180.], &
133            phys_out_lonmax = [   180., 180., 180., 180., 180., 180., 180., 180., 180., 180.], &
134            phys_out_latmin = [   -90., -90., -90., -90., -90., -90., -90., -90., -90., -90.], &
135            phys_out_latmax = [    90., 90., 90., 90., 90., 90., 90., 90., 90., 90.]
136    REAL, DIMENSION(klev, 2) :: Ahyb_bounds, Bhyb_bounds
137    REAL, DIMENSION(klev + 1) :: lev_index
138
139    ! ug Variables utilis\'ees pour r\'ecup\'erer le calendrier pour xios
140    INTEGER :: x_an, x_mois, x_jour
141    REAL :: x_heure
142    INTEGER :: ini_an, ini_mois, ini_jour
143    REAL :: ini_heure
144
145    INTEGER :: ISW
146    REAL, DIMENSION(NSW) :: wl1_sun, wl2_sun !wavelength bounds (in um) for SW
147    REAL, DIMENSION(NSW) :: wn1_sun, wn2_sun !wavenumber bounds (in m-1) for SW
148    REAL, DIMENSION(NSW) :: spectband  !mean wavenumb. of each sp.band
149    REAL, DIMENSION(NSW, 2) :: spbnds_sun !bounds of spectband
150
151    WRITE(lunout, *) 'Debut phys_output_mod.F90'
152    ! Initialisations (Valeurs par defaut
153
154    DO ilev = 1, klev
155      Ahyb_bounds(ilev, 1) = ap(ilev)
156      Ahyb_bounds(ilev, 2) = ap(ilev + 1)
157      Bhyb_bounds(ilev, 1) = bp(ilev)
158      Bhyb_bounds(ilev, 2) = bp(ilev + 1)
159      lev_index(ilev) = REAL(ilev)
160    END DO
161    lev_index(klev + 1) = REAL(klev + 1)
162
163    IF (.NOT. ALLOCATED(o_trac)) ALLOCATE(o_trac(nqtot))
164    IF (.NOT. ALLOCATED(o_trac_cum)) ALLOCATE(o_trac_cum(nqtot))
165    ALLOCATE(o_dtr_the(nqtot), o_dtr_con(nqtot), o_dtr_lessi_impa(nqtot))
166    ALLOCATE(o_dtr_lessi_nucl(nqtot), o_dtr_insc(nqtot), o_dtr_bcscav(nqtot))
167    ALLOCATE(o_dtr_evapls(nqtot), o_dtr_ls(nqtot), o_dtr_trsp(nqtot))
168    ALLOCATE(o_dtr_sscav(nqtot), o_dtr_sat(nqtot), o_dtr_uscav(nqtot))
169    ALLOCATE(o_dtr_dry(nqtot), o_dtr_vdf(nqtot))
170    IF (CPPKEY_STRATAER) THEN
171      ALLOCATE(o_nd_mode(nbtr_bin), o_sulfmmr_mode(nbtr_bin))
172    END IF
173#ifdef ISO
174    ALLOCATE(o_xtprecip(ntraciso))
175    ALLOCATE(o_xtplul(ntraciso))
176    ALLOCATE(o_xtpluc(ntraciso))
177    ALLOCATE(o_xtevap(ntraciso))
178    ALLOCATE(o_xtevap_srf(ntraciso,4))
179    ALLOCATE(o_xtovap(ntraciso))
180    ALLOCATE(o_xtoliq(ntraciso))
181    ALLOCATE(o_xtcond(ntraciso))
182    ALLOCATE(o_xtrunoff_diag(ntraciso))
183    ALLOCATE(o_dxtdyn(ntraciso))
184    ALLOCATE(o_dxtldyn(ntraciso))
185    ALLOCATE(o_dxtcon(ntraciso))
186    ALLOCATE(o_dxtlsc(ntraciso))
187    ALLOCATE(o_dxteva(ntraciso))
188    ALLOCATE(o_dxtajs(ntraciso))
189    ALLOCATE(o_dxtvdf(ntraciso))
190    ALLOCATE(o_dxtthe(ntraciso))
191    ALLOCATE(o_dxtch4(ntraciso))
192    IF (iso_HTO.gt.0) THEN
193      ALLOCATE(o_dxtprod_nucl(ntraciso))
194      ALLOCATE(o_dxtcosmo(ntraciso))
195      ALLOCATE(o_dxtdecroiss(ntraciso))
196    endif
197#endif
198
199    levmax = [klev, klev, klev, klev, klev, klev, nlevSTD, nlevSTD, nlevSTD, klev]
200
201    phys_out_filenames(1) = 'histmth'
202    phys_out_filenames(2) = 'histday'
203    phys_out_filenames(3) = 'histhf6h'
204    phys_out_filenames(4) = 'histhf3h'
205    phys_out_filenames(5) = 'histhf3hm'
206    phys_out_filenames(6) = 'histstn'
207    phys_out_filenames(7) = 'histmthNMC'
208    phys_out_filenames(8) = 'histdayNMC'
209    phys_out_filenames(9) = 'histhfNMC'
210    phys_out_filenames(10) = 'histstrataer'
211
212    type_ecri(1) = 'ave(X)'
213    type_ecri(2) = 'ave(X)'
214    type_ecri(3) = 'inst(X)'
215    type_ecri(4) = 'inst(X)'
216    type_ecri(5) = 'ave(X)'
217    type_ecri(6) = 'inst(X)'
218    type_ecri(7) = 'inst(X)'
219    type_ecri(8) = 'inst(X)'
220    type_ecri(9) = 'inst(X)'
221    type_ecri(10) = 'ave(X)'
222
223    clef_files(1) = ok_mensuel
224    clef_files(2) = ok_journe
225    clef_files(3) = ok_hf
226    clef_files(4) = ok_instan
227    clef_files(5) = ok_LES
228    clef_files(6) = ok_instan
229    clef_files(7) = ok_histNMC(1)
230    clef_files(8) = ok_histNMC(2)
231    clef_files(9) = ok_histNMC(3)
232    clef_files(10) = CPPKEY_STRATAER
233
234    !sortir des fichiers "stations" si clef_stations(:)=.TRUE.
235    clef_stations(1) = .FALSE.
236    clef_stations(2) = .FALSE.
237    clef_stations(3) = .FALSE.
238    clef_stations(4) = .FALSE.
239    clef_stations(5) = .FALSE.
240    clef_stations(6) = .FALSE.
241    clef_stations(7) = .FALSE.
242    clef_stations(8) = .FALSE.
243    clef_stations(9) = .FALSE.
244    clef_stations(10) = .FALSE.
245
246    lev_files(1) = lev_histmth
247    lev_files(2) = lev_histday
248    lev_files(3) = lev_histhf
249    lev_files(4) = lev_histins
250    lev_files(5) = lev_histLES
251    lev_files(6) = lev_histins
252    lev_files(7) = levout_histNMC(1)
253    lev_files(8) = levout_histNMC(2)
254    lev_files(9) = levout_histNMC(3)
255    lev_files(10) = 5
256
257    ecrit_files(1) = ecrit_mth
258    ecrit_files(2) = ecrit_day
259    ecrit_files(3) = ecrit_hf
260    ecrit_files(4) = ecrit_ins
261    ecrit_files(5) = ecrit_LES
262    ecrit_files(6) = ecrit_ins
263    ecrit_files(7) = freq_outNMC(1)
264    ecrit_files(8) = freq_outNMC(2)
265    ecrit_files(9) = freq_outNMC(3)
266    ecrit_files(10) = ecrit_mth
267
268    !! Lectures des parametres de sorties dans physiq.def
269
270    CALL getin('phys_out_regfkey', phys_out_regfkey)
271    CALL getin('phys_out_lonmin', phys_out_lonmin)
272    CALL getin('phys_out_lonmax', phys_out_lonmax)
273    CALL getin('phys_out_latmin', phys_out_latmin)
274    CALL getin('phys_out_latmax', phys_out_latmax)
275    phys_out_levmin(:) = levmin(:)
276    CALL getin('phys_out_levmin', levmin)
277    phys_out_levmax(:) = levmax(:)
278    CALL getin('phys_out_levmax', levmax)
279    CALL getin('phys_out_filenames', phys_out_filenames)
280    phys_out_filekeys(:) = clef_files(:)
281    CALL getin('phys_out_filekeys', clef_files)
282    phys_out_filestations(:) = clef_stations(:)
283    CALL getin('phys_out_filestations', clef_stations)
284    phys_out_filelevels(:) = lev_files(:)
285    CALL getin('phys_out_filelevels', lev_files)
286    CALL getin('phys_out_filetimesteps', chtimestep)
287    phys_out_filetypes(:) = type_ecri(:)
288    CALL getin('phys_out_filetypes', type_ecri)
289
290    type_ecri_files(:) = type_ecri(:)
291
292    !    if (ok_all_xml) phys_out_filelevels = 999
293
294    WRITE(lunout, *)'phys_out_lonmin=', phys_out_lonmin
295    WRITE(lunout, *)'phys_out_lonmax=', phys_out_lonmax
296    WRITE(lunout, *)'phys_out_latmin=', phys_out_latmin
297    WRITE(lunout, *)'phys_out_latmax=', phys_out_latmax
298    WRITE(lunout, *)'phys_out_filenames=', phys_out_filenames
299    WRITE(lunout, *)'phys_out_filetypes=', type_ecri
300    WRITE(lunout, *)'phys_out_filekeys=', clef_files
301    WRITE(lunout, *)'phys_out_filestations=', clef_stations
302    WRITE(lunout, *)'phys_out_filelevels=', lev_files
303    WRITE(lunout, *)'phys_out_regfkey=', phys_out_regfkey
304
305    ! A noter pour
306    ! l heure initiale - dans les fichiers histoire hist* - on met comme
307    ! heure de debut soit la vraie heure (pour le 1D) soit 0h (pour le 3D)
308    ! afin d avoir une seule sortie mensuelle par mois lorsque l on tourne
309    ! par annee (IM).
310
311    idayref = day_ref
312    IF (klon_glo==1) THEN
313      ! current_time (used to compute hour) is updated at the begining of
314      ! the physics; to set the correct outputs "initial time" we thus
315      ! have to use (hour-dtphys).
316      CALL ymds2ju(annee_ref, 1, idayref, hour - pdtphys, zjulian)
317      PRINT *, 'phys_output_mod: annee,iday,hour,zjulian=', annee_ref, idayref, hour, zjulian
318    ELSE
319      CALL ymds2ju(annee_ref, 1, idayref, 0.0, zjulian)
320      CALL ymds2ju(annee_ref, 1, day_ini, start_time * rday, zjulian_start)
321    ENDIF
322
323    IF (using_xios) THEN
324      ! ug R\'eglage du calendrier xios
325      !Temps julian => an, mois, jour, heure
326      CALL ju2ymds(zjulian, x_an, x_mois, x_jour, x_heure)
327      CALL ju2ymds(zjulian_start, ini_an, ini_mois, ini_jour, ini_heure)
328      CALL wxios_set_cal(dtime, calend, x_an, x_mois, x_jour, x_heure, ini_an, &
329              ini_mois, ini_jour, ini_heure)
330    ENDIF
331
332    !!!!!!!!!!!!!!!!!!!!!!! Boucle sur les fichiers !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
333    ! Appel de histbeg et histvert pour creer le fichier et les niveaux verticaux !!
334    ! Appel des histbeg pour definir les variables (nom, moy ou inst, freq de sortie ..
335    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
336
337    zdtime_moy = dtime         ! Frequence ou l on moyenne
338
339    ecrit_files(7) = ecrit_files(1)
340    ecrit_files(8) = ecrit_files(2)
341    ecrit_files(9) = ecrit_files(3)
342
343    DO iff = 1, nfiles
344
345      ! Calculate ecrit_files for all files
346      IF (chtimestep(iff)=='Default') THEN
347        ! Par defaut ecrit_files = (ecrit_mensuel ecrit_jour ecrit_hf
348        ! ...)*86400.
349        ecrit_files(iff) = ecrit_files(iff) * 86400.
350      ELSE IF (chtimestep(iff)=='-1') THEN
351        PRINT*, 'ecrit_files(', iff, ') < 0 so IOIPSL work on different'
352        PRINT*, 'months length'
353        ecrit_files(iff) = -1.
354      ELSE
355        CALL convers_timesteps(chtimestep(iff), dtime, ecrit_files(iff))
356      ENDIF
357
358      WRITE(lunout, *)'ecrit_files(', iff, ')= ', ecrit_files(iff)
359      zoutm(iff) = ecrit_files(iff) ! Frequence ou l on ecrit en seconde
360
361      IF (using_xios) THEN
362        !!! Ouverture de chaque fichier XIOS !!!!!!!!!!!
363        IF (.NOT. ok_all_xml) THEN
364          IF (prt_level >= 10) THEN
365            PRINT*, 'phys_output_open: CALL wxios_add_file with phys_out_filenames(iff)=', trim(phys_out_filenames(iff))
366          ENDIF
367          CALL wxios_add_file(phys_out_filenames(iff), chtimestep(iff), lev_files(iff))
368        ENDIF
369
370        !!! Declaration des axes verticaux de chaque fichier:
371        IF (prt_level >= 10) THEN
372          PRINT*, 'phys_output_open: Declare vertical axes for each file'
373        ENDIF
374
375        IF (iff<=6.OR.iff==10) THEN
376          CALL wxios_add_vaxis("presnivs", &
377                  levmax(iff) - levmin(iff) + 1, presnivs(levmin(iff):levmax(iff)))
378          CALL wxios_add_vaxis("presinter", &
379                  klev + 1, presinter(1:klev + 1))
380          CALL wxios_add_vaxis("Ahyb", &
381                  levmax(iff) - levmin(iff) + 1, aps(levmin(iff):levmax(iff)), positif = 'down', &
382                  bnds = Ahyb_bounds(levmin(iff):levmax(iff), :))
383          CALL wxios_add_vaxis("Bhyb", &
384                  levmax(iff) - levmin(iff) + 1, bps(levmin(iff):levmax(iff)), positif = 'down', &
385                  bnds = Bhyb_bounds(levmin(iff):levmax(iff), :))
386          CALL wxios_add_vaxis("klev", levmax(iff) - levmin(iff) + 1, &
387                  lev_index(levmin(iff):levmax(iff)))
388          CALL wxios_add_vaxis("klevp1", klev + 1, &
389                  lev_index(1:klev + 1))
390          CALL wxios_add_vaxis("bnds", 2, (/1., 2./))
391
392          CALL wxios_add_vaxis("Alt", &
393                  levmax(iff) - levmin(iff) + 1, pseudoalt)
394
395          ! wl1_sun/wl2_sun: minimum/maximum bound of wavelength (in um)
396          SELECT CASE(NSW)
397          CASE(6)
398            wl1_sun(1:6) = [0.180, 0.250, 0.440, 0.690, 1.190, 2.380]
399            wl2_sun(1:6) = [0.250, 0.440, 0.690, 1.190, 2.380, 4.000]
400          CASE(2)
401            wl1_sun(1:2) = [0.250, 0.690]
402            wl2_sun(1:2) = [0.690, 4.000]
403          END SELECT
404
405          DO ISW = 1, NSW
406            wn1_sun(ISW) = 1.e+6 / wl1_sun(ISW)
407            wn2_sun(ISW) = 1.e+6 / wl2_sun(ISW)
408            spbnds_sun(ISW, 1) = wn2_sun(ISW)
409            spbnds_sun(ISW, 2) = wn1_sun(ISW)
410            spectband(ISW) = (wn1_sun(ISW) + wn2_sun(ISW)) / 2
411          ENDDO
412
413          !!! ajout axe vertical spectband : solar band number
414          CALL wxios_add_vaxis("spectband", NSW, spectband, positif = 'down')
415        ELSE
416          ! NMC files
417          CALL wxios_add_vaxis("plev", &
418                  levmax(iff) - levmin(iff) + 1, rlevSTD(levmin(iff):levmax(iff)))
419        ENDIF
420      ENDIF !using_xios
421
422      IF (clef_files(iff)) THEN
423        !!!!!!!!!!!!!!!!! Traitement dans le cas ou l'on veut stocker sur un domaine limite !!
424        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
425        IF (phys_out_regfkey(iff)) THEN
426          imin_ins = 1
427          imax_ins = nbp_lon
428          jmin_ins = 1
429          jmax_ins = jjmp1
430
431          ! correction abderr
432          DO i = 1, nbp_lon
433            WRITE(lunout, *)'io_lon(i)=', io_lon(i)
434            IF (io_lon(i)<=phys_out_lonmin(iff)) imin_ins = i
435            IF (io_lon(i)<=phys_out_lonmax(iff)) imax_ins = i + 1
436          ENDDO
437
438          DO j = 1, jjmp1
439            WRITE(lunout, *)'io_lat(j)=', io_lat(j)
440            IF (io_lat(j)>=phys_out_latmin(iff)) jmax_ins = j + 1
441            IF (io_lat(j)>=phys_out_latmax(iff)) jmin_ins = j
442          ENDDO
443
444          WRITE(lunout, *)'On stoke le fichier histoire numero ', iff, ' sur ', &
445                  imin_ins, imax_ins, jmin_ins, jmax_ins
446          WRITE(lunout, *)'longitudes : ', &
447                  io_lon(imin_ins), io_lon(imax_ins), &
448                  'latitudes : ', &
449                  io_lat(jmax_ins), io_lat(jmin_ins)
450
451          CALL histbeg(phys_out_filenames(iff), nbp_lon, io_lon, jjmp1, io_lat, &
452                  imin_ins, imax_ins - imin_ins + 1, &
453                  jmin_ins, jmax_ins - jmin_ins + 1, &
454                  itau_phy, zjulian, dtime, nhorim(iff), nid_files(iff))
455          !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
456          !IM fichiers stations
457        ELSE IF (clef_stations(iff)) THEN
458
459          IF (prt_level >= 10) THEN
460            WRITE(lunout, *)'phys_output_open: iff=', iff, '  phys_out_filenames(iff)=', phys_out_filenames(iff)
461          ENDIF
462
463          CALL histbeg_phy_all(rlon, rlat, pim, tabij, ipt, jpt, plon, plat, plon_bounds, plat_bounds, &
464                  phys_out_filenames(iff), &
465                  itau_phy, zjulian, dtime, nhorim(iff), nid_files(iff))
466        ELSE
467
468          IF (prt_level >= 10) THEN
469            WRITE(lunout, *)'phys_output_open: iff=', iff, '  phys_out_filenames(iff)=', phys_out_filenames(iff)
470          ENDIF
471
472          CALL histbeg_phy_all(phys_out_filenames(iff), itau_phy, zjulian, &
473                  dtime, nhorim(iff), nid_files(iff))
474        ENDIF
475
476#ifndef CPP_IOIPSL_NO_OUTPUT
477        IF (iff<=6.OR.iff==10) THEN
478          CALL histvert(nid_files(iff), "presnivs", "Vertical levels", "Pa", &
479                  levmax(iff) - levmin(iff) + 1, &
480                  presnivs(levmin(iff):levmax(iff)), nvertm(iff), "down")
481          !!!! Composantes de la coordonnee sigma-hybride
482          CALL histvert(nid_files(iff), "Ahyb", "Ahyb comp of Hyb Cord ", "Pa", &
483                  levmax(iff) - levmin(iff) + 1, aps, nvertap(iff))
484
485          CALL histvert(nid_files(iff), "Bhyb", "Bhyb comp of Hyb Cord", " ", &
486                  levmax(iff) - levmin(iff) + 1, bps, nvertbp(iff))
487
488          CALL histvert(nid_files(iff), "Alt", "Height approx for scale heigh of 8km at levels", "Km", &
489                  levmax(iff) - levmin(iff) + 1, pseudoalt, nvertAlt(iff))
490
491        ELSE
492          ! NMC files
493          CALL histvert(nid_files(iff), "plev", "pressure", "Pa", &
494                  levmax(iff) - levmin(iff) + 1, &
495                  rlevSTD(levmin(iff):levmax(iff)), nvertm(iff), "down")
496        ENDIF
497#endif
498
499      ENDIF ! clef_files
500
501      itr = 0; itrb = 0
502      DO iq = 1, nqtot
503        IF(.NOT.(tracers(iq)%isAdvected .AND. tracers(iq)%isInPhysics)) CYCLE
504        itr = itr + 1
505        dn = 'd' // TRIM(tracers(iq)%name) // '_'
506
507        flag = [1, 5, 5, 5, 10, 10, 11, 11, 11, 11]
508        lnam = 'Tracer ' // TRIM(tracers(iq)%longName)
509        tnam = TRIM(tracers(iq)%name);  o_trac          (itr) = ctrl_out(flag, tnam, lnam, "-", [('', i = 1, nfiles)])
510
511        flag = [4, 7, 7, 7, 10, 10, 11, 11, 11, 11]
512        lnam = 'Tendance tracer ' // TRIM(tracers(iq)%longName)
513        tnam = TRIM(dn) // 'vdf';         o_dtr_vdf       (itr) = ctrl_out(flag, tnam, lnam, "-", [('', i = 1, nfiles)])
514
515        flag = [5, 7, 7, 7, 10, 10, 11, 11, 11, 11]
516        tnam = TRIM(dn) // 'the';         o_dtr_the       (itr) = ctrl_out(flag, tnam, lnam, "-", [('', i = 1, nfiles)])
517        tnam = TRIM(dn) // 'con';         o_dtr_con       (itr) = ctrl_out(flag, tnam, lnam, "-", [('', i = 1, nfiles)])
518
519        flag = [7, 7, 7, 7, 10, 10, 11, 11, 11, 11]
520        tnam = TRIM(dn) // 'lessi_impa';  o_dtr_lessi_impa(itr) = ctrl_out(flag, tnam, lnam, "-", [('', i = 1, nfiles)])
521        tnam = TRIM(dn) // 'lessi_nucl';  o_dtr_lessi_nucl(itr) = ctrl_out(flag, tnam, lnam, "-", [('', i = 1, nfiles)])
522        tnam = TRIM(dn) // 'insc';        o_dtr_insc      (itr) = ctrl_out(flag, tnam, lnam, "-", [('', i = 1, nfiles)])
523        tnam = TRIM(dn) // 'bcscav';      o_dtr_bcscav    (itr) = ctrl_out(flag, tnam, lnam, "-", [('', i = 1, nfiles)])
524        tnam = TRIM(dn) // 'evapls';      o_dtr_evapls    (itr) = ctrl_out(flag, tnam, lnam, "-", [('', i = 1, nfiles)])
525        tnam = TRIM(dn) // 'ls';          o_dtr_ls        (itr) = ctrl_out(flag, tnam, lnam, "-", [('', i = 1, nfiles)])
526        tnam = TRIM(dn) // 'trsp';        o_dtr_trsp      (itr) = ctrl_out(flag, tnam, lnam, "-", [('', i = 1, nfiles)])
527        tnam = TRIM(dn) // 'sscav';       o_dtr_sscav     (itr) = ctrl_out(flag, tnam, lnam, "-", [('', i = 1, nfiles)])
528        tnam = TRIM(dn) // 'sat';         o_dtr_sat       (itr) = ctrl_out(flag, tnam, lnam, "-", [('', i = 1, nfiles)])
529        tnam = TRIM(dn) // 'uscav';       o_dtr_uscav     (itr) = ctrl_out(flag, tnam, lnam, "-", [('', i = 1, nfiles)])
530
531        lnam = 'tracer tendency dry deposition' // TRIM(tracers(iq)%longName)
532        tnam = 'cum' // TRIM(dn) // 'dry';  o_dtr_dry       (itr) = ctrl_out(flag, tnam, lnam, "-", [('', i = 1, nfiles)])
533
534        flag = [1, 4, 10, 10, 10, 10, 11, 11, 11, 11]
535        lnam = 'Cumulated tracer ' // TRIM(tracers(iq)%longName)
536        tnam = 'cum' // TRIM(tracers(iq)%name); o_trac_cum(itr) = ctrl_out(flag, tnam, lnam, "-", [('', i = 1, nfiles)])
537
538        IF (CPPKEY_STRATAER) THEN
539          IF(tracers(iq)%name(1:3)=='BIN') THEN
540            itrb = itrb + 1
541            flag = [11, 11, 11, 11, 11, 11, 11, 11, 11, 1]
542            lnam = 'Dry particle concentration in ' // TRIM(tracers(iq)%longName)
543            tnam = TRIM(tracers(iq)%name) // '_nd_mode';     o_nd_mode       (itrb) = ctrl_out(flag, tnam, lnam, "part/m3", [('', i = 1, nfiles)])
544            lnam = 'Sulfate MMR in ' // TRIM(tracers(iq)%longName)
545            tnam = TRIM(tracers(iq)%name) // '_sulfmmr_mode';o_sulfmmr_mode  (itrb) = ctrl_out(flag, tnam, lnam, "kg(H2SO4)/kg(air)", [('', i = 1, nfiles)])
546          endif
547        END IF
548      ENDDO
549
550    ENDDO !  iff
551
552#ifdef ISO
553    WRITE(*,*) 'phys_output_mid 589'
554    DO ixt=1,ntraciso
555      outiso = TRIM(isoName(ixt))
556      i = INDEX(outiso, '_', .TRUE.)
557      outiso = outiso(1:i-1)//outiso(i+1:LEN_TRIM(outiso))
558
559      flag = [1,  1,  1, 10,  5, 10, 11, 11, 11, 11]; unit = 'kg/(s*m2)'
560      o_xtprecip(ixt)=ctrl_out(flag, 'precip'//TRIM(outiso), 'Precip Totale liq+sol', unit, [('',i=1,nfiles)])
561      o_xtpluc  (ixt)=ctrl_out(flag,   'pluc'//TRIM(outiso),    'Convective Precip.', unit, [('',i=1,nfiles)])
562
563      flag = [1,  1,  1, 10, 10, 10, 11, 11, 11, 11]
564      o_xtplul  (ixt)=ctrl_out(flag,   'plul'//TRIM(outiso),   'Large-scale Precip.', unit, [('',i=1,nfiles)])
565      o_xtevap  (ixt)=ctrl_out(flag,   'evap'//TRIM(outiso),             'Evaporat.', unit, [('',i=1,nfiles)])
566
567      ! ajout Camille 8 mai 2023
568      flag = [1, 6, 10, 10, 10, 10, 11, 11, 11, 11]
569      o_xtevap_srf (ixt,1)=ctrl_out(flag,   'evap_ter'//TRIM(outiso), 'Evap sfc'//clnsurf(1), unit, [('',i=1,nfiles)])
570      o_xtevap_srf (ixt,2)=ctrl_out(flag,   'evap_lic'//TRIM(outiso), 'Evap sfc'//clnsurf(2), unit, [('',i=1,nfiles)])
571      o_xtevap_srf (ixt,3)=ctrl_out(flag,   'evap_oce'//TRIM(outiso), 'Evap sfc'//clnsurf(3), unit, [('',i=1,nfiles)])
572      o_xtevap_srf (ixt,4)=ctrl_out(flag,   'evap_sic'//TRIM(outiso), 'Evap sfc'//clnsurf(4), unit, [('',i=1,nfiles)])
573
574      flag = [2,  3,  4, 10, 10, 10, 11, 11, 11, 11]; unit = 'kg/kg'
575      o_xtovap  (ixt)=ctrl_out(flag,   'ovap'//TRIM(outiso),     'Specific humidity', unit, [('',i=1,nfiles)])
576      o_xtoliq  (ixt)=ctrl_out(flag,   'oliq'//TRIM(outiso),          'Liquid water', unit, [('',i=1,nfiles)])
577      o_xtcond  (ixt)=ctrl_out(flag,  'ocond'//TRIM(outiso),       'Condensed water', unit, [('',i=1,nfiles)])
578
579      flag = [1,  1,  1, 10, 5, 10, 11, 11, 11, 11]; unit = 'kg/m2/s'
580      o_xtrunoff_diag  (ixt)=ctrl_out(flag, 'runoffland'//TRIM(outiso), 'Run-off rate land for bucket', unit, [('',i=1,nfiles)])
581
582      flag = [4, 10, 10, 10, 10, 10, 11, 11, 11, 11]; unit = '(kg/kg)/s'
583      o_dxtdyn  (ixt)=ctrl_out(flag,  'dqdyn'//TRIM(outiso),           'Dynamics dQ', unit, [('',i=1,nfiles)])
584      o_dxtldyn (ixt)=ctrl_out(flag, 'dqldyn'//TRIM(outiso),          'Dynamics dQL', unit, [('',i=1,nfiles)])
585      o_dxtcon  (ixt)=ctrl_out(flag,  'dqcon'//TRIM(outiso),         'Convection dQ', unit, [('',i=1,nfiles)])
586      o_dxteva  (ixt)=ctrl_out(flag,  'dqeva'//TRIM(outiso),      'Reevaporation dQ', unit, [('',i=1,nfiles)])
587      o_dxtlsc  (ixt)=ctrl_out(flag,  'dqlsc'//TRIM(outiso),       'Condensation dQ', unit, [('',i=1,nfiles)])
588      o_dxtajs  (ixt)=ctrl_out(flag,  'dqajs'//TRIM(outiso),        'Dry adjust. dQ', unit, [('',i=1,nfiles)])
589      o_dxtvdf  (ixt)=ctrl_out(flag,  'dqvdf'//TRIM(outiso),     'Boundary-layer dQ', unit, [('',i=1,nfiles)])
590      o_dxtthe  (ixt)=ctrl_out(flag,  'dqthe'//TRIM(outiso),            'Thermal dQ', unit, [('',i=1,nfiles)])
591
592      IF(ok_qch4) o_dxtch4(ixt)=ctrl_out(flag, 'dqch4'//TRIM(outiso), 'H2O due to CH4 oxidation & photolysis', &
593                                                                                      unit, [('',i=1,nfiles)])
594      IF(ixt == iso_HTO) THEN
595      o_dxtprod_nucl(ixt)=ctrl_out(flag, 'dqprodnucl'//TRIM(outiso), 'dHTO/dt due to nuclear production',      &
596                                                                                      unit, [('',i=1,nfiles)])
597      o_dxtcosmo    (ixt)=ctrl_out(flag,    'dqcosmo'//TRIM(outiso), 'dHTO/dt due to cosmogenic production',   &
598                                                                                      unit, [('',i=1,nfiles)])
599      o_dxtdecroiss (ixt)=ctrl_out(flag, 'dqdecroiss'//TRIM(outiso), 'dHTO/dt due to radiative destruction',   &
600                                                                                      unit, [('',i=1,nfiles)])
601      END IF
602    enddo !do ixt=1,niso
603    WRITE(*,*) 'phys_output_mid 596'
604#endif
605
606    ! Updated write frequencies due to phys_out_filetimesteps.
607    ! Write frequencies are now in seconds. 
608    ecrit_mth = ecrit_files(1)
609    ecrit_day = ecrit_files(2)
610    ecrit_hf = ecrit_files(3)
611    ecrit_ins = ecrit_files(4)
612    ecrit_LES = ecrit_files(5)
613    ecrit_ins = ecrit_files(6)
614
615    IF (prt_level >= 10) THEN
616      WRITE(lunout, *)'swaerofree_diag=', swaerofree_diag
617      WRITE(lunout, *)'swaero_diag=', swaero_diag
618      WRITE(lunout, *)'dryaod_diag=', dryaod_diag
619      WRITE(lunout, *)'ok_4xCO2atm=', ok_4xCO2atm
620      WRITE(lunout, *)'phys_output_open: ends here'
621    ENDIF
622
623    !  DO iq=1,nqtot
624    !    IF(.NOT.(tracers(iq)%isAdvected .AND. tracers(iq)%isInPhysics)) CYCLE
625    !    WRITE(*,'(a,i1,a,10i3)')'trac(',iq,')%flag = ',o_trac(iq)%flag
626    !    WRITE(*,'(a,i1,a)')'trac(',iq,')%name = '//TRIM(o_trac(iq)%name)
627    !    WRITE(*,'(a,i1,a)')'trac(',iq,')%description = '//TRIM(o_trac(iq)%description)
628    !  END DO
629
630  END SUBROUTINE phys_output_open
631
632  SUBROUTINE convers_timesteps(str, dtime, timestep)
633
634    USE ioipsl
635    USE phys_cal_mod
636    USE time_phylmdz_mod, ONLY: day_ref, annee_ref
637    USE lmdz_print_control, ONLY: lunout
638    USE lmdz_abort_physic, ONLY: abort_physic
639
640    IMPLICIT NONE
641
642    CHARACTER(LEN = 20) :: str
643    CHARACTER(LEN = 10) :: type
644    INTEGER :: ipos, il
645    REAL :: ttt, xxx, timestep, dayseconde, dtime
646    parameter (dayseconde = 86400.)
647
648    ipos = scan(str, '0123456789.', .TRUE.)
649
650    il = len_trim(str)
651    WRITE(lunout, *) "ipos = ", ipos
652    WRITE(lunout, *) "il = ", il
653    IF (ipos == 0) CALL abort_physic("convers_timesteps", "bad str", 1)
654    read(str(1:ipos), *) ttt
655    WRITE(lunout, *)ttt
656    type = str(ipos + 1:il)
657
658    IF (il == ipos) THEN
659      type = 'day'
660    ENDIF
661
662    IF (type == 'day'.OR.type == 'days'.OR.type == 'jours'.OR.type == 'jour') timestep = ttt * dayseconde
663    IF (type == 'mounths'.OR.type == 'mth'.OR.type == 'mois') THEN
664      WRITE(lunout, *)'annee_ref,day_ref mon_len', annee_ref, day_ref, mth_len
665      timestep = ttt * dayseconde * mth_len
666    ENDIF
667    IF (type == 'hours'.OR.type == 'hr'.OR.type == 'heurs') timestep = ttt * dayseconde / 24.
668    IF (type == 'mn'.OR.type == 'minutes') timestep = ttt * 60.
669    IF (type == 's'.OR.type == 'sec'.OR.type == 'secondes') timestep = ttt
670    IF (type == 'TS') timestep = ttt * dtime
671
672    WRITE(lunout, *)'type =      ', type
673    WRITE(lunout, *)'nb j/h/m =  ', ttt
674    WRITE(lunout, *)'timestep(s)=', timestep
675
676  END SUBROUTINE convers_timesteps
677
678END MODULE phys_output_mod
Note: See TracBrowser for help on using the repository browser.