GCC Code Coverage Report
Directory: ./ Exec Total Coverage
File: phylmd/traclmdz_mod.F90 Lines: 103 202 51.0 %
Date: 2023-06-30 12:51:15 Branches: 140 358 39.1 %

Line Branch Exec Source
1
!$Id $
2
!
3
MODULE 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
62
CONTAINS
63
64
1
  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



2
    ALLOCATE( trs(klon,nbtr), stat=ierr)
79
1
    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

1991
    trs(:,:) = trs_in(:,:)
83
84
1
  END SUBROUTINE traclmdz_from_restart
85
86
87
1
  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
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, iq, i, k
116
2
    REAL, DIMENSION(klon_glo,klev) :: varglo ! variable temporaire sur la grille global
117
2
    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
2
      real     beryllium(klon,klev)
124
! profil initial Pb210
125
      integer ilesfil2
126
      parameter (ilesfil2=1)
127
      integer  irr2,kradio2
128
2
      real     plomb(klon,klev)
129
!! RomP <<<
130
! --------------------------------------------
131
! Allocation
132
! --------------------------------------------
133


1
    ALLOCATE( scavtr(nbtr), stat=ierr)
134
1
    IF (ierr /= 0) CALL abort_physic('traclmdz_init', 'pb in allocation 9',1)
135
3
    scavtr(:)=1.
136
137


1
    ALLOCATE( radio(nbtr), stat=ierr)
138
1
    IF (ierr /= 0) CALL abort_physic('traclmdz_init', 'pb in allocation 11',1)
139
3
    radio(:) = .false.    ! Par defaut pas decroissance radioactive
140
141



2
    ALLOCATE( masktr(klon,nbtr), stat=ierr)
142
1
    IF (ierr /= 0) CALL abort_physic('traclmdz_init', 'pb in allocation 2',1)
143
144



2
    ALLOCATE( fshtr(klon,nbtr), stat=ierr)
145
1
    IF (ierr /= 0) CALL abort_physic('traclmdz_init', 'pb in allocation 3',1)
146
147


1
    ALLOCATE( hsoltr(nbtr), stat=ierr)
148
1
    IF (ierr /= 0) CALL abort_physic('traclmdz_init', 'pb in allocation 4',1)
149
150


1
    ALLOCATE( tautr(nbtr), stat=ierr)
151
1
    IF (ierr /= 0) CALL abort_physic('traclmdz_init', 'pb in allocation 5',1)
152
3
    tautr(:)  = 0.
153
154


1
    ALLOCATE( vdeptr(nbtr), stat=ierr)
155
1
    IF (ierr /= 0) CALL abort_physic('traclmdz_init', 'pb in allocation 6',1)
156
3
    vdeptr(:) = 0.
157
158
159
1
    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
3
    aerosol(:) = .FALSE.  ! Tous les traceurs sont des gaz par defaut
168
169
!
170
! Recherche des traceurs connus : Be7, O3, CO2,...
171
! --------------------------------------------
172
1
    id_rn=0; id_pb=0; id_aga=0; id_be=0; id_o3=0
173
1
    id_pcsat=0; id_pcocsat=0; id_pcq=0; id_pcs0=0; id_pcos0=0; id_pcq0=0
174
    it = 0
175
6
    DO iq = 1, nqtot
176
5
       IF(.NOT.(tracers(iq)%isInPhysics)) CYCLE
177
2
       it = it+1
178
1
       SELECT CASE(tracers(iq)%name)
179
1
         CASE("RN");                  id_rn     = it ! radon
180
1
         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
190
         CASE DEFAULT
191



2
           WRITE(lunout,*) 'This is an unknown tracer in LMDZ : ', trim(tracers(iq)%name)
192
       END SELECT
193
194
1
       SELECT CASE(tracers(iq)%name)
195
         CASE("PB")                        !--- RomP >>> profil initial de PB210
196
1
           OPEN(ilesfil2,file='prof.pb210',status='old',iostat=irr2)
