source: LMDZ6/trunk/libf/phylmdiso/phys_output_mod.F90 @ 4500

Last change on this file since 4500 was 4170, checked in by dcugnet, 2 years ago

The variable "types_trac" is the equivalent of "type_trac" in case multiple sections must be read
and used in "tracer.def" file.
Tests on the "type_trac" were replaced with tests on the vector "types_trac".
Most of the time, there are two components: 'lmdz' and a second one. The later has priority on 'lmdz'
and must be used for the tests. For more components, care must be taken to execute specific parts
of the code on the right tracers ; the tracers(:)%component has been created in that respect.

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