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

Last change on this file since 3721 was 3581, checked in by oboucher, 5 years ago

Big update to the interactive carbon cycle
from Patricia's code

  • 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.7 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
99! Input variables
100    REAL,DIMENSION(klon,nbsrf),INTENT(IN)     :: pctsrf ! Pourcentage de sol f(nature du sol)
101    REAL,DIMENSION(klon),INTENT(IN)           :: xlat   ! latitudes en degres pour chaque point
102    REAL,DIMENSION(klon),INTENT(IN)           :: xlon   ! longitudes en degres pour chaque point
103    REAL,DIMENSION(klon,nbsrf),INTENT(IN)     :: ftsol  ! Temperature du sol (surf)(Kelvin)
104    REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT) :: tr_seri! Concentration Traceur [U/KgA] 
105    REAL,DIMENSION(klon,klev),INTENT(IN)   :: t_seri  ! Temperature
106    REAL,DIMENSION(klon,klev),INTENT(IN)   :: pplay   ! pression pour le mileu de chaque couche (en Pa)
107    REAL,DIMENSION(klon,klev),INTENT(IN)   :: sh      ! humidite specifique
108    REAL,INTENT(IN)                        :: pdtphys ! Pas d'integration pour la physique (seconde) 
109
110! Output variables
111    LOGICAL,DIMENSION(nbtr), INTENT(OUT) :: aerosol
112    LOGICAL,INTENT(OUT)                  :: lessivage
113       
114! Local variables   
115    INTEGER :: ierr, it, iiq, i, k
116    REAL, DIMENSION(klon_glo,klev) :: varglo ! variable temporaire sur la grille global   
117    REAL, DIMENSION(klev)          :: mintmp, maxtmp
118    LOGICAL                        :: zero
119! RomP >>> profil initial Be7
120      integer ilesfil
121      parameter (ilesfil=1)
122      integer  irr,kradio
123      real     beryllium(klon,klev)
124! profil initial Pb210
125      integer ilesfil2
126      parameter (ilesfil2=1)
127      integer  irr2,kradio2
128      real     plomb(klon,klev)
129!! RomP <<<
130! --------------------------------------------
131! Allocation
132! --------------------------------------------
133    ALLOCATE( scavtr(nbtr), stat=ierr)
134    IF (ierr /= 0) CALL abort_physic('traclmdz_init', 'pb in allocation 9',1)
135    scavtr(:)=1.
136   
137    ALLOCATE( radio(nbtr), stat=ierr)
138    IF (ierr /= 0) CALL abort_physic('traclmdz_init', 'pb in allocation 11',1)
139    radio(:) = .false.    ! Par defaut pas decroissance radioactive
140   
141    ALLOCATE( masktr(klon,nbtr), stat=ierr)
142    IF (ierr /= 0) CALL abort_physic('traclmdz_init', 'pb in allocation 2',1)
143   
144    ALLOCATE( fshtr(klon,nbtr), stat=ierr)
145    IF (ierr /= 0) CALL abort_physic('traclmdz_init', 'pb in allocation 3',1)
146   
147    ALLOCATE( hsoltr(nbtr), stat=ierr)
148    IF (ierr /= 0) CALL abort_physic('traclmdz_init', 'pb in allocation 4',1)
149   
150    ALLOCATE( tautr(nbtr), stat=ierr)
151    IF (ierr /= 0) CALL abort_physic('traclmdz_init', 'pb in allocation 5',1)
152    tautr(:)  = 0.
153   
154    ALLOCATE( vdeptr(nbtr), stat=ierr)
155    IF (ierr /= 0) CALL abort_physic('traclmdz_init', 'pb in allocation 6',1)
156    vdeptr(:) = 0.
157
158
159    lessivage  = .TRUE.
160!!jyg(20130206) : le choix d activation du lessivage est fait dans phytrac avec iflag_lscav
161!!    call getin('lessivage',lessivage)
162!!    if(lessivage) then
163!!     print*,'lessivage lsc ON'
164!!    else
165!!     print*,'lessivage lsc OFF'
166!!    endif
167    aerosol(:) = .FALSE.  ! Tous les traceurs sont des gaz par defaut
168   
169!
170! Recherche des traceurs connus : Be7, O3, CO2,...
171! --------------------------------------------
172    id_rn=0; id_pb=0; id_aga=0; id_be=0; id_o3=0
173    id_pcsat=0; id_pcocsat=0; id_pcq=0; id_pcs0=0; id_pcos0=0; id_pcq0=0
174    DO it=1,nbtr
175!!       iiq=niadv(it+2)                                                            ! jyg
176       iiq=niadv(it+nqo)                                                            ! jyg
177       IF ( tname(iiq) == "RN" ) THEN                                               
178          id_rn=it ! radon
179       ELSE IF ( tname(iiq) == "PB") THEN
180          id_pb=it ! plomb
181! RomP >>> profil initial de PB210
182     open (ilesfil2,file='prof.pb210',status='old',iostat=irr2)
183     IF (irr2 == 0) THEN
184      read(ilesfil2,*) kradio2
185      print*,'number of levels for pb210 profile ',kradio2
186      do k=kradio2,1,-1
187       read (ilesfil2,*) plomb(:,k)
188      enddo
189      close(ilesfil2)
190      do k=1,klev
191       do i=1,klon
192         tr_seri(i,k,id_pb)=plomb(i,k)
193!!        print*, 'tr_seri',i,k,tr_seri(i,k,id_pb)
194        enddo
195      enddo
196     ELSE
197       print *, 'Prof.pb210 does not exist: use restart values'
198     ENDIF
199! RomP <<<
200       ELSE IF ( tname(iiq) == "Aga" .OR. tname(iiq) == "AGA" ) THEN
201          ! Age of stratospheric air
202          id_aga=it
203          radio(id_aga) = .FALSE.
204          aerosol(id_aga) = .FALSE.
205          pbl_flg(id_aga) = 0
206         
207          ! Find the first model layer above 1.5km from the surface
208          IF (klev>=30) THEN
209             lev_1p5km=6   ! NB! This value is for klev=39
210          ELSE IF (klev>=10) THEN
211             lev_1p5km=5   ! NB! This value is for klev=19
212          ELSE
213             lev_1p5km=klev/2
214          END IF
215       ELSE IF ( tname(iiq) == "BE" .OR. tname(iiq) == "Be" .OR.  &
216            tname(iiq) == "BE7" .OR. tname(iiq) == "Be7" ) THEN 
217          ! Recherche du Beryllium 7
218          id_be=it
219          ALLOCATE( srcbe(klon,klev) )
220          radio(id_be) = .TRUE.
221          aerosol(id_be) = .TRUE. ! le Be est un aerosol
222!jyg le 13/03/2013 ; ajout de pplay en argument de init_be
223!!!          CALL init_be(pctsrf,masktr(:,id_be),tautr(id_be),vdeptr(id_be),scavtr(id_be),srcbe)
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       print*,'number of levels for Be7 profile ',kradio
231       do k=kradio,1,-1
232        read (ilesfil,*) beryllium(:,k)
233       enddo
234       close(ilesfil)
235       do k=1,klev
236         do i=1,klon
237          tr_seri(i,k,id_be)=beryllium(i,k)
238!!        print*, 'tr_seri',i,k,tr_seri(i,k,id_be)
239         enddo
240       enddo
241     ELSE
242       print *, 'Prof.Be7 does not exist: use restart values'
243     ENDIF
244! RomP <<<
245       ELSE IF (tname(iiq)=="O3" .OR. tname(iiq)=="o3") THEN
246          ! Recherche de l'ozone : parametrization de la chimie par Cariolle
247          id_o3=it
248          CALL alloc_coefoz   ! allocate ozone coefficients
249          CALL press_coefoz   ! read input pressure levels
250       ELSE IF ( tname(iiq) == "pcsat" .OR. tname(iiq) == "Pcsat" ) THEN
251          id_pcsat=it
252       ELSE IF ( tname(iiq) == "pcocsat" .OR. tname(iiq) == "Pcocsat" ) THEN
253          id_pcocsat=it
254       ELSE IF ( tname(iiq) == "pcq" .OR. tname(iiq) == "Pcq" ) THEN
255          id_pcq=it
256       ELSE IF ( tname(iiq) == "pcs0" .OR. tname(iiq) == "Pcs0" ) THEN
257          id_pcs0=it
258          conv_flg(it)=0 ! No transport by convection for this tracer
259       ELSE IF ( tname(iiq) == "pcos0" .OR. tname(iiq) == "Pcos0" ) THEN
260          id_pcos0=it
261          conv_flg(it)=0 ! No transport by convection for this tracer
262       ELSE IF ( tname(iiq) == "pcq0" .OR. tname(iiq) == "Pcq0" ) THEN
263          id_pcq0=it
264          conv_flg(it)=0 ! No transport by convection for this tracer
265       ELSE
266          WRITE(lunout,*) 'This is an unknown tracer in LMDZ : ', trim(tname(iiq))
267       END IF
268    END DO
269
270!
271! Valeurs specifiques pour les traceurs Rn222 et Pb210
272! ----------------------------------------------
273    IF ( id_rn/=0 .AND. id_pb/=0 ) THEN
274       rnpb = .TRUE.
275       radio(id_rn)= .TRUE.
276       radio(id_pb)= .TRUE.
277       pbl_flg(id_rn) = 0 ! au lieu de clsol=true ! CL au sol calcule
278       pbl_flg(id_pb) = 0 ! au lieu de clsol=true
279       aerosol(id_rn) = .FALSE.
280       aerosol(id_pb) = .TRUE. ! le Pb est un aerosol
281       
282       CALL initrrnpb (ftsol,pctsrf,masktr,fshtr,hsoltr,tautr,vdeptr,scavtr)
283    END IF
284
285!
286! Check if all tracers have restart values
287! ----------------------------------------------
288    DO it=1,nbtr
289!!       iiq=niadv(it+2)                                                            ! jyg
290       iiq=niadv(it+nqo)                                                            ! jyg
291       ! Test if tracer is zero everywhere.
292       ! Done by master process MPI and master thread OpenMP
293       CALL gather(tr_seri(:,:,it),varglo)
294       IF (is_mpi_root .AND. is_omp_root) THEN
295          mintmp=MINVAL(varglo,dim=1)
296          maxtmp=MAXVAL(varglo,dim=1)
297          IF (MINVAL(mintmp,dim=1)==0. .AND. MAXVAL(maxtmp,dim=1)==0.) THEN
298             ! Tracer is zero everywhere
299             zero=.TRUE.
300          ELSE
301             zero=.FALSE.
302          END IF
303       END IF
304
305       ! Distribute variable at all processes
306       CALL bcast(zero)
307
308       ! Initalize tracer that was not found in restart file.
309       IF (zero) THEN
310          ! The tracer was not found in restart file or it was equal zero everywhere.
311          WRITE(lunout,*) "The tracer ",trim(tname(iiq))," will be initialized"
312          IF (it==id_pcsat .OR. it==id_pcq .OR. &
313               it==id_pcs0 .OR. it==id_pcq0) THEN
314             tr_seri(:,:,it) = 100.
315          ELSE IF (it==id_pcocsat .OR. it==id_pcos0) THEN
316             DO i = 1, klon
317                IF ( pctsrf (i, is_oce) == 0. ) THEN
318                   tr_seri(i,:,it) = 0.
319                ELSE
320                   tr_seri(i,:,it) = 100.
321                END IF
322             END DO
323          ELSE
324             ! No specific initialization exist for this tracer
325             tr_seri(:,:,it) = 0.
326          END IF
327       END IF
328    END DO
329
330  END SUBROUTINE traclmdz_init
331
332  SUBROUTINE traclmdz(nstep, julien, gmtime, pdtphys, t_seri, paprs, pplay, &
333       cdragh,  coefh, yu1, yv1, ftsol, pctsrf, xlat, xlon, couchelimite, sh, &
334       rh, pphi, ustar, wstar, ale_bl, ale_wake,  zu10m, zv10m, &
335       tr_seri, source, d_tr_cl,d_tr_dec, zmasse)               !RomP
336   
337    USE dimphy
338    USE infotrac_phy
339    USE regr_pr_comb_coefoz_m, ONLY: regr_pr_comb_coefoz
340    USE o3_chem_m, ONLY: o3_chem
341    USE indice_sol_mod
342
343    INCLUDE "YOMCST.h"
344
345!==========================================================================
346!                   -- DESCRIPTION DES ARGUMENTS --
347!==========================================================================
348
349! Input arguments
350!
351!Configuration grille,temps:
352    INTEGER,INTENT(IN) :: nstep      ! nombre d'appels de la physiq
353    INTEGER,INTENT(IN) :: julien     ! Jour julien
354    REAL,INTENT(IN)    :: gmtime
355    REAL,INTENT(IN)    :: pdtphys    ! Pas d'integration pour la physique (seconde) 
356    REAL,DIMENSION(klon),INTENT(IN) :: xlat    ! latitudes pour chaque point
357    REAL, INTENT(IN):: xlon(:) ! dim(klon) longitude
358
359!
360!Physique:
361!--------
362    REAL,DIMENSION(klon,klev),INTENT(IN)   :: t_seri  ! Temperature
363    REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs   ! pression pour chaque inter-couche (en Pa)
364    REAL,DIMENSION(klon,klev),INTENT(IN)   :: pplay   ! pression pour le mileu de chaque couche (en Pa)
365    REAL,intent(in):: zmasse (:, :)   ! dim(klon,klev) density of air, in kg/m2
366
367
368!Couche limite:
369!--------------
370!
371    REAL,DIMENSION(klon),INTENT(IN)      :: cdragh     ! coeff drag pour T et Q
372    REAL,DIMENSION(klon,klev),INTENT(IN) :: coefh      ! coeff melange CL (m**2/s)
373    REAL,DIMENSION(klon),INTENT(IN)      :: yu1        ! vents au premier niveau
374    REAL,DIMENSION(klon),INTENT(IN)      :: yv1        ! vents au premier niveau
375    LOGICAL,INTENT(IN)                   :: couchelimite
376    REAL,DIMENSION(klon,klev),INTENT(IN) :: sh         ! humidite specifique
377    REAL,DIMENSION(klon,klev),INTENT(IN) :: rh      ! Humidite relative
378    REAL,DIMENSION(klon,klev),INTENT(IN) :: pphi    ! geopotentie
379    REAL,DIMENSION(klon),INTENT(IN)      :: ustar   ! ustar (m/s)
380    REAL,DIMENSION(klon),INTENT(IN)      :: wstar,ale_bl,ale_wake   ! wstar (m/s) and Avail. Lifti. Energ.
381    REAL,DIMENSION(klon),INTENT(IN)      :: zu10m   ! vent zonal 10m (m/s)
382    REAL,DIMENSION(klon),INTENT(IN)      :: zv10m   ! vent zonal 10m (m/s)
383
384! Arguments necessaires pour les sources et puits de traceur:
385    REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: ftsol  ! Temperature du sol (surf)(Kelvin)
386    REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: pctsrf ! Pourcentage de sol f(nature du sol)
387
388! InOutput argument
389    REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT)  :: tr_seri ! Concentration Traceur [U/KgA] 
390
391! Output argument
392    REAL,DIMENSION(klon,nbtr), INTENT(OUT)        :: source  ! a voir lorsque le flux de surface est prescrit
393    REAL,DIMENSION(klon,klev,nbtr), INTENT(OUT)   :: d_tr_cl ! Td couche limite/traceur
394
395!=======================================================================================
396!                        -- VARIABLES LOCALES TRACEURS --
397!=======================================================================================
398
399    INTEGER :: i, k, it
400    INTEGER :: lmt_pas ! number of time steps of "physics" per day
401
402    REAL,DIMENSION(klon)           :: d_trs    ! Td dans le reservoir
403    REAL,DIMENSION(klon,klev)      :: qsat     ! pression de la vapeur a saturation
404    REAL,DIMENSION(klon,klev,nbtr) :: d_tr_dec ! Td radioactive
405    REAL                           :: zrho     ! Masse Volumique de l'air KgA/m3
406    REAL                           :: amn, amx
407!
408!=================================================================
409!  Ajout de la production en  Be7 (Beryllium) srcbe U/s/kgA
410!=================================================================
411!
412    IF ( id_be /= 0 ) THEN
413       DO k = 1, klev
414          DO i = 1, klon
415             tr_seri(i,k,id_be) = tr_seri(i,k,id_be)+srcbe(i,k)*pdtphys
416          END DO
417       END DO
418       WRITE(*,*) 'Ajout srcBe dans tr_seri: OK'
419    END IF
420 
421
422!=================================================================
423! Update pseudo-vapor tracers
424!=================================================================
425
426    CALL q_sat(klon*klev,t_seri,pplay,qsat)
427
428    IF ( id_pcsat /= 0 ) THEN
429     DO k = 1, klev
430      DO i = 1, klon
431         IF ( pplay(i,k).GE.85000.) THEN
432            tr_seri(i,k,id_pcsat) = qsat(i,k)
433         ELSE
434            tr_seri(i,k,id_pcsat) = MIN (qsat(i,k), tr_seri(i,k,id_pcsat))           
435         END IF
436      END DO
437     END DO
438    END IF
439
440    IF ( id_pcocsat /= 0 ) THEN
441     DO k = 1, klev
442      DO i = 1, klon
443         IF ( pplay(i,k).GE.85000.) THEN
444            IF ( pctsrf (i, is_oce) > 0. ) THEN
445               tr_seri(i,k,id_pcocsat) = qsat(i,k)
446            ELSE
447               tr_seri(i,k,id_pcocsat) = 0.
448          END IF
449       ELSE
450          tr_seri(i,k,id_pcocsat) = MIN (qsat(i,k), tr_seri(i,k,id_pcocsat))
451       END IF
452      END DO
453     END DO
454    END IF
455
456    IF ( id_pcq /= 0 ) THEN
457     DO k = 1, klev
458      DO i = 1, klon
459         IF ( pplay(i,k).GE.85000.) THEN
460            tr_seri(i,k,id_pcq) = sh(i,k)
461         ELSE
462            tr_seri(i,k,id_pcq) = MIN (qsat(i,k), tr_seri(i,k,id_pcq))   
463         END IF
464      END DO
465     END DO
466    END IF
467
468
469    IF ( id_pcs0 /= 0 ) THEN
470     DO k = 1, klev
471      DO i = 1, klon
472         IF ( pplay(i,k).GE.85000.) THEN
473            tr_seri(i,k,id_pcs0) = qsat(i,k)
474         ELSE
475            tr_seri(i,k,id_pcs0) = MIN (qsat(i,k), tr_seri(i,k,id_pcs0))   
476         END IF
477      END DO
478     END DO
479    END IF
480
481
482    IF ( id_pcos0 /= 0 ) THEN
483     DO k = 1, klev
484      DO i = 1, klon
485         IF ( pplay(i,k).GE.85000.) THEN
486            IF ( pctsrf (i, is_oce) > 0. ) THEN
487               tr_seri(i,k,id_pcos0) = qsat(i,k)
488            ELSE
489               tr_seri(i,k,id_pcos0) = 0.
490            END IF
491         ELSE
492            tr_seri(i,k,id_pcos0) = MIN (qsat(i,k), tr_seri(i,k,id_pcos0))
493         END IF
494      END DO
495     END DO
496    END IF
497
498
499    IF ( id_pcq0 /= 0 ) THEN
500     DO k = 1, klev
501      DO i = 1, klon
502         IF ( pplay(i,k).GE.85000.) THEN
503            tr_seri(i,k,id_pcq0) = sh(i,k)
504         ELSE
505            tr_seri(i,k,id_pcq0) = MIN (qsat(i,k), tr_seri(i,k,id_pcq0))
506         END IF
507      END DO
508     END DO
509    END IF
510
511!=================================================================
512! Update tracer : Age of stratospheric air
513!=================================================================
514    IF (id_aga/=0) THEN
515       
516       ! Bottom layers
517       DO k = 1, lev_1p5km
518          tr_seri(:,k,id_aga) = 0.0
519       END DO
520       
521       ! Layers above 1.5km
522       DO k = lev_1p5km+1,klev-1
523          tr_seri(:,k,id_aga) = tr_seri(:,k,id_aga) + pdtphys
524       END DO
525       
526       ! Top layer
527       tr_seri(:,klev,id_aga) = tr_seri(:,klev-1,id_aga)
528       
529    END IF
530
531!======================================================================
532!     -- Calcul de l'effet de la couche limite --
533!======================================================================
534
535    IF (couchelimite) THEN             
536       source(:,:) = 0.0
537
538       IF (id_be /=0) THEN
539          DO i=1, klon
540             zrho = pplay(i,1)/t_seri(i,1)/RD
541             source(i,id_be) = - vdeptr(id_be)*tr_seri(i,1,id_be)*zrho
542          END DO
543       END IF
544
545    END IF
546   
547    DO it=1, nbtr
548       IF (couchelimite .AND. pbl_flg(it) == 0 .AND. (it==id_rn .OR. it==id_pb)) THEN
549          ! couche limite avec quantite dans le sol calculee
550          CALL cltracrn(it, pdtphys, yu1, yv1,     &
551               cdragh, coefh,t_seri,ftsol,pctsrf,  &
552               tr_seri(:,:,it),trs(:,it),          &
553               paprs, pplay, zmasse * rg, &
554               masktr(:,it),fshtr(:,it),hsoltr(it),&
555               tautr(it),vdeptr(it),               &
556               xlat,d_tr_cl(:,:,it),d_trs)
557         
558          DO k = 1, klev
559             DO i = 1, klon
560                tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cl(i,k,it)
561             END DO
562          END DO
563       
564          ! Traceur dans le reservoir sol
565          DO i = 1, klon
566             trs(i,it) = trs(i,it) + d_trs(i)
567          END DO
568       END IF
569    END DO
570         
571
572!======================================================================
573!   Calcul de l'effet du puits radioactif
574!======================================================================
575    CALL radio_decay (radio,rnpb,pdtphys,tautr,tr_seri,d_tr_dec)
576
577    DO it=1,nbtr
578       WRITE(solsym(it),'(i2)') it
579    END DO
580
581    DO it=1,nbtr
582       IF(radio(it)) then     
583          DO k = 1, klev
584             DO i = 1, klon
585                tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_dec(i,k,it)
586             END DO
587          END DO
588          CALL minmaxqfi(tr_seri(:,:,it),0.,1.e33,'puits rn it='//solsym(it))
589       END IF
590    END DO
591
592!======================================================================
593!   Parameterization of ozone chemistry
594!======================================================================
595
596    IF (id_o3 /= 0) then
597       lmt_pas = NINT(86400./pdtphys)
598       IF (MOD(nstep - 1, lmt_pas) == 0) THEN
599          ! Once per day, update the coefficients for ozone chemistry:
600          CALL regr_pr_comb_coefoz(julien, xlat, paprs, pplay)
601       END IF
602       CALL o3_chem(julien, gmtime, t_seri, zmasse, pdtphys, xlat, &
603            xlon, tr_seri(:, :, id_o3))
604    END IF
605
606  END SUBROUTINE traclmdz
607
608
609  SUBROUTINE traclmdz_to_restart(trs_out)
610    ! This subroutine is called from phyredem.F where the module
611    ! variable trs is written to restart file (restartphy.nc)
612    USE dimphy
613    USE infotrac_phy
614   
615    REAL,DIMENSION(klon,nbtr), INTENT(OUT) :: trs_out
616    INTEGER :: ierr
617
618    IF ( ALLOCATED(trs) ) THEN
619       trs_out(:,:) = trs(:,:)
620    ELSE
621       ! No previous allocate of trs. This is the case for create_etat0_limit.
622       trs_out(:,:) = 0.0
623    END IF
624   
625  END SUBROUTINE traclmdz_to_restart
626 
627
628END MODULE traclmdz_mod
Note: See TracBrowser for help on using the repository browser.