197
1
           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
204
1
             WRITE(lunout,*)'prof.pb210 does not exist: use restart values'
205
           END IF
206
         CASE("Aga","AGA")
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
218
         CASE("Be","BE","Be7","BE7")
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
235
         CASE("o3","O3")                    !--- Parametrisation par la chimie de Cariolle
236
           CALL alloc_coefoz                !--- Allocate ozone coefficients
237
           CALL press_coefoz                !--- Read input pressure levels
238
         CASE("pcs0","Pcs0", "pcos0","Pcos0", "pcq0","Pcq0")
239

2
           conv_flg(it)=0                   !--- No transport by convection for this tracer
240
       END SELECT
241
    END DO
242
243
!
244
! Valeurs specifiques pour les traceurs Rn222 et Pb210
245
! ----------------------------------------------
246

1
    IF ( id_rn/=0 .AND. id_pb/=0 ) THEN
247
1
       rnpb = .TRUE.
248
1
       radio(id_rn)= .TRUE.
249
1
       radio(id_pb)= .TRUE.
250
1
       pbl_flg(id_rn) = 0 ! au lieu de clsol=true ! CL au sol calcule
251
1
       pbl_flg(id_pb) = 0 ! au lieu de clsol=true
252
1
       aerosol(id_rn) = .FALSE.
253
1
       aerosol(id_pb) = .TRUE. ! le Pb est un aerosol
254
255
1
       CALL initrrnpb (ftsol,pctsrf,masktr,fshtr,hsoltr,tautr,vdeptr,scavtr)
256
    END IF
257
258
!
259
! Check if all tracers have restart values
260
! ----------------------------------------------
261
    it = 0
262
6
    DO iq = 1, nqtot
263

5
       IF(.NOT.(tracers(iq)%isAdvected .AND. tracers(iq)%isInPhysics)) CYCLE
264
2
       it = it+1
265
       ! Test if tracer is zero everywhere.
266
       ! Done by master process MPI and master thread OpenMP
267
2
       CALL gather(tr_seri(:,:,it),varglo)
268

2
       IF (is_mpi_root .AND. is_omp_root) THEN
269
2
          mintmp=MINVAL(varglo,dim=1)
270
2
          maxtmp=MAXVAL(varglo,dim=1)
271





164
          IF (MINVAL(mintmp,dim=1)==0. .AND. MAXVAL(maxtmp,dim=1)==0.) THEN
272
             ! Tracer is zero everywhere
273
             zero=.TRUE.
274
          ELSE
275
2
             zero=.FALSE.
276
          END IF
277
       END IF
278
279
       ! Distribute variable at all processes
280
2
       CALL bcast(zero)
281
282
       ! Initalize tracer that was not found in restart file.
283
3
       IF (zero) THEN
284
          ! The tracer was not found in restart file or it was equal zero everywhere.
285
          WRITE(lunout,*) "The tracer ",trim(tracers(iq)%name)," will be initialized"
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
301
       END IF
302
    END DO
303
304
1
  END SUBROUTINE traclmdz_init
305
306
288
  SUBROUTINE traclmdz(nstep, julien, gmtime, pdtphys, t_seri, paprs, pplay, &
307
864
       cdragh,  coefh, yu1, yv1, ftsol, pctsrf, xlat, xlon, couchelimite, sh, &
308
       rh, pphi, ustar, wstar, ale_bl, ale_wake,  zu10m, zv10m, &
309

288
       tr_seri, source, d_tr_cl,d_tr_dec, zmasse)               !RomP
310
311
    USE dimphy
312
    USE infotrac_phy, ONLY: nbtr, pbl_flg
313
    USE strings_mod,  ONLY: int2str
314
    USE regr_pr_comb_coefoz_m, ONLY: regr_pr_comb_coefoz
315
    USE o3_chem_m, ONLY: o3_chem
316
    USE indice_sol_mod
