source: LMDZ6/trunk/libf/phylmd/traclmdz_mod.f90 @ 5452

Last change on this file since 5452 was 5285, checked in by abarral, 2 months ago

As discussed internally, remove generic ONLY: ... for new _mod_h modules

  • 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
[4325]176       IF(.NOT.(tracers(iq)%isInPhysics)) CYCLE
[4056]177       it = it+1
[4393]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
[4393]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
[4389]204             WRITE(lunout,*)'prof.pb210 does not exist: use restart values'
[4046]205           END IF
[4393]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
[4393]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
[4393]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
[4393]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   
[5285]311    USE yomcst_mod_h
[5274]312USE dimphy
[4124]313    USE infotrac_phy, ONLY: nbtr, pbl_flg
314    USE strings_mod,  ONLY: int2str
[1403]315    USE regr_pr_comb_coefoz_m, ONLY: regr_pr_comb_coefoz
316    USE o3_chem_m, ONLY: o3_chem
[1785]317    USE indice_sol_mod
318
[1191]319
[5274]320
[1191]321!==========================================================================
322!                   -- DESCRIPTION DES ARGUMENTS --
323!==========================================================================
324
325! Input arguments
326!
327!Configuration grille,temps:
[1212]328    INTEGER,INTENT(IN) :: nstep      ! nombre d'appels de la physiq
[1403]329    INTEGER,INTENT(IN) :: julien     ! Jour julien
330    REAL,INTENT(IN)    :: gmtime
[1191]331    REAL,INTENT(IN)    :: pdtphys    ! Pas d'integration pour la physique (seconde) 
332    REAL,DIMENSION(klon),INTENT(IN) :: xlat    ! latitudes pour chaque point
[1403]333    REAL, INTENT(IN):: xlon(:) ! dim(klon) longitude
[1191]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)
[1403]341    REAL,intent(in):: zmasse (:, :)   ! dim(klon,klev) density of air, in kg/m2
[1191]342
343
344!Couche limite:
345!--------------
346!
[1742]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
[1191]351    LOGICAL,INTENT(IN)                   :: couchelimite
[1742]352    REAL,DIMENSION(klon,klev),INTENT(IN) :: sh         ! humidite specifique
[1670]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)
[1813]356    REAL,DIMENSION(klon),INTENT(IN)      :: wstar,ale_bl,ale_wake   ! wstar (m/s) and Avail. Lifti. Energ.
[1670]357    REAL,DIMENSION(klon),INTENT(IN)      :: zu10m   ! vent zonal 10m (m/s)
358    REAL,DIMENSION(klon),INTENT(IN)      :: zv10m   ! vent zonal 10m (m/s)
[1191]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
[1421]376    INTEGER :: lmt_pas ! number of time steps of "physics" per day
[1191]377
378    REAL,DIMENSION(klon)           :: d_trs    ! Td dans le reservoir
[1421]379    REAL,DIMENSION(klon,klev)      :: qsat     ! pression de la vapeur a saturation
[1191]380    REAL,DIMENSION(klon,klev,nbtr) :: d_tr_dec ! Td radioactive
[1421]381    REAL                           :: zrho     ! Masse Volumique de l'air KgA/m3
382    REAL                           :: amn, amx
[1191]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 
[1421]397
398!=================================================================
399! Update pseudo-vapor tracers
400!=================================================================
401
402    CALL q_sat(klon*klev,t_seri,pplay,qsat)
403
[1403]404    IF ( id_pcsat /= 0 ) THEN
405     DO k = 1, klev
406      DO i = 1, klon
[1421]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
[1403]412      END DO
413     END DO
414    END IF
[1191]415
[1403]416    IF ( id_pcocsat /= 0 ) THEN
417     DO k = 1, klev
418      DO i = 1, klon
[1421]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))
[1403]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
[1421]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
[1403]440      END DO
441     END DO
442    END IF
443
[1421]444
[1403]445    IF ( id_pcs0 /= 0 ) THEN
446     DO k = 1, klev
447      DO i = 1, klon
[1421]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
[1403]453      END DO
454     END DO
455    END IF
456
[1421]457
[1403]458    IF ( id_pcos0 /= 0 ) THEN
459     DO k = 1, klev
460      DO i = 1, klon
[1421]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
[1403]470      END DO
471     END DO
472    END IF
473
[1421]474
[1403]475    IF ( id_pcq0 /= 0 ) THEN
476     DO k = 1, klev
477      DO i = 1, klon
[1421]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
[1403]483      END DO
484     END DO
485    END IF
486
[1409]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
[1191]507!======================================================================
508!     -- Calcul de l'effet de la couche limite --
509!======================================================================
510
511    IF (couchelimite) THEN             
[1212]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
[1191]521    END IF
522   
523    DO it=1, nbtr
[1409]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
[1191]526          CALL cltracrn(it, pdtphys, yu1, yv1,     &
527               cdragh, coefh,t_seri,ftsol,pctsrf,  &
528               tr_seri(:,:,it),trs(:,it),          &
[1403]529               paprs, pplay, zmasse * rg, &
[1191]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
[1403]546         
547
[1191]548!======================================================================
549!   Calcul de l'effet du puits radioactif
550!======================================================================
551    CALL radio_decay (radio,rnpb,pdtphys,tautr,tr_seri,d_tr_dec)
[1421]552
[1191]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
[4124]560          CALL minmaxqfi(tr_seri(:,:,it),0.,1.e33,'puits rn it='//TRIM(int2str(it)))
[1191]561       END IF
562    END DO
563
[1227]564!======================================================================
[1403]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
[1191]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
[4050]585    USE infotrac_phy, ONLY: nbtr
[1191]586   
587    REAL,DIMENSION(klon,nbtr), INTENT(OUT) :: trs_out
[1203]588    INTEGER :: ierr
[1212]589
590    IF ( ALLOCATED(trs) ) THEN
591       trs_out(:,:) = trs(:,:)
[1203]592    ELSE
[1212]593       ! No previous allocate of trs. This is the case for create_etat0_limit.
594       trs_out(:,:) = 0.0
[1203]595    END IF
[1212]596   
[1191]597  END SUBROUTINE traclmdz_to_restart
598 
599
600END MODULE traclmdz_mod
Note: See TracBrowser for help on using the repository browser.