source: LMDZ6/trunk/libf/phylmd/traclmdz_mod.F90 @ 4046

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

First commit for new tracers.

  • parser routines readTracFiles, strings_mod and tracer_types added in misc using revision 4 of https://svn.lmd.jussieu.fr/tracers-parser
  • tested in sequential and parallel mode using ioipsl.
  • for now, only two fields of "tracers(:)" derived type vector are used: "name" and "longName".
  • 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
File size: 22.3 KB
Line 
1!$Id $
2!
3MODULE traclmdz_mod
4
5!
6! In this module all tracers specific to LMDZ are treated. This module is used
7! only if running without any other chemestry model as INCA or REPROBUS. 
8!
9  IMPLICIT NONE
10
11  REAL,DIMENSION(:,:),ALLOCATABLE,SAVE :: masktr   ! Masque reservoir de sol traceur
12!$OMP THREADPRIVATE(masktr)                        ! Masque de l'echange avec la surface (1 = reservoir) ou (possible >= 1 )
13  REAL,DIMENSION(:,:),ALLOCATABLE,SAVE :: fshtr    ! Flux surfacique dans le reservoir de sol
14!$OMP THREADPRIVATE(fshtr)
15  REAL,DIMENSION(:),ALLOCATABLE,SAVE   :: hsoltr   ! Epaisseur equivalente du reservoir de sol
16!$OMP THREADPRIVATE(hsoltr)
17!
18!Radioelements:
19!--------------
20!
21  REAL,DIMENSION(:),ALLOCATABLE,SAVE   :: tautr    ! Constante de decroissance radioactive
22!$OMP THREADPRIVATE(tautr)
23  REAL,DIMENSION(:),ALLOCATABLE,SAVE   :: vdeptr   ! Vitesse de depot sec dans la couche Brownienne
24!$OMP THREADPRIVATE(vdeptr)
25  REAL,DIMENSION(:),ALLOCATABLE,SAVE   :: scavtr   ! Coefficient de lessivage
26!$OMP THREADPRIVATE(scavtr)
27  REAL,DIMENSION(:,:),ALLOCATABLE,SAVE :: srcbe    ! Production du beryllium7 dans l atmosphere (U/s/kgA)
28!$OMP THREADPRIVATE(srcbe)
29
30  LOGICAL,DIMENSION(:),ALLOCATABLE,SAVE :: radio    ! radio(it)   = true  => decroisssance radioactive
31!$OMP THREADPRIVATE(radio) 
32
33  REAL,DIMENSION(:,:),ALLOCATABLE,SAVE :: trs     ! Conc. radon ds le sol
34!$OMP THREADPRIVATE(trs)
35
36  INTEGER,SAVE :: id_aga      ! Identification number for tracer : Age of stratospheric air
37!$OMP THREADPRIVATE(id_aga)
38  INTEGER,SAVE :: lev_1p5km   ! Approximative vertical layer number at 1.5km above surface, used for calculation of the age of air. The result shouldn't be that sensible to the exactness of this value as long as it is in the lower troposphere.
39!$OMP THREADPRIVATE(lev_1p5km)
40
41  INTEGER,SAVE :: id_rn, id_pb ! Identification number for tracer : radon (Rn222), lead (Pb210)
42!$OMP THREADPRIVATE(id_rn, id_pb)
43
44  INTEGER,SAVE :: id_be       ! Activation et position du traceur Be7 [ id_be=0 -> desactive ]
45!$OMP THREADPRIVATE(id_be)
46
47  INTEGER,SAVE :: id_pcsat, id_pcocsat, id_pcq ! traceurs pseudo-vapeur CL qsat, qsat_oc, q
48!$OMP THREADPRIVATE(id_pcsat, id_pcocsat, id_pcq)
49  INTEGER,SAVE :: id_pcs0, id_pcos0, id_pcq0   ! traceurs pseudo-vapeur CL qsat, qsat_oc, q
50!                                              ! qui ne sont pas transportes par la convection
51!$OMP THREADPRIVATE(id_pcs0, id_pcos0, id_pcq0)
52
53  INTEGER, SAVE:: id_o3
54!$OMP THREADPRIVATE(id_o3)
55! index of ozone tracer with Cariolle parameterization
56! 0 means no ozone tracer
57
58  LOGICAL,SAVE :: rnpb=.FALSE. ! Presence du couple Rn222, Pb210
59!$OMP THREADPRIVATE(rnpb)
60
61
62CONTAINS
63
64  SUBROUTINE traclmdz_from_restart(trs_in)
65    ! This subroutine initialize the module saved variable trs with values from restart file (startphy.nc).
66    ! This subroutine is called from phyetat0 after the field trs_in has been read.
67   
68    USE dimphy
69    USE infotrac_phy
70   
71    ! Input argument
72    REAL,DIMENSION(klon,nbtr), INTENT(IN) :: trs_in
73   
74    ! Local variables
75    INTEGER :: ierr
76   
77    ! Allocate restart variables trs
78    ALLOCATE( trs(klon,nbtr), stat=ierr)
79    IF (ierr /= 0) CALL abort_physic('traclmdz_from_restart', 'pb in allocation 1',1)
80   
81    ! Initialize trs with values read from restart file
82    trs(:,:) = trs_in(:,:)
83   
84  END SUBROUTINE traclmdz_from_restart
85
86
87  SUBROUTINE traclmdz_init(pctsrf, xlat, xlon, ftsol, tr_seri, t_seri, pplay, sh, pdtphys, aerosol, lessivage)
88    ! This subroutine allocates and initialize module variables and control variables.
89    ! Initialization of the tracers should be done here only for those not found in the restart file.
90    USE dimphy
91    USE infotrac_phy
92    USE regr_pr_comb_coefoz_m, ONLY: alloc_coefoz
93    USE press_coefoz_m, ONLY: press_coefoz
94    USE mod_grid_phy_lmdz
95    USE mod_phys_lmdz_para
96    USE indice_sol_mod
97    USE print_control_mod, ONLY: lunout
98    USE strings_mod, ONLY: strLower
99
100! Input variables
101    REAL,DIMENSION(klon,nbsrf),INTENT(IN)     :: pctsrf ! Pourcentage de sol f(nature du sol)
102    REAL,DIMENSION(klon),INTENT(IN)           :: xlat   ! latitudes en degres pour chaque point
103    REAL,DIMENSION(klon),INTENT(IN)           :: xlon   ! longitudes en degres pour chaque point
104    REAL,DIMENSION(klon,nbsrf),INTENT(IN)     :: ftsol  ! Temperature du sol (surf)(Kelvin)
105    REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT) :: tr_seri! Concentration Traceur [U/KgA] 
106    REAL,DIMENSION(klon,klev),INTENT(IN)   :: t_seri  ! Temperature
107    REAL,DIMENSION(klon,klev),INTENT(IN)   :: pplay   ! pression pour le mileu de chaque couche (en Pa)
108    REAL,DIMENSION(klon,klev),INTENT(IN)   :: sh      ! humidite specifique
109    REAL,INTENT(IN)                        :: pdtphys ! Pas d'integration pour la physique (seconde) 
110
111! Output variables
112    LOGICAL,DIMENSION(nbtr), INTENT(OUT) :: aerosol
113    LOGICAL,INTENT(OUT)                  :: lessivage
114       
115! Local variables   
116    INTEGER :: ierr, it, iiq, i, k
117    REAL, DIMENSION(klon_glo,klev) :: varglo ! variable temporaire sur la grille global   
118    REAL, DIMENSION(klev)          :: mintmp, maxtmp
119    LOGICAL                        :: zero
120! RomP >>> profil initial Be7
121      integer ilesfil
122      parameter (ilesfil=1)
123      integer  irr,kradio
124      real     beryllium(klon,klev)
125! profil initial Pb210
126      integer ilesfil2
127      parameter (ilesfil2=1)
128      integer  irr2,kradio2
129      real     plomb(klon,klev)
130!! RomP <<<
131! --------------------------------------------
132! Allocation
133! --------------------------------------------
134    ALLOCATE( scavtr(nbtr), stat=ierr)
135    IF (ierr /= 0) CALL abort_physic('traclmdz_init', 'pb in allocation 9',1)
136    scavtr(:)=1.
137   
138    ALLOCATE( radio(nbtr), stat=ierr)
139    IF (ierr /= 0) CALL abort_physic('traclmdz_init', 'pb in allocation 11',1)
140    radio(:) = .false.    ! Par defaut pas decroissance radioactive
141   
142    ALLOCATE( masktr(klon,nbtr), stat=ierr)
143    IF (ierr /= 0) CALL abort_physic('traclmdz_init', 'pb in allocation 2',1)
144   
145    ALLOCATE( fshtr(klon,nbtr), stat=ierr)
146    IF (ierr /= 0) CALL abort_physic('traclmdz_init', 'pb in allocation 3',1)
147   
148    ALLOCATE( hsoltr(nbtr), stat=ierr)
149    IF (ierr /= 0) CALL abort_physic('traclmdz_init', 'pb in allocation 4',1)
150   
151    ALLOCATE( tautr(nbtr), stat=ierr)
152    IF (ierr /= 0) CALL abort_physic('traclmdz_init', 'pb in allocation 5',1)
153    tautr(:)  = 0.
154   
155    ALLOCATE( vdeptr(nbtr), stat=ierr)
156    IF (ierr /= 0) CALL abort_physic('traclmdz_init', 'pb in allocation 6',1)
157    vdeptr(:) = 0.
158
159
160    lessivage  = .TRUE.
161!!jyg(20130206) : le choix d activation du lessivage est fait dans phytrac avec iflag_lscav
162!!    call getin('lessivage',lessivage)
163!!    if(lessivage) then
164!!     print*,'lessivage lsc ON'
165!!    else
166!!     print*,'lessivage lsc OFF'
167!!    endif
168    aerosol(:) = .FALSE.  ! Tous les traceurs sont des gaz par defaut
169   
170!
171! Recherche des traceurs connus : Be7, O3, CO2,...
172! --------------------------------------------
173    id_rn=0; id_pb=0; id_aga=0; id_be=0; id_o3=0
174    id_pcsat=0; id_pcocsat=0; id_pcq=0; id_pcs0=0; id_pcos0=0; id_pcq0=0
175    DO it=1,nbtr
176!!       iiq=niadv(it+2)                                                            ! jyg
177       iiq=niadv(it+nqo)                                                            ! jyg
178print*,'###'//TRIM(tracers(iiq)%name)//'###'
179print*,'###'//TRIM(strLower(tracers(iiq)%name))//'###'
180       SELECT CASE(strLower(tracers(iiq)%name))
181         CASE("rn");      id_rn     = it ! radon
182         CASE("pb");      id_pb     = it ! plomb
183         CASE("aga");     id_aga    = it ! Age of stratospheric air
184         CASE("be","be7");id_be     = it ! Recherche du Beryllium 7
185         CASE("o3");      id_o3     = it ! Recherche de l'ozone
186         CASE("pcsat");   id_pcsat  = it
187         CASE("pcocsat"); id_pcocsat= it
188         CASE("pcq");     id_pcq    = it
189         CASE("pcs0");    id_pcs0   = it
190         CASE("pcos0");   id_pcos0  = it
191         CASE("pcq0");    id_pcq0   = it
192         CASE DEFAULT
193           WRITE(lunout,*) 'This is an unknown tracer in LMDZ : ', trim(tracers(iiq)%name)
194       END SELECT
195
196       SELECT CASE(strLower(tracers(iiq)%name))
197         CASE("pb")                        !--- RomP >>> profil initial de PB210
198           OPEN(ilesfil2,file='prof.pb210',status='old',iostat=irr2)
199           IF(irr2 == 0) THEN
200             READ(ilesfil2,*) kradio2
201             WRITE(lunout,*)'number of levels for pb210 profile ',kradio2
202             DO k=kradio2,1,-1; READ (ilesfil2,*) plomb(:,k); END DO
203             CLOSE(ilesfil2)
204             tr_seri(:,:,id_pb) = plomb(:,:)
205           ELSE
206             WRITE(lunout,*)'Prof. Pb210 does not exist: use restart values'
207           END IF
208         CASE("aga")
209           radio(id_aga) = .FALSE.
210           aerosol(id_aga) = .FALSE.
211           pbl_flg(id_aga) = 0
212           ! Find the first model layer above 1.5km from the surface
213           IF (klev>=30) THEN
214              lev_1p5km=6                  !--- NB: This value is for klev=39
215           ELSE IF (klev>=10) THEN
216              lev_1p5km=5                  !--- NB: This value is for klev=19
217           ELSE
218              lev_1p5km=klev/2
219           END IF
220         CASE("be","be7")
221           ALLOCATE( srcbe(klon,klev) )
222           radio(id_be) = .TRUE.
223           aerosol(id_be) = .TRUE.         !--- Le Be est un aerosol
224           CALL init_be(pctsrf,pplay,masktr(:,id_be),tautr(id_be),vdeptr(id_be),scavtr(id_be),srcbe)
225           WRITE(lunout,*) 'Initialisation srcBe: OK'
226                                           !--- RomP >>> profil initial de Be7
227           OPEN(ilesfil,file='prof.be7',status='old',iostat=irr)
228           IF(irr == 0) THEN
229             READ(ilesfil,*) kradio
230             WRITE(lunout,*)'number of levels for Be7 profile ',kradio
231             DO k=kradio,1,-1; READ(ilesfil,*) beryllium(:,k); END DO
232             CLOSE(ilesfil)
233             tr_seri(:,:,id_be)=beryllium(:,:)
234           ELSE
235             WRITE(lunout,*)'Prof. Be7 does not exist: use restart values'
236           END IF
237         CASE("o3")                         !--- Parametrisation par la chimie de Cariolle
238           CALL alloc_coefoz                !--- Allocate ozone coefficients
239           CALL press_coefoz                !--- Read input pressure levels
240         CASE("pcs0","pcos0","pcq0")
241           conv_flg(it)=0                   !--- No transport by convection for this tracer
242       END SELECT
243    END DO
244
245!
246! Valeurs specifiques pour les traceurs Rn222 et Pb210
247! ----------------------------------------------
248    IF ( id_rn/=0 .AND. id_pb/=0 ) THEN
249       rnpb = .TRUE.
250       radio(id_rn)= .TRUE.
251       radio(id_pb)= .TRUE.
252       pbl_flg(id_rn) = 0 ! au lieu de clsol=true ! CL au sol calcule
253       pbl_flg(id_pb) = 0 ! au lieu de clsol=true
254       aerosol(id_rn) = .FALSE.
255       aerosol(id_pb) = .TRUE. ! le Pb est un aerosol
256       
257       CALL initrrnpb (ftsol,pctsrf,masktr,fshtr,hsoltr,tautr,vdeptr,scavtr)
258    END IF
259
260!
261! Check if all tracers have restart values
262! ----------------------------------------------
263    DO it=1,nbtr
264!!       iiq=niadv(it+2)                                                            ! jyg
265       iiq=niadv(it+nqo)                                                            ! jyg
266       ! Test if tracer is zero everywhere.
267       ! Done by master process MPI and master thread OpenMP
268       CALL gather(tr_seri(:,:,it),varglo)
269       IF (is_mpi_root .AND. is_omp_root) THEN
270          mintmp=MINVAL(varglo,dim=1)
271          maxtmp=MAXVAL(varglo,dim=1)
272          IF (MINVAL(mintmp,dim=1)==0. .AND. MAXVAL(maxtmp,dim=1)==0.) THEN
273             ! Tracer is zero everywhere
274             zero=.TRUE.
275          ELSE
276             zero=.FALSE.
277          END IF
278       END IF
279
280       ! Distribute variable at all processes
281       CALL bcast(zero)
282
283       ! Initalize tracer that was not found in restart file.
284       IF (zero) THEN
285          ! The tracer was not found in restart file or it was equal zero everywhere.
286          WRITE(lunout,*) "The tracer ",trim(tracers(iiq)%name)," will be initialized"
287          IF (it==id_pcsat .OR. it==id_pcq .OR. &
288               it==id_pcs0 .OR. it==id_pcq0) THEN
289             tr_seri(:,:,it) = 100.
290          ELSE IF (it==id_pcocsat .OR. it==id_pcos0) THEN
291             DO i = 1, klon
292                IF ( pctsrf (i, is_oce) == 0. ) THEN
293                   tr_seri(i,:,it) = 0.
294                ELSE
295                   tr_seri(i,:,it) = 100.
296                END IF
297             END DO
298          ELSE
299             ! No specific initialization exist for this tracer
300             tr_seri(:,:,it) = 0.
301          END IF
302       END IF
303    END DO
304
305  END SUBROUTINE traclmdz_init
306
307  SUBROUTINE traclmdz(nstep, julien, gmtime, pdtphys, t_seri, paprs, pplay, &
308       cdragh,  coefh, yu1, yv1, ftsol, pctsrf, xlat, xlon, couchelimite, sh, &
309       rh, pphi, ustar, wstar, ale_bl, ale_wake,  zu10m, zv10m, &
310       tr_seri, source, d_tr_cl,d_tr_dec, zmasse)               !RomP
311   
312    USE dimphy
313    USE infotrac_phy
314    USE regr_pr_comb_coefoz_m, ONLY: regr_pr_comb_coefoz
315    USE o3_chem_m, ONLY: o3_chem
316    USE indice_sol_mod
317
318    INCLUDE "YOMCST.h"
319
320!==========================================================================
321!                   -- DESCRIPTION DES ARGUMENTS --
322!==========================================================================
323
324! Input arguments
325!
326!Configuration grille,temps:
327    INTEGER,INTENT(IN) :: nstep      ! nombre d'appels de la physiq
328    INTEGER,INTENT(IN) :: julien     ! Jour julien
329    REAL,INTENT(IN)    :: gmtime
330    REAL,INTENT(IN)    :: pdtphys    ! Pas d'integration pour la physique (seconde) 
331    REAL,DIMENSION(klon),INTENT(IN) :: xlat    ! latitudes pour chaque point
332    REAL, INTENT(IN):: xlon(:) ! dim(klon) longitude
333
334!
335!Physique:
336!--------
337    REAL,DIMENSION(klon,klev),INTENT(IN)   :: t_seri  ! Temperature
338    REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs   ! pression pour chaque inter-couche (en Pa)
339    REAL,DIMENSION(klon,klev),INTENT(IN)   :: pplay   ! pression pour le mileu de chaque couche (en Pa)
340    REAL,intent(in):: zmasse (:, :)   ! dim(klon,klev) density of air, in kg/m2
341
342
343!Couche limite:
344!--------------
345!
346    REAL,DIMENSION(klon),INTENT(IN)      :: cdragh     ! coeff drag pour T et Q
347    REAL,DIMENSION(klon,klev),INTENT(IN) :: coefh      ! coeff melange CL (m**2/s)
348    REAL,DIMENSION(klon),INTENT(IN)      :: yu1        ! vents au premier niveau
349    REAL,DIMENSION(klon),INTENT(IN)      :: yv1        ! vents au premier niveau
350    LOGICAL,INTENT(IN)                   :: couchelimite
351    REAL,DIMENSION(klon,klev),INTENT(IN) :: sh         ! humidite specifique
352    REAL,DIMENSION(klon,klev),INTENT(IN) :: rh      ! Humidite relative
353    REAL,DIMENSION(klon,klev),INTENT(IN) :: pphi    ! geopotentie
354    REAL,DIMENSION(klon),INTENT(IN)      :: ustar   ! ustar (m/s)
355    REAL,DIMENSION(klon),INTENT(IN)      :: wstar,ale_bl,ale_wake   ! wstar (m/s) and Avail. Lifti. Energ.
356    REAL,DIMENSION(klon),INTENT(IN)      :: zu10m   ! vent zonal 10m (m/s)
357    REAL,DIMENSION(klon),INTENT(IN)      :: zv10m   ! vent zonal 10m (m/s)
358
359! Arguments necessaires pour les sources et puits de traceur:
360    REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: ftsol  ! Temperature du sol (surf)(Kelvin)
361    REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: pctsrf ! Pourcentage de sol f(nature du sol)
362
363! InOutput argument
364    REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT)  :: tr_seri ! Concentration Traceur [U/KgA] 
365
366! Output argument
367    REAL,DIMENSION(klon,nbtr), INTENT(OUT)        :: source  ! a voir lorsque le flux de surface est prescrit
368    REAL,DIMENSION(klon,klev,nbtr), INTENT(OUT)   :: d_tr_cl ! Td couche limite/traceur
369
370!=======================================================================================
371!                        -- VARIABLES LOCALES TRACEURS --
372!=======================================================================================
373
374    INTEGER :: i, k, it
375    INTEGER :: lmt_pas ! number of time steps of "physics" per day
376
377    REAL,DIMENSION(klon)           :: d_trs    ! Td dans le reservoir
378    REAL,DIMENSION(klon,klev)      :: qsat     ! pression de la vapeur a saturation
379    REAL,DIMENSION(klon,klev,nbtr) :: d_tr_dec ! Td radioactive
380    REAL                           :: zrho     ! Masse Volumique de l'air KgA/m3
381    REAL                           :: amn, amx
382!
383!=================================================================
384!  Ajout de la production en  Be7 (Beryllium) srcbe U/s/kgA
385!=================================================================
386!
387    IF ( id_be /= 0 ) THEN
388       DO k = 1, klev
389          DO i = 1, klon
390             tr_seri(i,k,id_be) = tr_seri(i,k,id_be)+srcbe(i,k)*pdtphys
391          END DO
392       END DO
393       WRITE(*,*) 'Ajout srcBe dans tr_seri: OK'
394    END IF
395 
396
397!=================================================================
398! Update pseudo-vapor tracers
399!=================================================================
400
401    CALL q_sat(klon*klev,t_seri,pplay,qsat)
402
403    IF ( id_pcsat /= 0 ) THEN
404     DO k = 1, klev
405      DO i = 1, klon
406         IF ( pplay(i,k).GE.85000.) THEN
407            tr_seri(i,k,id_pcsat) = qsat(i,k)
408         ELSE
409            tr_seri(i,k,id_pcsat) = MIN (qsat(i,k), tr_seri(i,k,id_pcsat))           
410         END IF
411      END DO
412     END DO
413    END IF
414
415    IF ( id_pcocsat /= 0 ) THEN
416     DO k = 1, klev
417      DO i = 1, klon
418         IF ( pplay(i,k).GE.85000.) THEN
419            IF ( pctsrf (i, is_oce) > 0. ) THEN
420               tr_seri(i,k,id_pcocsat) = qsat(i,k)
421            ELSE
422               tr_seri(i,k,id_pcocsat) = 0.
423          END IF
424       ELSE
425          tr_seri(i,k,id_pcocsat) = MIN (qsat(i,k), tr_seri(i,k,id_pcocsat))
426       END IF
427      END DO
428     END DO
429    END IF
430
431    IF ( id_pcq /= 0 ) THEN
432     DO k = 1, klev
433      DO i = 1, klon
434         IF ( pplay(i,k).GE.85000.) THEN
435            tr_seri(i,k,id_pcq) = sh(i,k)
436         ELSE
437            tr_seri(i,k,id_pcq) = MIN (qsat(i,k), tr_seri(i,k,id_pcq))   
438         END IF
439      END DO
440     END DO
441    END IF
442
443
444    IF ( id_pcs0 /= 0 ) THEN
445     DO k = 1, klev
446      DO i = 1, klon
447         IF ( pplay(i,k).GE.85000.) THEN
448            tr_seri(i,k,id_pcs0) = qsat(i,k)
449         ELSE
450            tr_seri(i,k,id_pcs0) = MIN (qsat(i,k), tr_seri(i,k,id_pcs0))   
451         END IF
452      END DO
453     END DO
454    END IF
455
456
457    IF ( id_pcos0 /= 0 ) THEN
458     DO k = 1, klev
459      DO i = 1, klon
460         IF ( pplay(i,k).GE.85000.) THEN
461            IF ( pctsrf (i, is_oce) > 0. ) THEN
462               tr_seri(i,k,id_pcos0) = qsat(i,k)
463            ELSE
464               tr_seri(i,k,id_pcos0) = 0.
465            END IF
466         ELSE
467            tr_seri(i,k,id_pcos0) = MIN (qsat(i,k), tr_seri(i,k,id_pcos0))
468         END IF
469      END DO
470     END DO
471    END IF
472
473
474    IF ( id_pcq0 /= 0 ) THEN
475     DO k = 1, klev
476      DO i = 1, klon
477         IF ( pplay(i,k).GE.85000.) THEN
478            tr_seri(i,k,id_pcq0) = sh(i,k)
479         ELSE
480            tr_seri(i,k,id_pcq0) = MIN (qsat(i,k), tr_seri(i,k,id_pcq0))
481         END IF
482      END DO
483     END DO
484    END IF
485
486!=================================================================
487! Update tracer : Age of stratospheric air
488!=================================================================
489    IF (id_aga/=0) THEN
490       
491       ! Bottom layers
492       DO k = 1, lev_1p5km
493          tr_seri(:,k,id_aga) = 0.0
494       END DO
495       
496       ! Layers above 1.5km
497       DO k = lev_1p5km+1,klev-1
498          tr_seri(:,k,id_aga) = tr_seri(:,k,id_aga) + pdtphys
499       END DO
500       
501       ! Top layer
502       tr_seri(:,klev,id_aga) = tr_seri(:,klev-1,id_aga)
503       
504    END IF
505
506!======================================================================
507!     -- Calcul de l'effet de la couche limite --
508!======================================================================
509
510    IF (couchelimite) THEN             
511       source(:,:) = 0.0
512
513       IF (id_be /=0) THEN
514          DO i=1, klon
515             zrho = pplay(i,1)/t_seri(i,1)/RD
516             source(i,id_be) = - vdeptr(id_be)*tr_seri(i,1,id_be)*zrho
517          END DO
518       END IF
519
520    END IF
521   
522    DO it=1, nbtr
523       IF (couchelimite .AND. pbl_flg(it) == 0 .AND. (it==id_rn .OR. it==id_pb)) THEN
524          ! couche limite avec quantite dans le sol calculee
525          CALL cltracrn(it, pdtphys, yu1, yv1,     &
526               cdragh, coefh,t_seri,ftsol,pctsrf,  &
527               tr_seri(:,:,it),trs(:,it),          &
528               paprs, pplay, zmasse * rg, &
529               masktr(:,it),fshtr(:,it),hsoltr(it),&
530               tautr(it),vdeptr(it),               &
531               xlat,d_tr_cl(:,:,it),d_trs)
532         
533          DO k = 1, klev
534             DO i = 1, klon
535                tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cl(i,k,it)
536             END DO
537          END DO
538       
539          ! Traceur dans le reservoir sol
540          DO i = 1, klon
541             trs(i,it) = trs(i,it) + d_trs(i)
542          END DO
543       END IF
544    END DO
545         
546
547!======================================================================
548!   Calcul de l'effet du puits radioactif
549!======================================================================
550    CALL radio_decay (radio,rnpb,pdtphys,tautr,tr_seri,d_tr_dec)
551
552    DO it=1,nbtr
553       WRITE(solsym(it),'(i2)') it
554    END DO
555
556    DO it=1,nbtr
557       IF(radio(it)) then     
558          DO k = 1, klev
559             DO i = 1, klon
560                tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_dec(i,k,it)
561             END DO
562          END DO
563          CALL minmaxqfi(tr_seri(:,:,it),0.,1.e33,'puits rn it='//solsym(it))
564       END IF
565    END DO
566
567!======================================================================
568!   Parameterization of ozone chemistry
569!======================================================================
570
571    IF (id_o3 /= 0) then
572       lmt_pas = NINT(86400./pdtphys)
573       IF (MOD(nstep - 1, lmt_pas) == 0) THEN
574          ! Once per day, update the coefficients for ozone chemistry:
575          CALL regr_pr_comb_coefoz(julien, xlat, paprs, pplay)
576       END IF
577       CALL o3_chem(julien, gmtime, t_seri, zmasse, pdtphys, xlat, &
578            xlon, tr_seri(:, :, id_o3))
579    END IF
580
581  END SUBROUTINE traclmdz
582
583
584  SUBROUTINE traclmdz_to_restart(trs_out)
585    ! This subroutine is called from phyredem.F where the module
586    ! variable trs is written to restart file (restartphy.nc)
587    USE dimphy
588    USE infotrac_phy
589   
590    REAL,DIMENSION(klon,nbtr), INTENT(OUT) :: trs_out
591    INTEGER :: ierr
592
593    IF ( ALLOCATED(trs) ) THEN
594       trs_out(:,:) = trs(:,:)
595    ELSE
596       ! No previous allocate of trs. This is the case for create_etat0_limit.
597       trs_out(:,:) = 0.0
598    END IF
599   
600  END SUBROUTINE traclmdz_to_restart
601 
602
603END MODULE traclmdz_mod
Note: See TracBrowser for help on using the repository browser.