source: LMDZ6/branches/Ocean_skin/libf/phylmd/traclmdz_mod.F90 @ 5434

Last change on this file since 5434 was 4368, checked in by lguez, 2 years ago

Sync latest trunk changes to Ocean_skin

  • 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.1 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, ONLY: nbtr
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, ONLY: nbtr, nqtot, tracers, pbl_flg, conv_flg
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, iq, 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    it = 0
176    DO iq = 1, nqtot
177       IF(.NOT.(tracers(iq)%isInPhysics)) CYCLE
178       it = it+1
179       SELECT CASE(strLower(tracers(iq)%name))
180         CASE("rn");      id_rn     = it ! radon
181         CASE("pb");      id_pb     = it ! plomb
182         CASE("aga");     id_aga    = it ! Age of stratospheric air
183         CASE("be","be7");id_be     = it ! Recherche du Beryllium 7
184         CASE("o3");      id_o3     = it ! Recherche de l'ozone
185         CASE("pcsat");   id_pcsat  = it
186         CASE("pcocsat"); id_pcocsat= it
187         CASE("pcq");     id_pcq    = it
188         CASE("pcs0");    id_pcs0   = it
189         CASE("pcos0");   id_pcos0  = it
190         CASE("pcq0");    id_pcq0   = it
191         CASE DEFAULT
192           WRITE(lunout,*) 'This is an unknown tracer in LMDZ : ', trim(tracers(iq)%name)
193       END SELECT
194
195       SELECT CASE(strLower(tracers(iq)%name))
196         CASE("pb")                        !--- RomP >>> profil initial de PB210
197           OPEN(ilesfil2,file='prof.pb210',status='old',iostat=irr2)
198           IF(irr2 == 0) THEN
199             READ(ilesfil2,*) kradio2
200             WRITE(lunout,*)'number of levels for pb210 profile ',kradio2
201             DO k=kradio2,1,-1; READ (ilesfil2,*) plomb(:,k); END DO
202             CLOSE(ilesfil2)
203             tr_seri(:,:,id_pb) = plomb(:,:)
204           ELSE
205             WRITE(lunout,*)'Prof. Pb210 does not exist: use restart values'
206           END IF
207         CASE("aga")
208           radio(id_aga) = .FALSE.
209           aerosol(id_aga) = .FALSE.
210           pbl_flg(id_aga) = 0
211           ! Find the first model layer above 1.5km from the surface
212           IF (klev>=30) THEN
213              lev_1p5km=6                  !--- NB: This value is for klev=39
214           ELSE IF (klev>=10) THEN
215              lev_1p5km=5                  !--- NB: This value is for klev=19
216           ELSE
217              lev_1p5km=klev/2
218           END IF
219         CASE("be","be7")
220           ALLOCATE( srcbe(klon,klev) )
221           radio(id_be) = .TRUE.
222           aerosol(id_be) = .TRUE.         !--- Le Be est un aerosol
223           CALL init_be(pctsrf,pplay,masktr(:,id_be),tautr(id_be),vdeptr(id_be),scavtr(id_be),srcbe)
224           WRITE(lunout,*) 'Initialisation srcBe: OK'
225                                           !--- RomP >>> profil initial de Be7
226           OPEN(ilesfil,file='prof.be7',status='old',iostat=irr)
227           IF(irr == 0) THEN
228             READ(ilesfil,*) kradio
229             WRITE(lunout,*)'number of levels for Be7 profile ',kradio
230             DO k=kradio,1,-1; READ(ilesfil,*) beryllium(:,k); END DO
231             CLOSE(ilesfil)
232             tr_seri(:,:,id_be)=beryllium(:,:)
233           ELSE
234             WRITE(lunout,*)'Prof. Be7 does not exist: use restart values'
235           END IF
236         CASE("o3")                         !--- Parametrisation par la chimie de Cariolle
237           CALL alloc_coefoz                !--- Allocate ozone coefficients
238           CALL press_coefoz                !--- Read input pressure levels
239         CASE("pcs0","pcos0","pcq0")
240           conv_flg(it)=0                   !--- No transport by convection for this tracer
241       END SELECT
242    END DO
243
244!
245! Valeurs specifiques pour les traceurs Rn222 et Pb210
246! ----------------------------------------------
247    IF ( id_rn/=0 .AND. id_pb/=0 ) THEN
248       rnpb = .TRUE.
249       radio(id_rn)= .TRUE.
250       radio(id_pb)= .TRUE.
251       pbl_flg(id_rn) = 0 ! au lieu de clsol=true ! CL au sol calcule
252       pbl_flg(id_pb) = 0 ! au lieu de clsol=true
253       aerosol(id_rn) = .FALSE.
254       aerosol(id_pb) = .TRUE. ! le Pb est un aerosol
255       
256       CALL initrrnpb (ftsol,pctsrf,masktr,fshtr,hsoltr,tautr,vdeptr,scavtr)
257    END IF
258
259!
260! Check if all tracers have restart values
261! ----------------------------------------------
262    it = 0
263    DO iq = 1, nqtot
264       IF(.NOT.(tracers(iq)%isAdvected .AND. tracers(iq)%isInPhysics)) CYCLE
265       it = it+1
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(iq)%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, ONLY: nbtr, pbl_flg
314    USE strings_mod,  ONLY: int2str
315    USE regr_pr_comb_coefoz_m, ONLY: regr_pr_comb_coefoz
316    USE o3_chem_m, ONLY: o3_chem
317    USE indice_sol_mod
318
319    INCLUDE "YOMCST.h"
320
321!==========================================================================
322!                   -- DESCRIPTION DES ARGUMENTS --
323!==========================================================================
324
325! Input arguments
326!
327!Configuration grille,temps:
328    INTEGER,INTENT(IN) :: nstep      ! nombre d'appels de la physiq
329    INTEGER,INTENT(IN) :: julien     ! Jour julien
330    REAL,INTENT(IN)    :: gmtime
331    REAL,INTENT(IN)    :: pdtphys    ! Pas d'integration pour la physique (seconde) 
332    REAL,DIMENSION(klon),INTENT(IN) :: xlat    ! latitudes pour chaque point
333    REAL, INTENT(IN):: xlon(:) ! dim(klon) longitude
334
335!
336!Physique:
337!--------
338    REAL,DIMENSION(klon,klev),INTENT(IN)   :: t_seri  ! Temperature
339    REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs   ! pression pour chaque inter-couche (en Pa)
340    REAL,DIMENSION(klon,klev),INTENT(IN)   :: pplay   ! pression pour le mileu de chaque couche (en Pa)
341    REAL,intent(in):: zmasse (:, :)   ! dim(klon,klev) density of air, in kg/m2
342
343
344!Couche limite:
345!--------------
346!
347    REAL,DIMENSION(klon),INTENT(IN)      :: cdragh     ! coeff drag pour T et Q
348    REAL,DIMENSION(klon,klev),INTENT(IN) :: coefh      ! coeff melange CL (m**2/s)
349    REAL,DIMENSION(klon),INTENT(IN)      :: yu1        ! vents au premier niveau
350    REAL,DIMENSION(klon),INTENT(IN)      :: yv1        ! vents au premier niveau
351    LOGICAL,INTENT(IN)                   :: couchelimite
352    REAL,DIMENSION(klon,klev),INTENT(IN) :: sh         ! humidite specifique
353    REAL,DIMENSION(klon,klev),INTENT(IN) :: rh      ! Humidite relative
354    REAL,DIMENSION(klon,klev),INTENT(IN) :: pphi    ! geopotentie
355    REAL,DIMENSION(klon),INTENT(IN)      :: ustar   ! ustar (m/s)
356    REAL,DIMENSION(klon),INTENT(IN)      :: wstar,ale_bl,ale_wake   ! wstar (m/s) and Avail. Lifti. Energ.
357    REAL,DIMENSION(klon),INTENT(IN)      :: zu10m   ! vent zonal 10m (m/s)
358    REAL,DIMENSION(klon),INTENT(IN)      :: zv10m   ! vent zonal 10m (m/s)
359
360! Arguments necessaires pour les sources et puits de traceur:
361    REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: ftsol  ! Temperature du sol (surf)(Kelvin)
362    REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: pctsrf ! Pourcentage de sol f(nature du sol)
363
364! InOutput argument
365    REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT)  :: tr_seri ! Concentration Traceur [U/KgA] 
366
367! Output argument
368    REAL,DIMENSION(klon,nbtr), INTENT(OUT)        :: source  ! a voir lorsque le flux de surface est prescrit
369    REAL,DIMENSION(klon,klev,nbtr), INTENT(OUT)   :: d_tr_cl ! Td couche limite/traceur
370
371!=======================================================================================
372!                        -- VARIABLES LOCALES TRACEURS --
373!=======================================================================================
374
375    INTEGER :: i, k, it
376    INTEGER :: lmt_pas ! number of time steps of "physics" per day
377
378    REAL,DIMENSION(klon)           :: d_trs    ! Td dans le reservoir
379    REAL,DIMENSION(klon,klev)      :: qsat     ! pression de la vapeur a saturation
380    REAL,DIMENSION(klon,klev,nbtr) :: d_tr_dec ! Td radioactive
381    REAL                           :: zrho     ! Masse Volumique de l'air KgA/m3
382    REAL                           :: amn, amx
383!
384!=================================================================
385!  Ajout de la production en  Be7 (Beryllium) srcbe U/s/kgA
386!=================================================================
387!
388    IF ( id_be /= 0 ) THEN
389       DO k = 1, klev
390          DO i = 1, klon
391             tr_seri(i,k,id_be) = tr_seri(i,k,id_be)+srcbe(i,k)*pdtphys
392          END DO
393       END DO
394       WRITE(*,*) 'Ajout srcBe dans tr_seri: OK'
395    END IF
396 
397
398!=================================================================
399! Update pseudo-vapor tracers
400!=================================================================
401
402    CALL q_sat(klon*klev,t_seri,pplay,qsat)
403
404    IF ( id_pcsat /= 0 ) THEN
405     DO k = 1, klev
406      DO i = 1, klon
407         IF ( pplay(i,k).GE.85000.) THEN
408            tr_seri(i,k,id_pcsat) = qsat(i,k)
409         ELSE
410            tr_seri(i,k,id_pcsat) = MIN (qsat(i,k), tr_seri(i,k,id_pcsat))           
411         END IF
412      END DO
413     END DO
414    END IF
415
416    IF ( id_pcocsat /= 0 ) THEN
417     DO k = 1, klev
418      DO i = 1, klon
419         IF ( pplay(i,k).GE.85000.) THEN
420            IF ( pctsrf (i, is_oce) > 0. ) THEN
421               tr_seri(i,k,id_pcocsat) = qsat(i,k)
422            ELSE
423               tr_seri(i,k,id_pcocsat) = 0.
424          END IF
425       ELSE
426          tr_seri(i,k,id_pcocsat) = MIN (qsat(i,k), tr_seri(i,k,id_pcocsat))
427       END IF
428      END DO
429     END DO
430    END IF
431
432    IF ( id_pcq /= 0 ) THEN
433     DO k = 1, klev
434      DO i = 1, klon
435         IF ( pplay(i,k).GE.85000.) THEN
436            tr_seri(i,k,id_pcq) = sh(i,k)
437         ELSE
438            tr_seri(i,k,id_pcq) = MIN (qsat(i,k), tr_seri(i,k,id_pcq))   
439         END IF
440      END DO
441     END DO
442    END IF
443
444
445    IF ( id_pcs0 /= 0 ) THEN
446     DO k = 1, klev
447      DO i = 1, klon
448         IF ( pplay(i,k).GE.85000.) THEN
449            tr_seri(i,k,id_pcs0) = qsat(i,k)
450         ELSE
451            tr_seri(i,k,id_pcs0) = MIN (qsat(i,k), tr_seri(i,k,id_pcs0))   
452         END IF
453      END DO
454     END DO
455    END IF
456
457
458    IF ( id_pcos0 /= 0 ) THEN
459     DO k = 1, klev
460      DO i = 1, klon
461         IF ( pplay(i,k).GE.85000.) THEN
462            IF ( pctsrf (i, is_oce) > 0. ) THEN
463               tr_seri(i,k,id_pcos0) = qsat(i,k)
464            ELSE
465               tr_seri(i,k,id_pcos0) = 0.
466            END IF
467         ELSE
468            tr_seri(i,k,id_pcos0) = MIN (qsat(i,k), tr_seri(i,k,id_pcos0))
469         END IF
470      END DO
471     END DO
472    END IF
473
474
475    IF ( id_pcq0 /= 0 ) THEN
476     DO k = 1, klev
477      DO i = 1, klon
478         IF ( pplay(i,k).GE.85000.) THEN
479            tr_seri(i,k,id_pcq0) = sh(i,k)
480         ELSE
481            tr_seri(i,k,id_pcq0) = MIN (qsat(i,k), tr_seri(i,k,id_pcq0))
482         END IF
483      END DO
484     END DO
485    END IF
486
487!=================================================================
488! Update tracer : Age of stratospheric air
489!=================================================================
490    IF (id_aga/=0) THEN
491       
492       ! Bottom layers
493       DO k = 1, lev_1p5km
494          tr_seri(:,k,id_aga) = 0.0
495       END DO
496       
497       ! Layers above 1.5km
498       DO k = lev_1p5km+1,klev-1
499          tr_seri(:,k,id_aga) = tr_seri(:,k,id_aga) + pdtphys
500       END DO
501       
502       ! Top layer
503       tr_seri(:,klev,id_aga) = tr_seri(:,klev-1,id_aga)
504       
505    END IF
506
507!======================================================================
508!     -- Calcul de l'effet de la couche limite --
509!======================================================================
510
511    IF (couchelimite) THEN             
512       source(:,:) = 0.0
513
514       IF (id_be /=0) THEN
515          DO i=1, klon
516             zrho = pplay(i,1)/t_seri(i,1)/RD
517             source(i,id_be) = - vdeptr(id_be)*tr_seri(i,1,id_be)*zrho
518          END DO
519       END IF
520
521    END IF
522   
523    DO it=1, nbtr
524       IF (couchelimite .AND. pbl_flg(it) == 0 .AND. (it==id_rn .OR. it==id_pb)) THEN
525          ! couche limite avec quantite dans le sol calculee
526          CALL cltracrn(it, pdtphys, yu1, yv1,     &
527               cdragh, coefh,t_seri,ftsol,pctsrf,  &
528               tr_seri(:,:,it),trs(:,it),          &
529               paprs, pplay, zmasse * rg, &
530               masktr(:,it),fshtr(:,it),hsoltr(it),&
531               tautr(it),vdeptr(it),               &
532               xlat,d_tr_cl(:,:,it),d_trs)
533         
534          DO k = 1, klev
535             DO i = 1, klon
536                tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cl(i,k,it)
537             END DO
538          END DO
539       
540          ! Traceur dans le reservoir sol
541          DO i = 1, klon
542             trs(i,it) = trs(i,it) + d_trs(i)
543          END DO
544       END IF
545    END DO
546         
547
548!======================================================================
549!   Calcul de l'effet du puits radioactif
550!======================================================================
551    CALL radio_decay (radio,rnpb,pdtphys,tautr,tr_seri,d_tr_dec)
552
553    DO it=1,nbtr
554       IF(radio(it)) then     
555          DO k = 1, klev
556             DO i = 1, klon
557                tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_dec(i,k,it)
558             END DO
559          END DO
560          CALL minmaxqfi(tr_seri(:,:,it),0.,1.e33,'puits rn it='//TRIM(int2str(it)))
561       END IF
562    END DO
563
564!======================================================================
565!   Parameterization of ozone chemistry
566!======================================================================
567
568    IF (id_o3 /= 0) then
569       lmt_pas = NINT(86400./pdtphys)
570       IF (MOD(nstep - 1, lmt_pas) == 0) THEN
571          ! Once per day, update the coefficients for ozone chemistry:
572          CALL regr_pr_comb_coefoz(julien, xlat, paprs, pplay)
573       END IF
574       CALL o3_chem(julien, gmtime, t_seri, zmasse, pdtphys, xlat, &
575            xlon, tr_seri(:, :, id_o3))
576    END IF
577
578  END SUBROUTINE traclmdz
579
580
581  SUBROUTINE traclmdz_to_restart(trs_out)
582    ! This subroutine is called from phyredem.F where the module
583    ! variable trs is written to restart file (restartphy.nc)
584    USE dimphy
585    USE infotrac_phy, ONLY: nbtr
586   
587    REAL,DIMENSION(klon,nbtr), INTENT(OUT) :: trs_out
588    INTEGER :: ierr
589
590    IF ( ALLOCATED(trs) ) THEN
591       trs_out(:,:) = trs(:,:)
592    ELSE
593       ! No previous allocate of trs. This is the case for create_etat0_limit.
594       trs_out(:,:) = 0.0
595    END IF
596   
597  END SUBROUTINE traclmdz_to_restart
598 
599
600END MODULE traclmdz_mod
Note: See TracBrowser for help on using the repository browser.