317
318
    INCLUDE "YOMCST.h"
319
320
!==========================================================================
321
!                   -- DESCRIPTION DES ARGUMENTS --
322
!==========================================================================
323
324
! Input arguments
325
!
326
!Configuration grille,temps:
327
    INTEGER,INTENT(IN) :: nstep      ! nombre d'appels de la physiq
328
    INTEGER,INTENT(IN) :: julien     ! Jour julien
329
    REAL,INTENT(IN)    :: gmtime
330
    REAL,INTENT(IN)    :: pdtphys    ! Pas d'integration pour la physique (seconde)
331
    REAL,DIMENSION(klon),INTENT(IN) :: xlat    ! latitudes pour chaque point
332
    REAL, INTENT(IN):: xlon(:) ! dim(klon) longitude
333
334
!
335
!Physique:
336
!--------
337
    REAL,DIMENSION(klon,klev),INTENT(IN)   :: t_seri  ! Temperature
338
    REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs   ! pression pour chaque inter-couche (en Pa)
339
    REAL,DIMENSION(klon,klev),INTENT(IN)   :: pplay   ! pression pour le mileu de chaque couche (en Pa)
340
    REAL,intent(in):: zmasse (:, :)   ! dim(klon,klev) density of air, in kg/m2
341
342
343
!Couche limite:
344
!--------------
345
!
346
    REAL,DIMENSION(klon),INTENT(IN)      :: cdragh     ! coeff drag pour T et Q
347
    REAL,DIMENSION(klon,klev),INTENT(IN) :: coefh      ! coeff melange CL (m**2/s)
348
    REAL,DIMENSION(klon),INTENT(IN)      :: yu1        ! vents au premier niveau
349
    REAL,DIMENSION(klon),INTENT(IN)      :: yv1        ! vents au premier niveau
350
    LOGICAL,INTENT(IN)                   :: couchelimite
351
    REAL,DIMENSION(klon,klev),INTENT(IN) :: sh         ! humidite specifique
352
    REAL,DIMENSION(klon,klev),INTENT(IN) :: rh      ! Humidite relative
353
    REAL,DIMENSION(klon,klev),INTENT(IN) :: pphi    ! geopotentie
354
    REAL,DIMENSION(klon),INTENT(IN)      :: ustar   ! ustar (m/s)
355
    REAL,DIMENSION(klon),INTENT(IN)      :: wstar,ale_bl,ale_wake   ! wstar (m/s) and Avail. Lifti. Energ.
356
    REAL,DIMENSION(klon),INTENT(IN)      :: zu10m   ! vent zonal 10m (m/s)
357
    REAL,DIMENSION(klon),INTENT(IN)      :: zv10m   ! vent zonal 10m (m/s)
358
359
! Arguments necessaires pour les sources et puits de traceur:
360
    REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: ftsol  ! Temperature du sol (surf)(Kelvin)
361
    REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: pctsrf ! Pourcentage de sol f(nature du sol)
362
363
! InOutput argument
364
    REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT)  :: tr_seri ! Concentration Traceur [U/KgA]
365
366
! Output argument
367
    REAL,DIMENSION(klon,nbtr), INTENT(OUT)        :: source  ! a voir lorsque le flux de surface est prescrit
368
    REAL,DIMENSION(klon,klev,nbtr), INTENT(OUT)   :: d_tr_cl ! Td couche limite/traceur
369
370
!=======================================================================================
371
!                        -- VARIABLES LOCALES TRACEURS --
372
!=======================================================================================
373
374
    INTEGER :: i, k, it
375
    INTEGER :: lmt_pas ! number of time steps of "physics" per day
376
377
576
    REAL,DIMENSION(klon)           :: d_trs    ! Td dans le reservoir
378
576
    REAL,DIMENSION(klon,klev)      :: qsat     ! pression de la vapeur a saturation
379
    REAL,DIMENSION(klon,klev,nbtr) :: d_tr_dec ! Td radioactive
