source: LMDZ6/branches/Amaury_dev/libf/phylmd/traclmdz_mod.F90 @ 5136

Last change on this file since 5136 was 5117, checked in by abarral, 4 months ago

rename modules properly lmdz_*
move some unused files to obsolete/
(lint) uppercase fortran keywords

  • 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  ! In this module all tracers specific to LMDZ are treated. This module is used
6  ! only if running without any other chemestry model as INCA or REPROBUS.
7  USE lmdz_abort_physic, ONLY: abort_physic
8  IMPLICIT NONE
9
10  REAL, DIMENSION(:, :), ALLOCATABLE, SAVE :: masktr   ! Masque reservoir de sol traceur
11  !$OMP THREADPRIVATE(masktr)                        ! Masque de l'echange avec la surface (1 = reservoir) ou (possible >= 1 )
12  REAL, DIMENSION(:, :), ALLOCATABLE, SAVE :: fshtr    ! Flux surfacique dans le reservoir de sol
13  !$OMP THREADPRIVATE(fshtr)
14  REAL, DIMENSION(:), ALLOCATABLE, SAVE :: hsoltr   ! Epaisseur equivalente du reservoir de sol
15  !$OMP THREADPRIVATE(hsoltr)
16
17  !Radioelements:
18  !--------------
19
20  REAL, DIMENSION(:), ALLOCATABLE, SAVE :: tautr    ! Constante de decroissance radioactive
21  !$OMP THREADPRIVATE(tautr)
22  REAL, DIMENSION(:), ALLOCATABLE, SAVE :: vdeptr   ! Vitesse de depot sec dans la couche Brownienne
23  !$OMP THREADPRIVATE(vdeptr)
24  REAL, DIMENSION(:), ALLOCATABLE, SAVE :: scavtr   ! Coefficient de lessivage
25  !$OMP THREADPRIVATE(scavtr)
26  REAL, DIMENSION(:, :), ALLOCATABLE, SAVE :: srcbe    ! Production du beryllium7 dans l atmosphere (U/s/kgA)
27  !$OMP THREADPRIVATE(srcbe)
28
29  LOGICAL, DIMENSION(:), ALLOCATABLE, SAVE :: radio    ! radio(it)   = true  => decroisssance radioactive
30  !$OMP THREADPRIVATE(radio)
31
32  REAL, DIMENSION(:, :), ALLOCATABLE, SAVE :: trs     ! Conc. radon ds le sol
33  !$OMP THREADPRIVATE(trs)
34
35  INTEGER, SAVE :: id_aga      ! Identification number for tracer : Age of stratospheric air
36  !$OMP THREADPRIVATE(id_aga)
37  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.
38  !$OMP THREADPRIVATE(lev_1p5km)
39
40  INTEGER, SAVE :: id_rn, id_pb ! Identification number for tracer : radon (Rn222), lead (Pb210)
41  !$OMP THREADPRIVATE(id_rn, id_pb)
42
43  INTEGER, SAVE :: id_be       ! Activation et position du traceur Be7 [ id_be=0 -> desactive ]
44  !$OMP THREADPRIVATE(id_be)
45
46  INTEGER, SAVE :: id_pcsat, id_pcocsat, id_pcq ! traceurs pseudo-vapeur CL qsat, qsat_oc, q
47  !$OMP THREADPRIVATE(id_pcsat, id_pcocsat, id_pcq)
48  INTEGER, SAVE :: id_pcs0, id_pcos0, id_pcq0   ! traceurs pseudo-vapeur CL qsat, qsat_oc, q
49  !                                              ! qui ne sont pas transportes par la convection
50  !$OMP THREADPRIVATE(id_pcs0, id_pcos0, id_pcq0)
51
52  INTEGER, SAVE :: id_o3
53  !$OMP THREADPRIVATE(id_o3)
54  ! index of ozone tracer with Cariolle parameterization
55  ! 0 means no ozone tracer
56
57  LOGICAL, SAVE :: rnpb = .FALSE. ! Presence du couple Rn222, Pb210
58  !$OMP THREADPRIVATE(rnpb)
59
60
61CONTAINS
62
63  SUBROUTINE traclmdz_from_restart(trs_in)
64    ! This SUBROUTINE initialize the module saved variable trs with values from restart file (startphy.nc).
65    ! This SUBROUTINE is called from phyetat0 after the field trs_in has been read.
66
67    USE dimphy
68    USE infotrac_phy, ONLY: nbtr
69
70    ! Input argument
71    REAL, DIMENSION(klon, nbtr), INTENT(IN) :: trs_in
72
73    ! Local variables
74    INTEGER :: ierr
75
76    ! Allocate restart variables trs
77    ALLOCATE(trs(klon, nbtr), stat = ierr)
78    IF (ierr /= 0) CALL abort_physic('traclmdz_from_restart', 'pb in allocation 1', 1)
79
80    ! Initialize trs with values read from restart file
81    trs(:, :) = trs_in(:, :)
82
83  END SUBROUTINE traclmdz_from_restart
84
85
86  SUBROUTINE traclmdz_init(pctsrf, xlat, xlon, ftsol, tr_seri, t_seri, pplay, sh, pdtphys, aerosol, lessivage)
87    ! This SUBROUTINE allocates and initialize module variables and control variables.
88    ! Initialization of the tracers should be done here only for those not found in the restart file.
89    USE dimphy
90    USE infotrac_phy, ONLY: nbtr, nqtot, tracers, pbl_flg, conv_flg
91    USE regr_pr_comb_coefoz_m, ONLY: alloc_coefoz
92    USE press_coefoz_m, ONLY: press_coefoz
93    USE lmdz_grid_phy
94    USE lmdz_phys_para
95    USE indice_sol_mod
96    USE lmdz_print_control, ONLY: lunout
97
98    ! Input variables
99    REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: pctsrf ! Pourcentage de sol f(nature du sol)
100    REAL, DIMENSION(klon), INTENT(IN) :: xlat   ! latitudes en degres pour chaque point
101    REAL, DIMENSION(klon), INTENT(IN) :: xlon   ! longitudes en degres pour chaque point
102    REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: ftsol  ! Temperature du sol (surf)(Kelvin)
103    REAL, DIMENSION(klon, klev, nbtr), INTENT(INOUT) :: tr_seri! Concentration Traceur [U/KgA]
104    REAL, DIMENSION(klon, klev), INTENT(IN) :: t_seri  ! Temperature
105    REAL, DIMENSION(klon, klev), INTENT(IN) :: pplay   ! pression pour le mileu de chaque couche (en Pa)
106    REAL, DIMENSION(klon, klev), INTENT(IN) :: sh      ! humidite specifique
107    REAL, INTENT(IN) :: pdtphys ! Pas d'integration pour la physique (seconde)
108
109    ! Output variables
110    LOGICAL, DIMENSION(nbtr), INTENT(OUT) :: aerosol
111    LOGICAL, INTENT(OUT) :: lessivage
112
113    ! Local variables
114    INTEGER :: ierr, it, iq, i, k
115    REAL, DIMENSION(klon_glo, klev) :: varglo ! variable temporaire sur la grille global
116    REAL, DIMENSION(klev) :: mintmp, maxtmp
117    LOGICAL :: zero
118    ! RomP >>> profil initial Be7
119    INTEGER ilesfil
120    parameter (ilesfil = 1)
121    INTEGER  irr, kradio
122    real     beryllium(klon, klev)
123    ! profil initial Pb210
124    INTEGER ilesfil2
125    parameter (ilesfil2 = 1)
126    INTEGER  irr2, kradio2
127    real     plomb(klon, klev)
128    !! RomP <<<
129    ! --------------------------------------------
130    ! Allocation
131    ! --------------------------------------------
132    ALLOCATE(scavtr(nbtr), stat = ierr)
133    IF (ierr /= 0) CALL abort_physic('traclmdz_init', 'pb in allocation 9', 1)
134    scavtr(:) = 1.
135
136    ALLOCATE(radio(nbtr), stat = ierr)
137    IF (ierr /= 0) CALL abort_physic('traclmdz_init', 'pb in allocation 11', 1)
138    radio(:) = .FALSE.    ! Par defaut pas decroissance radioactive
139
140    ALLOCATE(masktr(klon, nbtr), stat = ierr)
141    IF (ierr /= 0) CALL abort_physic('traclmdz_init', 'pb in allocation 2', 1)
142
143    ALLOCATE(fshtr(klon, nbtr), stat = ierr)
144    IF (ierr /= 0) CALL abort_physic('traclmdz_init', 'pb in allocation 3', 1)
145
146    ALLOCATE(hsoltr(nbtr), stat = ierr)
147    IF (ierr /= 0) CALL abort_physic('traclmdz_init', 'pb in allocation 4', 1)
148
149    ALLOCATE(tautr(nbtr), stat = ierr)
150    IF (ierr /= 0) CALL abort_physic('traclmdz_init', 'pb in allocation 5', 1)
151    tautr(:) = 0.
152
153    ALLOCATE(vdeptr(nbtr), stat = ierr)
154    IF (ierr /= 0) CALL abort_physic('traclmdz_init', 'pb in allocation 6', 1)
155    vdeptr(:) = 0.
156
157    lessivage = .TRUE.
158    !!jyg(20130206) : le choix d activation du lessivage est fait dans phytrac avec iflag_lscav
159    !!    CALL getin('lessivage',lessivage)
160    !!    IF(lessivage) THEN
161    !!     PRINT*,'lessivage lsc ON'
162    !!    else
163    !!     PRINT*,'lessivage lsc OFF'
164    !!    endif
165    aerosol(:) = .FALSE.  ! Tous les traceurs sont des gaz par defaut
166
167    ! Recherche des traceurs connus : Be7, O3, CO2,...
168    ! --------------------------------------------
169    id_rn = 0; id_pb = 0; id_aga = 0; id_be = 0; id_o3 = 0
170    id_pcsat = 0; id_pcocsat = 0; id_pcq = 0; id_pcs0 = 0; id_pcos0 = 0; id_pcq0 = 0
171    it = 0
172    DO iq = 1, nqtot
173      IF(.NOT.(tracers(iq)%isInPhysics)) CYCLE
174      it = it + 1
175      SELECT CASE(tracers(iq)%name)
176      CASE("RN");                  id_rn = it ! radon
177      CASE("PB");                  id_pb = it ! plomb
178      CASE("Aga", "AGA");           id_aga = it ! Age of stratospheric air
179      CASE("Be", "BE", "Be7", "BE7"); id_be = it ! Recherche du Beryllium 7
180      CASE("o3", "O3");             id_o3 = it ! Recherche de l'ozone
181      CASE("pcsat", "Pcsat");     id_pcsat = it
182      CASE("pcocsat", "Pcocsat");   id_pcocsat = it
183      CASE("pcq", "Pcq");       id_pcq = it
184      CASE("pcs0", "Pcs0");      id_pcs0 = it
185      CASE("pcos0", "Pcos0");     id_pcos0 = it
186      CASE("pcq0", "Pcq0");      id_pcq0 = it
187      CASE DEFAULT
188        WRITE(lunout, *) 'This is an unknown tracer in LMDZ : ', trim(tracers(iq)%name)
189      END SELECT
190
191      SELECT CASE(tracers(iq)%name)
192      CASE("PB")                        !--- RomP >>> profil initial de PB210
193        OPEN(ilesfil2, file = 'prof.pb210', status = 'old', iostat = irr2)
194        IF(irr2 == 0) THEN
195          READ(ilesfil2, *) kradio2
196          WRITE(lunout, *)'number of levels for pb210 profile ', kradio2
197          DO k = kradio2, 1, -1; READ (ilesfil2, *) plomb(:, k);
198          END DO
199          CLOSE(ilesfil2)
200          tr_seri(:, :, id_pb) = plomb(:, :)
201        ELSE
202          WRITE(lunout, *)'prof.pb210 does not exist: use restart values'
203        END IF
204      CASE("Aga", "AGA")
205        radio(id_aga) = .FALSE.
206        aerosol(id_aga) = .FALSE.
207        pbl_flg(id_aga) = 0
208        ! Find the first model layer above 1.5km from the surface
209        IF (klev>=30) THEN
210          lev_1p5km = 6                  !--- NB: This value is for klev=39
211        ELSE IF (klev>=10) THEN
212          lev_1p5km = 5                  !--- NB: This value is for klev=19
213        ELSE
214          lev_1p5km = klev / 2
215        END IF
216      CASE("Be", "BE", "Be7", "BE7")
217        ALLOCATE(srcbe(klon, klev))
218        radio(id_be) = .TRUE.
219        aerosol(id_be) = .TRUE.         !--- Le Be est un aerosol
220        CALL init_be(pctsrf, pplay, masktr(:, id_be), tautr(id_be), vdeptr(id_be), scavtr(id_be), srcbe)
221        WRITE(lunout, *) 'Initialisation srcBe: OK'
222        !--- RomP >>> profil initial de Be7
223        OPEN(ilesfil, file = 'prof.be7', status = 'old', iostat = irr)
224        IF(irr == 0) THEN
225          READ(ilesfil, *) kradio
226          WRITE(lunout, *)'number of levels for Be7 profile ', kradio
227          DO k = kradio, 1, -1; READ(ilesfil, *) beryllium(:, k);
228          END DO
229          CLOSE(ilesfil)
230          tr_seri(:, :, id_be) = beryllium(:, :)
231        ELSE
232          WRITE(lunout, *)'Prof. Be7 does not exist: use restart values'
233        END IF
234      CASE("o3", "O3")                    !--- Parametrisation par la chimie de Cariolle
235        CALL alloc_coefoz                !--- Allocate ozone coefficients
236        CALL press_coefoz                !--- Read input pressure levels
237      CASE("pcs0", "Pcs0", "pcos0", "Pcos0", "pcq0", "Pcq0")
238        conv_flg(it) = 0                   !--- No transport by convection for this tracer
239      END SELECT
240    END DO
241
242    ! Valeurs specifiques pour les traceurs Rn222 et Pb210
243    ! ----------------------------------------------
244    IF (id_rn/=0 .AND. id_pb/=0) THEN
245      rnpb = .TRUE.
246      radio(id_rn) = .TRUE.
247      radio(id_pb) = .TRUE.
248      pbl_flg(id_rn) = 0 ! au lieu de clsol=true ! CL au sol calcule
249      pbl_flg(id_pb) = 0 ! au lieu de clsol=true
250      aerosol(id_rn) = .FALSE.
251      aerosol(id_pb) = .TRUE. ! le Pb est un aerosol
252
253      CALL initrrnpb (ftsol, pctsrf, masktr, fshtr, hsoltr, tautr, vdeptr, scavtr)
254    END IF
255
256    ! Check if all tracers have restart values
257    ! ----------------------------------------------
258    it = 0
259    DO iq = 1, nqtot
260      IF(.NOT.(tracers(iq)%isAdvected .AND. tracers(iq)%isInPhysics)) CYCLE
261      it = it + 1
262      ! Test if tracer is zero everywhere.
263      ! Done by master process MPI and master thread OpenMP
264      CALL gather(tr_seri(:, :, it), varglo)
265      IF (is_mpi_root .AND. is_omp_root) THEN
266        mintmp = MINVAL(varglo, dim = 1)
267        maxtmp = MAXVAL(varglo, dim = 1)
268        IF (MINVAL(mintmp, dim = 1)==0. .AND. MAXVAL(maxtmp, dim = 1)==0.) THEN
269          ! Tracer is zero everywhere
270          zero = .TRUE.
271        ELSE
272          zero = .FALSE.
273        END IF
274      END IF
275
276      ! Distribute variable at all processes
277      CALL bcast(zero)
278
279      ! Initalize tracer that was not found in restart file.
280      IF (zero) THEN
281        ! The tracer was not found in restart file or it was equal zero everywhere.
282        WRITE(lunout, *) "The tracer ", trim(tracers(iq)%name), " will be initialized"
283        IF (it==id_pcsat .OR. it==id_pcq .OR. &
284                it==id_pcs0 .OR. it==id_pcq0) THEN
285          tr_seri(:, :, it) = 100.
286        ELSE IF (it==id_pcocsat .OR. it==id_pcos0) THEN
287          DO i = 1, klon
288            IF (pctsrf (i, is_oce) == 0.) THEN
289              tr_seri(i, :, it) = 0.
290            ELSE
291              tr_seri(i, :, it) = 100.
292            END IF
293          END DO
294        ELSE
295          ! No specific initialization exist for this tracer
296          tr_seri(:, :, it) = 0.
297        END IF
298      END IF
299    END DO
300
301  END SUBROUTINE traclmdz_init
302
303  SUBROUTINE traclmdz(nstep, julien, gmtime, pdtphys, t_seri, paprs, pplay, &
304          cdragh, coefh, yu1, yv1, ftsol, pctsrf, xlat, xlon, couchelimite, sh, &
305          rh, pphi, ustar, wstar, ale_bl, ale_wake, zu10m, zv10m, &
306          tr_seri, source, d_tr_cl, d_tr_dec, zmasse)               !RomP
307
308    USE dimphy
309    USE infotrac_phy, ONLY: nbtr, pbl_flg
310    USE lmdz_strings, ONLY: int2str
311    USE regr_pr_comb_coefoz_m, ONLY: regr_pr_comb_coefoz
312    USE o3_chem_m, ONLY: o3_chem
313    USE indice_sol_mod
314    USE lmdz_q_sat, ONLY: q_sat
315
316    INCLUDE "YOMCST.h"
317
318    !==========================================================================
319    !                   -- DESCRIPTION DES ARGUMENTS --
320    !==========================================================================
321
322    ! Input arguments
323
324    !Configuration grille,temps:
325    INTEGER, INTENT(IN) :: nstep      ! nombre d'appels de la physiq
326    INTEGER, INTENT(IN) :: julien     ! Jour julien
327    REAL, INTENT(IN) :: gmtime
328    REAL, INTENT(IN) :: pdtphys    ! Pas d'integration pour la physique (seconde)
329    REAL, DIMENSION(klon), INTENT(IN) :: xlat    ! latitudes pour chaque point
330    REAL, INTENT(IN) :: xlon(:) ! dim(klon) longitude
331
332    !Physique:
333    !--------
334    REAL, DIMENSION(klon, klev), INTENT(IN) :: t_seri  ! Temperature
335    REAL, DIMENSION(klon, klev + 1), INTENT(IN) :: paprs   ! pression pour chaque inter-couche (en Pa)
336    REAL, DIMENSION(klon, klev), INTENT(IN) :: pplay   ! pression pour le mileu de chaque couche (en Pa)
337    REAL, INTENT(IN) :: zmasse (:, :)   ! dim(klon,klev) density of air, in kg/m2
338
339
340    !Couche limite:
341    !--------------
342
343    REAL, DIMENSION(klon), INTENT(IN) :: cdragh     ! coeff drag pour T et Q
344    REAL, DIMENSION(klon, klev), INTENT(IN) :: coefh      ! coeff melange CL (m**2/s)
345    REAL, DIMENSION(klon), INTENT(IN) :: yu1        ! vents au premier niveau
346    REAL, DIMENSION(klon), INTENT(IN) :: yv1        ! vents au premier niveau
347    LOGICAL, INTENT(IN) :: couchelimite
348    REAL, DIMENSION(klon, klev), INTENT(IN) :: sh         ! humidite specifique
349    REAL, DIMENSION(klon, klev), INTENT(IN) :: rh      ! Humidite relative
350    REAL, DIMENSION(klon, klev), INTENT(IN) :: pphi    ! geopotentie
351    REAL, DIMENSION(klon), INTENT(IN) :: ustar   ! ustar (m/s)
352    REAL, DIMENSION(klon), INTENT(IN) :: wstar, ale_bl, ale_wake   ! wstar (m/s) and Avail. Lifti. Energ.
353    REAL, DIMENSION(klon), INTENT(IN) :: zu10m   ! vent zonal 10m (m/s)
354    REAL, DIMENSION(klon), INTENT(IN) :: zv10m   ! vent zonal 10m (m/s)
355
356    ! Arguments necessaires pour les sources et puits de traceur:
357    REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: ftsol  ! Temperature du sol (surf)(Kelvin)
358    REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: pctsrf ! Pourcentage de sol f(nature du sol)
359
360    ! InOutput argument
361    REAL, DIMENSION(klon, klev, nbtr), INTENT(INOUT) :: tr_seri ! Concentration Traceur [U/KgA]
362
363    ! Output argument
364    REAL, DIMENSION(klon, nbtr), INTENT(OUT) :: source  ! a voir lorsque le flux de surface est prescrit
365    REAL, DIMENSION(klon, klev, nbtr), INTENT(OUT) :: d_tr_cl ! Td couche limite/traceur
366
367    !=======================================================================================
368    !                        -- VARIABLES LOCALES TRACEURS --
369    !=======================================================================================
370
371    INTEGER :: i, k, it
372    INTEGER :: lmt_pas ! number of time steps of "physics" per day
373
374    REAL, DIMENSION(klon) :: d_trs    ! Td dans le reservoir
375    REAL, DIMENSION(klon, klev) :: qsat     ! pression de la vapeur a saturation
376    REAL, DIMENSION(klon, klev, nbtr) :: d_tr_dec ! Td radioactive
377    REAL :: zrho     ! Masse Volumique de l'air KgA/m3
378    REAL :: amn, amx
379
380    !=================================================================
381    !  Ajout de la production en  Be7 (Beryllium) srcbe U/s/kgA
382    !=================================================================
383
384    IF (id_be /= 0) THEN
385      DO k = 1, klev
386        DO i = 1, klon
387          tr_seri(i, k, id_be) = tr_seri(i, k, id_be) + srcbe(i, k) * pdtphys
388        END DO
389      END DO
390      WRITE(*, *) 'Ajout srcBe dans tr_seri: OK'
391    END IF
392
393
394    !=================================================================
395    ! Update pseudo-vapor tracers
396    !=================================================================
397
398    CALL q_sat(klon * klev, t_seri, pplay, qsat)
399
400    IF (id_pcsat /= 0) THEN
401      DO k = 1, klev
402        DO i = 1, klon
403          IF (pplay(i, k)>=85000.) THEN
404            tr_seri(i, k, id_pcsat) = qsat(i, k)
405          ELSE
406            tr_seri(i, k, id_pcsat) = MIN (qsat(i, k), tr_seri(i, k, id_pcsat))
407          END IF
408        END DO
409      END DO
410    END IF
411
412    IF (id_pcocsat /= 0) THEN
413      DO k = 1, klev
414        DO i = 1, klon
415          IF (pplay(i, k)>=85000.) THEN
416            IF (pctsrf (i, is_oce) > 0.) THEN
417              tr_seri(i, k, id_pcocsat) = qsat(i, k)
418            ELSE
419              tr_seri(i, k, id_pcocsat) = 0.
420            END IF
421          ELSE
422            tr_seri(i, k, id_pcocsat) = MIN (qsat(i, k), tr_seri(i, k, id_pcocsat))
423          END IF
424        END DO
425      END DO
426    END IF
427
428    IF (id_pcq /= 0) THEN
429      DO k = 1, klev
430        DO i = 1, klon
431          IF (pplay(i, k)>=85000.) THEN
432            tr_seri(i, k, id_pcq) = sh(i, k)
433          ELSE
434            tr_seri(i, k, id_pcq) = MIN (qsat(i, k), tr_seri(i, k, id_pcq))
435          END IF
436        END DO
437      END DO
438    END IF
439
440    IF (id_pcs0 /= 0) THEN
441      DO k = 1, klev
442        DO i = 1, klon
443          IF (pplay(i, k)>=85000.) THEN
444            tr_seri(i, k, id_pcs0) = qsat(i, k)
445          ELSE
446            tr_seri(i, k, id_pcs0) = MIN (qsat(i, k), tr_seri(i, k, id_pcs0))
447          END IF
448        END DO
449      END DO
450    END IF
451
452    IF (id_pcos0 /= 0) THEN
453      DO k = 1, klev
454        DO i = 1, klon
455          IF (pplay(i, k)>=85000.) THEN
456            IF (pctsrf (i, is_oce) > 0.) THEN
457              tr_seri(i, k, id_pcos0) = qsat(i, k)
458            ELSE
459              tr_seri(i, k, id_pcos0) = 0.
460            END IF
461          ELSE
462            tr_seri(i, k, id_pcos0) = MIN (qsat(i, k), tr_seri(i, k, id_pcos0))
463          END IF
464        END DO
465      END DO
466    END IF
467
468    IF (id_pcq0 /= 0) THEN
469      DO k = 1, klev
470        DO i = 1, klon
471          IF (pplay(i, k)>=85000.) THEN
472            tr_seri(i, k, id_pcq0) = sh(i, k)
473          ELSE
474            tr_seri(i, k, id_pcq0) = MIN (qsat(i, k), tr_seri(i, k, id_pcq0))
475          END IF
476        END DO
477      END DO
478    END IF
479
480    !=================================================================
481    ! Update tracer : Age of stratospheric air
482    !=================================================================
483    IF (id_aga/=0) THEN
484
485      ! Bottom layers
486      DO k = 1, lev_1p5km
487        tr_seri(:, k, id_aga) = 0.0
488      END DO
489
490      ! Layers above 1.5km
491      DO k = lev_1p5km + 1, klev - 1
492        tr_seri(:, k, id_aga) = tr_seri(:, k, id_aga) + pdtphys
493      END DO
494
495      ! Top layer
496      tr_seri(:, klev, id_aga) = tr_seri(:, klev - 1, id_aga)
497
498    END IF
499
500    !======================================================================
501    !     -- Calcul de l'effet de la couche limite --
502    !======================================================================
503
504    IF (couchelimite) THEN
505      source(:, :) = 0.0
506
507      IF (id_be /=0) THEN
508        DO i = 1, klon
509          zrho = pplay(i, 1) / t_seri(i, 1) / RD
510          source(i, id_be) = - vdeptr(id_be) * tr_seri(i, 1, id_be) * zrho
511        END DO
512      END IF
513
514    END IF
515
516    DO it = 1, nbtr
517      IF (couchelimite .AND. pbl_flg(it) == 0 .AND. (it==id_rn .OR. it==id_pb)) THEN
518        ! couche limite avec quantite dans le sol calculee
519        CALL cltracrn(it, pdtphys, yu1, yv1, &
520                cdragh, coefh, t_seri, ftsol, pctsrf, &
521                tr_seri(:, :, it), trs(:, it), &
522                paprs, pplay, zmasse * rg, &
523                masktr(:, it), fshtr(:, it), hsoltr(it), &
524                tautr(it), vdeptr(it), &
525                xlat, d_tr_cl(:, :, it), d_trs)
526
527        DO k = 1, klev
528          DO i = 1, klon
529            tr_seri(i, k, it) = tr_seri(i, k, it) + d_tr_cl(i, k, it)
530          END DO
531        END DO
532
533        ! Traceur dans le reservoir sol
534        DO i = 1, klon
535          trs(i, it) = trs(i, it) + d_trs(i)
536        END DO
537      END IF
538    END DO
539
540
541    !======================================================================
542    !   Calcul de l'effet du puits radioactif
543    !======================================================================
544    CALL radio_decay (radio, rnpb, pdtphys, tautr, tr_seri, d_tr_dec)
545
546    DO it = 1, nbtr
547      IF(radio(it)) THEN
548        DO k = 1, klev
549          DO i = 1, klon
550            tr_seri(i, k, it) = tr_seri(i, k, it) + d_tr_dec(i, k, it)
551          END DO
552        END DO
553        CALL minmaxqfi(tr_seri(:, :, it), 0., 1.e33, 'puits rn it=' // TRIM(int2str(it)))
554      END IF
555    END DO
556
557    !======================================================================
558    !   Parameterization of ozone chemistry
559    !======================================================================
560
561    IF (id_o3 /= 0) THEN
562      lmt_pas = NINT(86400. / pdtphys)
563      IF (MOD(nstep - 1, lmt_pas) == 0) THEN
564        ! Once per day, update the coefficients for ozone chemistry:
565        CALL regr_pr_comb_coefoz(julien, xlat, paprs, pplay)
566      END IF
567      CALL o3_chem(julien, gmtime, t_seri, zmasse, pdtphys, xlat, &
568              xlon, tr_seri(:, :, id_o3))
569    END IF
570
571  END SUBROUTINE traclmdz
572
573
574  SUBROUTINE traclmdz_to_restart(trs_out)
575    ! This SUBROUTINE is called from phyredem.F where the module
576    ! variable trs is written to restart file (restartphy.nc)
577    USE dimphy
578    USE infotrac_phy, ONLY: nbtr
579
580    REAL, DIMENSION(klon, nbtr), INTENT(OUT) :: trs_out
581    INTEGER :: ierr
582
583    IF (ALLOCATED(trs)) THEN
584      trs_out(:, :) = trs(:, :)
585    ELSE
586      ! No previous allocate of trs. This is the case for create_etat0_limit.
587      trs_out(:, :) = 0.0
588    END IF
589
590  END SUBROUTINE traclmdz_to_restart
591
592
593END MODULE traclmdz_mod
Note: See TracBrowser for help on using the repository browser.