source: LMDZ6/branches/LMDZ_ECRad/libf/phylmd/traclmdz_mod.F90 @ 5451

Last change on this file since 5451 was 4482, checked in by lguez, 21 months ago

Sync latest trunk changes to branch LMDZ_ECRad

  • 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.2 KB
RevLine 
[1191]1!$Id $
2!
3MODULE traclmdz_mod
[1742]4
[1191]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!
[1403]9  IMPLICIT NONE
10
[1191]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
[1409]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
[1191]44  INTEGER,SAVE :: id_be       ! Activation et position du traceur Be7 [ id_be=0 -> desactive ]
45!$OMP THREADPRIVATE(id_be)
46
[1403]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)
[1421]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
[1403]51!$OMP THREADPRIVATE(id_pcs0, id_pcos0, id_pcq0)
52
53  INTEGER, SAVE:: id_o3
[1421]54!$OMP THREADPRIVATE(id_o3)
55! index of ozone tracer with Cariolle parameterization
56! 0 means no ozone tracer
[1403]57
[1409]58  LOGICAL,SAVE :: rnpb=.FALSE. ! Presence du couple Rn222, Pb210
[1191]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
[4056]69    USE infotrac_phy, ONLY: nbtr
[1191]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)
[2311]79    IF (ierr /= 0) CALL abort_physic('traclmdz_from_restart', 'pb in allocation 1',1)
[1191]80   
81    ! Initialize trs with values read from restart file
82    trs(:,:) = trs_in(:,:)
83   
84  END SUBROUTINE traclmdz_from_restart
85
86
[1579]87  SUBROUTINE traclmdz_init(pctsrf, xlat, xlon, ftsol, tr_seri, t_seri, pplay, sh, pdtphys, aerosol, lessivage)
[1191]88    ! This subroutine allocates and initialize module variables and control variables.
[1421]89    ! Initialization of the tracers should be done here only for those not found in the restart file.
[1191]90    USE dimphy
[4056]91    USE infotrac_phy, ONLY: nbtr, nqtot, tracers, pbl_flg, conv_flg
[1403]92    USE regr_pr_comb_coefoz_m, ONLY: alloc_coefoz
93    USE press_coefoz_m, ONLY: press_coefoz
[1421]94    USE mod_grid_phy_lmdz
95    USE mod_phys_lmdz_para
[1785]96    USE indice_sol_mod
[2311]97    USE print_control_mod, ONLY: lunout
[1191]98
99! Input variables
[1227]100    REAL,DIMENSION(klon,nbsrf),INTENT(IN)     :: pctsrf ! Pourcentage de sol f(nature du sol)
[1579]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
[1227]103    REAL,DIMENSION(klon,nbsrf),INTENT(IN)     :: ftsol  ! Temperature du sol (surf)(Kelvin)
[1403]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
[1454]108    REAL,INTENT(IN)                        :: pdtphys ! Pas d'integration pour la physique (seconde) 
[1227]109
[1191]110! Output variables
111    LOGICAL,DIMENSION(nbtr), INTENT(OUT) :: aerosol
112    LOGICAL,INTENT(OUT)                  :: lessivage
113       
114! Local variables   
[4056]115    INTEGER :: ierr, it, iq, i, k
[1421]116    REAL, DIMENSION(klon_glo,klev) :: varglo ! variable temporaire sur la grille global   
117    REAL, DIMENSION(klev)          :: mintmp, maxtmp
118    LOGICAL                        :: zero
[1742]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 <<<
[1191]130! --------------------------------------------
131! Allocation
132! --------------------------------------------
133    ALLOCATE( scavtr(nbtr), stat=ierr)
[2311]134    IF (ierr /= 0) CALL abort_physic('traclmdz_init', 'pb in allocation 9',1)
[1191]135    scavtr(:)=1.
136   
137    ALLOCATE( radio(nbtr), stat=ierr)
[2311]138    IF (ierr /= 0) CALL abort_physic('traclmdz_init', 'pb in allocation 11',1)
[1191]139    radio(:) = .false.    ! Par defaut pas decroissance radioactive
140   
141    ALLOCATE( masktr(klon,nbtr), stat=ierr)
[2311]142    IF (ierr /= 0) CALL abort_physic('traclmdz_init', 'pb in allocation 2',1)
[1191]143   
144    ALLOCATE( fshtr(klon,nbtr), stat=ierr)
[2311]145    IF (ierr /= 0) CALL abort_physic('traclmdz_init', 'pb in allocation 3',1)
[1191]146   
147    ALLOCATE( hsoltr(nbtr), stat=ierr)
[2311]148    IF (ierr /= 0) CALL abort_physic('traclmdz_init', 'pb in allocation 4',1)
[1191]149   
150    ALLOCATE( tautr(nbtr), stat=ierr)
[2311]151    IF (ierr /= 0) CALL abort_physic('traclmdz_init', 'pb in allocation 5',1)
[1191]152    tautr(:)  = 0.
153   
154    ALLOCATE( vdeptr(nbtr), stat=ierr)
[2311]155    IF (ierr /= 0) CALL abort_physic('traclmdz_init', 'pb in allocation 6',1)
[1191]156    vdeptr(:) = 0.
157
158
159    lessivage  = .TRUE.
[1742]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
[1191]167    aerosol(:) = .FALSE.  ! Tous les traceurs sont des gaz par defaut
168   
169!
[1403]170! Recherche des traceurs connus : Be7, O3, CO2,...
[1191]171! --------------------------------------------
[1409]172    id_rn=0; id_pb=0; id_aga=0; id_be=0; id_o3=0
[1421]173    id_pcsat=0; id_pcocsat=0; id_pcq=0; id_pcs0=0; id_pcos0=0; id_pcq0=0
[4056]174    it = 0
175    DO iq = 1, nqtot
[4482]176       IF(.NOT.(tracers(iq)%isInPhysics)) CYCLE
[4056]177       it = it+1
[4482]178       SELECT CASE(tracers(iq)%name)
179         CASE("RN");                  id_rn     = it ! radon
180         CASE("PB");                  id_pb     = it ! plomb
181         CASE("Aga","AGA");           id_aga    = it ! Age of stratospheric air
182         CASE("Be","BE","Be7","BE7"); id_be     = it ! Recherche du Beryllium 7
183         CASE("o3","O3");             id_o3     = it ! Recherche de l'ozone
184         CASE("pcsat",  "Pcsat");     id_pcsat  = it
185         CASE("pcocsat","Pcocsat");   id_pcocsat= it
186         CASE("pcq",    "Pcq");       id_pcq    = it
187         CASE("pcs0",   "Pcs0");      id_pcs0   = it
188         CASE("pcos0",  "Pcos0");     id_pcos0  = it
189         CASE("pcq0",   "Pcq0");      id_pcq0   = it
[4046]190         CASE DEFAULT
[4056]191           WRITE(lunout,*) 'This is an unknown tracer in LMDZ : ', trim(tracers(iq)%name)
[4046]192       END SELECT
193
[4482]194       SELECT CASE(tracers(iq)%name)
195         CASE("PB")                        !--- RomP >>> profil initial de PB210
[4046]196           OPEN(ilesfil2,file='prof.pb210',status='old',iostat=irr2)
197           IF(irr2 == 0) THEN
198             READ(ilesfil2,*) kradio2
199             WRITE(lunout,*)'number of levels for pb210 profile ',kradio2
200             DO k=kradio2,1,-1; READ (ilesfil2,*) plomb(:,k); END DO
201             CLOSE(ilesfil2)
202             tr_seri(:,:,id_pb) = plomb(:,:)
203           ELSE
[4482]204             WRITE(lunout,*)'prof.pb210 does not exist: use restart values'
[4046]205           END IF
[4482]206         CASE("Aga","AGA")
[4046]207           radio(id_aga) = .FALSE.
208           aerosol(id_aga) = .FALSE.
209           pbl_flg(id_aga) = 0
210           ! Find the first model layer above 1.5km from the surface
211           IF (klev>=30) THEN
212              lev_1p5km=6                  !--- NB: This value is for klev=39
213           ELSE IF (klev>=10) THEN
214              lev_1p5km=5                  !--- NB: This value is for klev=19
215           ELSE
216              lev_1p5km=klev/2
217           END IF
[4482]218         CASE("Be","BE","Be7","BE7")
[4046]219           ALLOCATE( srcbe(klon,klev) )
220           radio(id_be) = .TRUE.
221           aerosol(id_be) = .TRUE.         !--- Le Be est un aerosol
222           CALL init_be(pctsrf,pplay,masktr(:,id_be),tautr(id_be),vdeptr(id_be),scavtr(id_be),srcbe)
223           WRITE(lunout,*) 'Initialisation srcBe: OK'
224                                           !--- RomP >>> profil initial de Be7
225           OPEN(ilesfil,file='prof.be7',status='old',iostat=irr)
226           IF(irr == 0) THEN
227             READ(ilesfil,*) kradio
228             WRITE(lunout,*)'number of levels for Be7 profile ',kradio
229             DO k=kradio,1,-1; READ(ilesfil,*) beryllium(:,k); END DO
230             CLOSE(ilesfil)
231             tr_seri(:,:,id_be)=beryllium(:,:)
232           ELSE
233             WRITE(lunout,*)'Prof. Be7 does not exist: use restart values'
234           END IF
[4482]235         CASE("o3","O3")                    !--- Parametrisation par la chimie de Cariolle
[4046]236           CALL alloc_coefoz                !--- Allocate ozone coefficients
237           CALL press_coefoz                !--- Read input pressure levels
[4482]238         CASE("pcs0","Pcs0", "pcos0","Pcos0", "pcq0","Pcq0")
[4046]239           conv_flg(it)=0                   !--- No transport by convection for this tracer
240       END SELECT
[1403]241    END DO
242
[1191]243!
244! Valeurs specifiques pour les traceurs Rn222 et Pb210
245! ----------------------------------------------
[1409]246    IF ( id_rn/=0 .AND. id_pb/=0 ) THEN
247       rnpb = .TRUE.
248       radio(id_rn)= .TRUE.
249       radio(id_pb)= .TRUE.
250       pbl_flg(id_rn) = 0 ! au lieu de clsol=true ! CL au sol calcule
251       pbl_flg(id_pb) = 0 ! au lieu de clsol=true
252       aerosol(id_rn) = .FALSE.
253       aerosol(id_pb) = .TRUE. ! le Pb est un aerosol
[1191]254       
255       CALL initrrnpb (ftsol,pctsrf,masktr,fshtr,hsoltr,tautr,vdeptr,scavtr)
256    END IF
257
[1227]258!
[1421]259! Check if all tracers have restart values
260! ----------------------------------------------
[4056]261    it = 0
262    DO iq = 1, nqtot
[4071]263       IF(.NOT.(tracers(iq)%isAdvected .AND. tracers(iq)%isInPhysics)) CYCLE
[4056]264       it = it+1
[1421]265       ! Test if tracer is zero everywhere.
266       ! Done by master process MPI and master thread OpenMP
267       CALL gather(tr_seri(:,:,it),varglo)
268       IF (is_mpi_root .AND. is_omp_root) THEN
269          mintmp=MINVAL(varglo,dim=1)
270          maxtmp=MAXVAL(varglo,dim=1)
271          IF (MINVAL(mintmp,dim=1)==0. .AND. MAXVAL(maxtmp,dim=1)==0.) THEN
272             ! Tracer is zero everywhere
273             zero=.TRUE.
274          ELSE
275             zero=.FALSE.
276          END IF
[1403]277       END IF
278
[1421]279       ! Distribute variable at all processes
280       CALL bcast(zero)
[1403]281
[1421]282       ! Initalize tracer that was not found in restart file.
283       IF (zero) THEN
284          ! The tracer was not found in restart file or it was equal zero everywhere.
[4056]285          WRITE(lunout,*) "The tracer ",trim(tracers(iq)%name)," will be initialized"
[1421]286          IF (it==id_pcsat .OR. it==id_pcq .OR. &
287               it==id_pcs0 .OR. it==id_pcq0) THEN
288             tr_seri(:,:,it) = 100.
289          ELSE IF (it==id_pcocsat .OR. it==id_pcos0) THEN
290             DO i = 1, klon
291                IF ( pctsrf (i, is_oce) == 0. ) THEN
292                   tr_seri(i,:,it) = 0.
293                ELSE
294                   tr_seri(i,:,it) = 100.
295                END IF
296             END DO
297          ELSE
298             ! No specific initialization exist for this tracer
299             tr_seri(:,:,it) = 0.
300          END IF
[1403]301       END IF
[1421]302    END DO
[1403]303
[1191]304  END SUBROUTINE traclmdz_init
305
[1403]306  SUBROUTINE traclmdz(nstep, julien, gmtime, pdtphys, t_seri, paprs, pplay, &
307       cdragh,  coefh, yu1, yv1, ftsol, pctsrf, xlat, xlon, couchelimite, sh, &
[1813]308       rh, pphi, ustar, wstar, ale_bl, ale_wake,  zu10m, zv10m, &
[2180]309       tr_seri, source, d_tr_cl,d_tr_dec, zmasse)               !RomP
[1191]310   
311    USE dimphy
[4124]312    USE infotrac_phy, ONLY: nbtr, pbl_flg
313    USE strings_mod,  ONLY: int2str
[1403]314    USE regr_pr_comb_coefoz_m, ONLY: regr_pr_comb_coefoz
315    USE o3_chem_m, ONLY: o3_chem
[1785]316    USE indice_sol_mod
317
[1191]318    INCLUDE "YOMCST.h"
319
320!==========================================================================
321!                   -- DESCRIPTION DES ARGUMENTS --
322!==========================================================================
323
324! Input arguments
325!
326!Configuration grille,temps:
[1212]327    INTEGER,INTENT(IN) :: nstep      ! nombre d'appels de la physiq
[1403]328    INTEGER,INTENT(IN) :: julien     ! Jour julien
329    REAL,INTENT(IN)    :: gmtime
[1191]330    REAL,INTENT(IN)    :: pdtphys    ! Pas d'integration pour la physique (seconde) 
331    REAL,DIMENSION(klon),INTENT(IN) :: xlat    ! latitudes pour chaque point
[1403]332    REAL, INTENT(IN):: xlon(:) ! dim(klon) longitude
[1191]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)
[1403]340    REAL,intent(in):: zmasse (:, :)   ! dim(klon,klev) density of air, in kg/m2
[1191]341
342
343!Couche limite:
344!--------------
345!
[1742]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
[1191]350    LOGICAL,INTENT(IN)                   :: couchelimite
[1742]351    REAL,DIMENSION(klon,klev),INTENT(IN) :: sh         ! humidite specifique
[1670]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)
[1813]355    REAL,DIMENSION(klon),INTENT(IN)      :: wstar,ale_bl,ale_wake   ! wstar (m/s) and Avail. Lifti. Energ.
[1670]356    REAL,DIMENSION(klon),INTENT(IN)      :: zu10m   ! vent zonal 10m (m/s)
357    REAL,DIMENSION(klon),INTENT(IN)      :: zv10m   ! vent zonal 10m (m/s)
[1191]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
[1421]375    INTEGER :: lmt_pas ! number of time steps of "physics" per day
[1191]376
377    REAL,DIMENSION(klon)           :: d_trs    ! Td dans le reservoir
[1421]378    REAL,DIMENSION(klon,klev)      :: qsat     ! pression de la vapeur a saturation
[1191]379    REAL,DIMENSION(klon,klev,nbtr) :: d_tr_dec ! Td radioactive
[1421]380    REAL                           :: zrho     ! Masse Volumique de l'air KgA/m3
381    REAL                           :: amn, amx
[1191]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 
[1421]396
397!=================================================================
398! Update pseudo-vapor tracers
399!=================================================================
400
401    CALL q_sat(klon*klev,t_seri,pplay,qsat)
402
[1403]403    IF ( id_pcsat /= 0 ) THEN
404     DO k = 1, klev
405      DO i = 1, klon
[1421]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
[1403]411      END DO
412     END DO
413    END IF
[1191]414
[1403]415    IF ( id_pcocsat /= 0 ) THEN
416     DO k = 1, klev
417      DO i = 1, klon
[1421]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))
[1403]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
[1421]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
[1403]439      END DO
440     END DO
441    END IF
442
[1421]443
[1403]444    IF ( id_pcs0 /= 0 ) THEN
445     DO k = 1, klev
446      DO i = 1, klon
[1421]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
[1403]452      END DO
453     END DO
454    END IF
455
[1421]456
[1403]457    IF ( id_pcos0 /= 0 ) THEN
458     DO k = 1, klev
459      DO i = 1, klon
[1421]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
[1403]469      END DO
470     END DO
471    END IF
472
[1421]473
[1403]474    IF ( id_pcq0 /= 0 ) THEN
475     DO k = 1, klev
476      DO i = 1, klon
[1421]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
[1403]482      END DO
483     END DO
484    END IF
485
[1409]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
[1191]506!======================================================================
507!     -- Calcul de l'effet de la couche limite --
508!======================================================================
509
510    IF (couchelimite) THEN             
[1212]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
[1191]520    END IF
521   
522    DO it=1, nbtr
[1409]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
[1191]525          CALL cltracrn(it, pdtphys, yu1, yv1,     &
526               cdragh, coefh,t_seri,ftsol,pctsrf,  &
527               tr_seri(:,:,it),trs(:,it),          &
[1403]528               paprs, pplay, zmasse * rg, &
[1191]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
[1403]545         
546
[1191]547!======================================================================
548!   Calcul de l'effet du puits radioactif
549!======================================================================
550    CALL radio_decay (radio,rnpb,pdtphys,tautr,tr_seri,d_tr_dec)
[1421]551
[1191]552    DO it=1,nbtr
553       IF(radio(it)) then     
554          DO k = 1, klev
555             DO i = 1, klon
556                tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_dec(i,k,it)
557             END DO
558          END DO
[4124]559          CALL minmaxqfi(tr_seri(:,:,it),0.,1.e33,'puits rn it='//TRIM(int2str(it)))
[1191]560       END IF
561    END DO
562
[1227]563!======================================================================
[1403]564!   Parameterization of ozone chemistry
565!======================================================================
566
567    IF (id_o3 /= 0) then
568       lmt_pas = NINT(86400./pdtphys)
569       IF (MOD(nstep - 1, lmt_pas) == 0) THEN
570          ! Once per day, update the coefficients for ozone chemistry:
571          CALL regr_pr_comb_coefoz(julien, xlat, paprs, pplay)
572       END IF
573       CALL o3_chem(julien, gmtime, t_seri, zmasse, pdtphys, xlat, &
574            xlon, tr_seri(:, :, id_o3))
575    END IF
576
[1191]577  END SUBROUTINE traclmdz
578
579
580  SUBROUTINE traclmdz_to_restart(trs_out)
581    ! This subroutine is called from phyredem.F where the module
582    ! variable trs is written to restart file (restartphy.nc)
583    USE dimphy
[4050]584    USE infotrac_phy, ONLY: nbtr
[1191]585   
586    REAL,DIMENSION(klon,nbtr), INTENT(OUT) :: trs_out
[1203]587    INTEGER :: ierr
[1212]588
589    IF ( ALLOCATED(trs) ) THEN
590       trs_out(:,:) = trs(:,:)
[1203]591    ELSE
[1212]592       ! No previous allocate of trs. This is the case for create_etat0_limit.
593       trs_out(:,:) = 0.0
[1203]594    END IF
[1212]595   
[1191]596  END SUBROUTINE traclmdz_to_restart
597 
598
599END MODULE traclmdz_mod
Note: See TracBrowser for help on using the repository browser.