380
    REAL                           :: zrho     ! Masse Volumique de l'air KgA/m3
381
    REAL                           :: amn, amx
382
!
383
!=================================================================
384
!  Ajout de la production en  Be7 (Beryllium) srcbe U/s/kgA
385
!=================================================================
386
!
387
288
    IF ( id_be /= 0 ) THEN
388
       DO k = 1, klev
389
          DO i = 1, klon
390
             tr_seri(i,k,id_be) = tr_seri(i,k,id_be)+srcbe(i,k)*pdtphys
391
          END DO
392
       END DO
393
       WRITE(*,*) 'Ajout srcBe dans tr_seri: OK'
394
    END IF
395
396
397
!=================================================================
398
! Update pseudo-vapor tracers
399
!=================================================================
400
401
288
    CALL q_sat(klon*klev,t_seri,pplay,qsat)
402
403
288
    IF ( id_pcsat /= 0 ) THEN
404
     DO k = 1, klev
405
      DO i = 1, klon
406
         IF ( pplay(i,k).GE.85000.) THEN
407
            tr_seri(i,k,id_pcsat) = qsat(i,k)
408
         ELSE
409
            tr_seri(i,k,id_pcsat) = MIN (qsat(i,k), tr_seri(i,k,id_pcsat))
410
         END IF
411
      END DO
412
     END DO
413
    END IF
414
415
288
    IF ( id_pcocsat /= 0 ) THEN
416
     DO k = 1, klev
417
      DO i = 1, klon
418
         IF ( pplay(i,k).GE.85000.) THEN
419
            IF ( pctsrf (i, is_oce) > 0. ) THEN
420
               tr_seri(i,k,id_pcocsat) = qsat(i,k)
421
            ELSE
422
               tr_seri(i,k,id_pcocsat) = 0.
423
          END IF
424
       ELSE
425
          tr_seri(i,k,id_pcocsat) = MIN (qsat(i,k), tr_seri(i,k,id_pcocsat))
426
       END IF
427
      END DO
428
     END DO
429
    END IF
430
431
288
    IF ( id_pcq /= 0 ) THEN
432
     DO k = 1, klev
433
      DO i = 1, klon
434
         IF ( pplay(i,k).GE.85000.) THEN
435
            tr_seri(i,k,id_pcq) = sh(i,k)
436
         ELSE
437
            tr_seri(i,k,id_pcq) = MIN (qsat(i,k), tr_seri(i,k,id_pcq))
438
         END IF
439
      END DO
440
     END DO
441
    END IF
442
443
444
288
    IF ( id_pcs0 /= 0 ) THEN
445
     DO k = 1, klev
446
      DO i = 1, klon
447
         IF ( pplay(i,k).GE.85000.) THEN
448
            tr_seri(i,k,id_pcs0) = qsat(i,k)
449
         ELSE
450
            tr_seri(i,k,id_pcs0) = MIN (qsat(i,k), tr_seri(i,k,id_pcs0))
451
         END IF
452
      END DO
453
     END DO
454
    END IF
455
456
457
288
    IF ( id_pcos0 /= 0 ) THEN
458
     DO k = 1, klev
459
      DO i = 1, klon
460
         IF ( pplay(i,k).GE.85000.) THEN
461
            IF ( pctsrf (i, is_oce) > 0. ) THEN
462
               tr_seri(i,k,id_pcos0) = qsat(i,k)
463
            ELSE
464
               tr_seri(i,k,id_pcos0) = 0.
465
            END IF
466
         ELSE
467
            tr_seri(i,k,id_pcos0) = MIN (qsat(i,k), tr_seri(i,k,id_pcos0))
468
         END IF
469
      END DO
470
     END DO
471
    END IF
472
473
474
288
    IF ( id_pcq0 /= 0 ) THEN
475
     DO k = 1, klev
476
      DO i = 1, klon
477
         IF ( pplay(i,k).GE.85000.) THEN
478
            tr_seri(i,k,id_pcq0) = sh(i,k)
479
         ELSE
480
            tr_seri(i,k,id_pcq0) = MIN (qsat(i,k), tr_seri(i,k,id_pcq0))
481
         END IF
482
      END DO
483
     END DO
484
    END IF
485
486
!=================================================================
487
! Update tracer : Age of stratospheric air
488
!=================================================================
489
288
    IF (id_aga/=0) THEN
490
491
       ! Bottom layers
492
       DO k = 1, lev_1p5km
493
          tr_seri(:,k,id_aga) = 0.0
494
       END DO
495
496
       ! Layers above 1.5km
497
       DO k = lev_1p5km+1,klev-1
498
          tr_seri(:,k,id_aga) = tr_seri(:,k,id_aga) + pdtphys
499
       END DO
500
501
       ! Top layer
502
       tr_seri(:,klev,id_aga) = tr_seri(:,klev-1,id_aga)
503
504
    END IF
505
506
!======================================================================
507
!     -- Calcul de l'effet de la couche limite --
508
!======================================================================
509
510
288
    IF (couchelimite) THEN
511

573408
       source(:,:) = 0.0
512
513
288
       IF (id_be /=0) THEN
514
          DO i=1, klon
515
             zrho = pplay(i,1)/t_seri(i,1)/RD
516
             source(i,id_be) = - vdeptr(id_be)*tr_seri(i,1,id_be)*zrho
517
          END DO
518
       END IF
519
520
    END IF
521
522
864
    DO it=1, nbtr
523


864
       IF (couchelimite .AND. pbl_flg(it) == 0 .AND. (it==id_rn .OR. it==id_pb)) THEN
524
          ! couche limite avec quantite dans le sol calculee
525
          CALL cltracrn(it, pdtphys, yu1, yv1,     &
526
               cdragh, coefh,t_seri,ftsol,pctsrf,  &
527
               tr_seri(:,:,it),trs(:,it),          &
528
               paprs, pplay, zmasse * rg, &
529
               masktr(:,it),fshtr(:,it),hsoltr(it),&
530
               tautr(it),vdeptr(it),               &
531

22352256
               xlat,d_tr_cl(:,:,it),d_trs)
532
533
23040
          DO k = 1, klev
534
22352256
             DO i = 1, klon
535
22351680
                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
573120
          DO i = 1, klon
541
573120
             trs(i,it) = trs(i,it) + d_trs(i)
542
          END DO
543
       END IF
544
    END DO
545
546
547
!======================================================================
548
!   Calcul de l'effet du puits radioactif
549
!======================================================================
550
288
    CALL radio_decay (radio,rnpb,pdtphys,tautr,tr_seri,d_tr_dec)
551
552
864
    DO it=1,nbtr
553
864
       IF(radio(it)) then
554
23040
          DO k = 1, klev
555
22352256
             DO i = 1, klon
556
22351680
                tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_dec(i,k,it)
557
             END DO
558
          END DO
559
576
          CALL minmaxqfi(tr_seri(:,:,it),0.,1.e33,'puits rn it='//TRIM(int2str(it)))
560
       END IF
561
    END DO
562
563
!======================================================================
564
!   Parameterization of ozone chemistry
565
!======================================================================
566
567
288
    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
577
288
  END SUBROUTINE traclmdz
578
579
580
2
  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
584
    USE infotrac_phy, ONLY: nbtr
585
586
    REAL,DIMENSION(klon,nbtr), INTENT(OUT) :: trs_out
587
    INTEGER :: ierr
588
589
2
    IF ( ALLOCATED(trs) ) THEN
590

3982
       trs_out(:,:) = trs(:,:)
591
    ELSE
592
       ! No previous allocate of trs. This is the case for create_etat0_limit.
593
       trs_out(:,:) = 0.0
594
    END IF
595
596
2
  END SUBROUTINE traclmdz_to_restart
597
598
599
END MODULE traclmdz_mod