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

Last change on this file since 5103 was 5103, checked in by abarral, 8 weeks ago

Handle CPP_INLANDSIS in lmdz_cppkeys_wrapper.F90
Remove obsolete key wrgrads_thermcell, _ADV_HALO, _ADV_HALLO, isminmax
Remove redundant uses of CPPKEY_INCA (thanks acozic)
Remove obsolete misc/write_field.F90
Remove unused ioipsl_* wrappers
Remove calls to WriteField_u with wrong signature
Convert .F -> .[fF]90
(lint) uppercase fortran operators
[note: 1d and iso still broken - working on it]

  • 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
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 mod_grid_phy_lmdz
94    USE mod_phys_lmdz_para
95    USE indice_sol_mod
96    USE print_control_mod, 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
158    lessivage  = .TRUE.
159!!jyg(20130206) : le choix d activation du lessivage est fait dans phytrac avec iflag_lscav
160!!    CALL getin('lessivage',lessivage)
161!!    if(lessivage) then
162!!     PRINT*,'lessivage lsc ON'
163!!    else
164!!     PRINT*,'lessivage lsc OFF'
165!!    endif
166    aerosol(:) = .FALSE.  ! Tous les traceurs sont des gaz par defaut
167
168! Recherche des traceurs connus : Be7, O3, CO2,...
169! --------------------------------------------
170    id_rn=0; id_pb=0; id_aga=0; id_be=0; id_o3=0
171    id_pcsat=0; id_pcocsat=0; id_pcq=0; id_pcs0=0; id_pcos0=0; id_pcq0=0
172    it = 0
173    DO iq = 1, nqtot
174       IF(.NOT.(tracers(iq)%isInPhysics)) CYCLE
175       it = it+1
176       SELECT CASE(tracers(iq)%name)
177         CASE("RN");                  id_rn     = it ! radon
178         CASE("PB");                  id_pb     = it ! plomb
179         CASE("Aga","AGA");           id_aga    = it ! Age of stratospheric air
180         CASE("Be","BE","Be7","BE7"); id_be     = it ! Recherche du Beryllium 7
181         CASE("o3","O3");             id_o3     = it ! Recherche de l'ozone
182         CASE("pcsat",  "Pcsat");     id_pcsat  = it
183         CASE("pcocsat","Pcocsat");   id_pcocsat= it
184         CASE("pcq",    "Pcq");       id_pcq    = it
185         CASE("pcs0",   "Pcs0");      id_pcs0   = it
186         CASE("pcos0",  "Pcos0");     id_pcos0  = it
187         CASE("pcq0",   "Pcq0");      id_pcq0   = it
188         CASE DEFAULT
189           WRITE(lunout,*) 'This is an unknown tracer in LMDZ : ', trim(tracers(iq)%name)
190       END SELECT
191
192       SELECT CASE(tracers(iq)%name)
193         CASE("PB")                        !--- RomP >>> profil initial de PB210
194           OPEN(ilesfil2,file='prof.pb210',status='old',iostat=irr2)
195           IF(irr2 == 0) THEN
196             READ(ilesfil2,*) kradio2
197             WRITE(lunout,*)'number of levels for pb210 profile ',kradio2
198             DO k=kradio2,1,-1; READ (ilesfil2,*) plomb(:,k); 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); END DO
228             CLOSE(ilesfil)
229             tr_seri(:,:,id_be)=beryllium(:,:)
230           ELSE
231             WRITE(lunout,*)'Prof. Be7 does not exist: use restart values'
232           END IF
233         CASE("o3","O3")                    !--- Parametrisation par la chimie de Cariolle
234           CALL alloc_coefoz                !--- Allocate ozone coefficients
235           CALL press_coefoz                !--- Read input pressure levels
236         CASE("pcs0","Pcs0", "pcos0","Pcos0", "pcq0","Pcq0")
237           conv_flg(it)=0                   !--- No transport by convection for this tracer
238       END SELECT
239    END DO
240
241! Valeurs specifiques pour les traceurs Rn222 et Pb210
242! ----------------------------------------------
243    IF ( id_rn/=0 .AND. id_pb/=0 ) THEN
244       rnpb = .TRUE.
245       radio(id_rn)= .TRUE.
246       radio(id_pb)= .TRUE.
247       pbl_flg(id_rn) = 0 ! au lieu de clsol=true ! CL au sol calcule
248       pbl_flg(id_pb) = 0 ! au lieu de clsol=true
249       aerosol(id_rn) = .FALSE.
250       aerosol(id_pb) = .TRUE. ! le Pb est un aerosol
251       
252       CALL initrrnpb (ftsol,pctsrf,masktr,fshtr,hsoltr,tautr,vdeptr,scavtr)
253    END IF
254
255! Check if all tracers have restart values
256! ----------------------------------------------
257    it = 0
258    DO iq = 1, nqtot
259       IF(.NOT.(tracers(iq)%isAdvected .AND. tracers(iq)%isInPhysics)) CYCLE
260       it = it+1
261       ! Test if tracer is zero everywhere.
262       ! Done by master process MPI and master thread OpenMP
263       CALL gather(tr_seri(:,:,it),varglo)
264       IF (is_mpi_root .AND. is_omp_root) THEN
265          mintmp=MINVAL(varglo,dim=1)
266          maxtmp=MAXVAL(varglo,dim=1)
267          IF (MINVAL(mintmp,dim=1)==0. .AND. MAXVAL(maxtmp,dim=1)==0.) THEN
268             ! Tracer is zero everywhere
269             zero=.TRUE.
270          ELSE
271             zero=.FALSE.
272          END IF
273       END IF
274
275       ! Distribute variable at all processes
276       CALL bcast(zero)
277
278       ! Initalize tracer that was not found in restart file.
279       IF (zero) THEN
280          ! The tracer was not found in restart file or it was equal zero everywhere.
281          WRITE(lunout,*) "The tracer ",trim(tracers(iq)%name)," will be initialized"
282          IF (it==id_pcsat .OR. it==id_pcq .OR. &
283               it==id_pcs0 .OR. it==id_pcq0) THEN
284             tr_seri(:,:,it) = 100.
285          ELSE IF (it==id_pcocsat .OR. it==id_pcos0) THEN
286             DO i = 1, klon
287                IF ( pctsrf (i, is_oce) == 0. ) THEN
288                   tr_seri(i,:,it) = 0.
289                ELSE
290                   tr_seri(i,:,it) = 100.
291                END IF
292             END DO
293          ELSE
294             ! No specific initialization exist for this tracer
295             tr_seri(:,:,it) = 0.
296          END IF
297       END IF
298    END DO
299
300  END SUBROUTINE traclmdz_init
301
302  SUBROUTINE traclmdz(nstep, julien, gmtime, pdtphys, t_seri, paprs, pplay, &
303       cdragh,  coefh, yu1, yv1, ftsol, pctsrf, xlat, xlon, couchelimite, sh, &
304       rh, pphi, ustar, wstar, ale_bl, ale_wake,  zu10m, zv10m, &
305       tr_seri, source, d_tr_cl,d_tr_dec, zmasse)               !RomP
306   
307    USE dimphy
308    USE infotrac_phy, ONLY: nbtr, pbl_flg
309    USE strings_mod,  ONLY: int2str
310    USE regr_pr_comb_coefoz_m, ONLY: regr_pr_comb_coefoz
311    USE o3_chem_m, ONLY: o3_chem
312    USE indice_sol_mod
313
314    INCLUDE "YOMCST.h"
315
316!==========================================================================
317!                   -- DESCRIPTION DES ARGUMENTS --
318!==========================================================================
319
320! Input arguments
321
322!Configuration grille,temps:
323    INTEGER,INTENT(IN) :: nstep      ! nombre d'appels de la physiq
324    INTEGER,INTENT(IN) :: julien     ! Jour julien
325    REAL,INTENT(IN)    :: gmtime
326    REAL,INTENT(IN)    :: pdtphys    ! Pas d'integration pour la physique (seconde) 
327    REAL,DIMENSION(klon),INTENT(IN) :: xlat    ! latitudes pour chaque point
328    REAL, INTENT(IN):: xlon(:) ! dim(klon) longitude
329
330!Physique:
331!--------
332    REAL,DIMENSION(klon,klev),INTENT(IN)   :: t_seri  ! Temperature
333    REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs   ! pression pour chaque inter-couche (en Pa)
334    REAL,DIMENSION(klon,klev),INTENT(IN)   :: pplay   ! pression pour le mileu de chaque couche (en Pa)
335    REAL,intent(in):: zmasse (:, :)   ! dim(klon,klev) density of air, in kg/m2
336
337
338!Couche limite:
339!--------------
340
341    REAL,DIMENSION(klon),INTENT(IN)      :: cdragh     ! coeff drag pour T et Q
342    REAL,DIMENSION(klon,klev),INTENT(IN) :: coefh      ! coeff melange CL (m**2/s)
343    REAL,DIMENSION(klon),INTENT(IN)      :: yu1        ! vents au premier niveau
344    REAL,DIMENSION(klon),INTENT(IN)      :: yv1        ! vents au premier niveau
345    LOGICAL,INTENT(IN)                   :: couchelimite
346    REAL,DIMENSION(klon,klev),INTENT(IN) :: sh         ! humidite specifique
347    REAL,DIMENSION(klon,klev),INTENT(IN) :: rh      ! Humidite relative
348    REAL,DIMENSION(klon,klev),INTENT(IN) :: pphi    ! geopotentie
349    REAL,DIMENSION(klon),INTENT(IN)      :: ustar   ! ustar (m/s)
350    REAL,DIMENSION(klon),INTENT(IN)      :: wstar,ale_bl,ale_wake   ! wstar (m/s) and Avail. Lifti. Energ.
351    REAL,DIMENSION(klon),INTENT(IN)      :: zu10m   ! vent zonal 10m (m/s)
352    REAL,DIMENSION(klon),INTENT(IN)      :: zv10m   ! vent zonal 10m (m/s)
353
354! Arguments necessaires pour les sources et puits de traceur:
355    REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: ftsol  ! Temperature du sol (surf)(Kelvin)
356    REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: pctsrf ! Pourcentage de sol f(nature du sol)
357
358! InOutput argument
359    REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT)  :: tr_seri ! Concentration Traceur [U/KgA] 
360
361! Output argument
362    REAL,DIMENSION(klon,nbtr), INTENT(OUT)        :: source  ! a voir lorsque le flux de surface est prescrit
363    REAL,DIMENSION(klon,klev,nbtr), INTENT(OUT)   :: d_tr_cl ! Td couche limite/traceur
364
365!=======================================================================================
366!                        -- VARIABLES LOCALES TRACEURS --
367!=======================================================================================
368
369    INTEGER :: i, k, it
370    INTEGER :: lmt_pas ! number of time steps of "physics" per day
371
372    REAL,DIMENSION(klon)           :: d_trs    ! Td dans le reservoir
373    REAL,DIMENSION(klon,klev)      :: qsat     ! pression de la vapeur a saturation
374    REAL,DIMENSION(klon,klev,nbtr) :: d_tr_dec ! Td radioactive
375    REAL                           :: zrho     ! Masse Volumique de l'air KgA/m3
376    REAL                           :: amn, amx
377
378!=================================================================
379!  Ajout de la production en  Be7 (Beryllium) srcbe U/s/kgA
380!=================================================================
381
382    IF ( id_be /= 0 ) THEN
383       DO k = 1, klev
384          DO i = 1, klon
385             tr_seri(i,k,id_be) = tr_seri(i,k,id_be)+srcbe(i,k)*pdtphys
386          END DO
387       END DO
388       WRITE(*,*) 'Ajout srcBe dans tr_seri: OK'
389    END IF
390 
391
392!=================================================================
393! Update pseudo-vapor tracers
394!=================================================================
395
396    CALL q_sat(klon*klev,t_seri,pplay,qsat)
397
398    IF ( id_pcsat /= 0 ) THEN
399     DO k = 1, klev
400      DO i = 1, klon
401         IF ( pplay(i,k)>=85000.) THEN
402            tr_seri(i,k,id_pcsat) = qsat(i,k)
403         ELSE
404            tr_seri(i,k,id_pcsat) = MIN (qsat(i,k), tr_seri(i,k,id_pcsat))           
405         END IF
406      END DO
407     END DO
408    END IF
409
410    IF ( id_pcocsat /= 0 ) THEN
411     DO k = 1, klev
412      DO i = 1, klon
413         IF ( pplay(i,k)>=85000.) THEN
414            IF ( pctsrf (i, is_oce) > 0. ) THEN
415               tr_seri(i,k,id_pcocsat) = qsat(i,k)
416            ELSE
417               tr_seri(i,k,id_pcocsat) = 0.
418          END IF
419       ELSE
420          tr_seri(i,k,id_pcocsat) = MIN (qsat(i,k), tr_seri(i,k,id_pcocsat))
421       END IF
422      END DO
423     END DO
424    END IF
425
426    IF ( id_pcq /= 0 ) THEN
427     DO k = 1, klev
428      DO i = 1, klon
429         IF ( pplay(i,k)>=85000.) THEN
430            tr_seri(i,k,id_pcq) = sh(i,k)
431         ELSE
432            tr_seri(i,k,id_pcq) = MIN (qsat(i,k), tr_seri(i,k,id_pcq))   
433         END IF
434      END DO
435     END DO
436    END IF
437
438
439    IF ( id_pcs0 /= 0 ) THEN
440     DO k = 1, klev
441      DO i = 1, klon
442         IF ( pplay(i,k)>=85000.) THEN
443            tr_seri(i,k,id_pcs0) = qsat(i,k)
444         ELSE
445            tr_seri(i,k,id_pcs0) = MIN (qsat(i,k), tr_seri(i,k,id_pcs0))   
446         END IF
447      END DO
448     END DO
449    END IF
450
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
469    IF ( id_pcq0 /= 0 ) THEN
470     DO k = 1, klev
471      DO i = 1, klon
472         IF ( pplay(i,k)>=85000.) THEN
473            tr_seri(i,k,id_pcq0) = sh(i,k)
474         ELSE
475            tr_seri(i,k,id_pcq0) = MIN (qsat(i,k), tr_seri(i,k,id_pcq0))
476         END IF
477      END DO
478     END DO
479    END IF
480
481!=================================================================
482! Update tracer : Age of stratospheric air
483!=================================================================
484    IF (id_aga/=0) THEN
485       
486       ! Bottom layers
487       DO k = 1, lev_1p5km
488          tr_seri(:,k,id_aga) = 0.0
489       END DO
490       
491       ! Layers above 1.5km
492       DO k = lev_1p5km+1,klev-1
493          tr_seri(:,k,id_aga) = tr_seri(:,k,id_aga) + pdtphys
494       END DO
495       
496       ! Top layer
497       tr_seri(:,klev,id_aga) = tr_seri(:,klev-1,id_aga)
498       
499    END IF
500
501!======================================================================
502!     -- Calcul de l'effet de la couche limite --
503!======================================================================
504
505    IF (couchelimite) THEN             
506       source(:,:) = 0.0
507
508       IF (id_be /=0) THEN
509          DO i=1, klon
510             zrho = pplay(i,1)/t_seri(i,1)/RD
511             source(i,id_be) = - vdeptr(id_be)*tr_seri(i,1,id_be)*zrho
512          END DO
513       END IF
514
515    END IF
516   
517    DO it=1, nbtr
518       IF (couchelimite .AND. pbl_flg(it) == 0 .AND. (it==id_rn .OR. it==id_pb)) THEN
519          ! couche limite avec quantite dans le sol calculee
520          CALL cltracrn(it, pdtphys, yu1, yv1,     &
521               cdragh, coefh,t_seri,ftsol,pctsrf,  &
522               tr_seri(:,:,it),trs(:,it),          &
523               paprs, pplay, zmasse * rg, &
524               masktr(:,it),fshtr(:,it),hsoltr(it),&
525               tautr(it),vdeptr(it),               &
526               xlat,d_tr_cl(:,:,it),d_trs)
527         
528          DO k = 1, klev
529             DO i = 1, klon
530                tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cl(i,k,it)
531             END DO
532          END DO
533       
534          ! Traceur dans le reservoir sol
535          DO i = 1, klon
536             trs(i,it) = trs(i,it) + d_trs(i)
537          END DO
538       END IF
539    END DO
540         
541
542!======================================================================
543!   Calcul de l'effet du puits radioactif
544!======================================================================
545    CALL radio_decay (radio,rnpb,pdtphys,tautr,tr_seri,d_tr_dec)
546
547    DO it=1,nbtr
548       IF(radio(it)) then     
549          DO k = 1, klev
550             DO i = 1, klon
551                tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_dec(i,k,it)
552             END DO
553          END DO
554          CALL minmaxqfi(tr_seri(:,:,it),0.,1.e33,'puits rn it='//TRIM(int2str(it)))
555       END IF
556    END DO
557
558!======================================================================
559!   Parameterization of ozone chemistry
560!======================================================================
561
562    IF (id_o3 /= 0) then
563       lmt_pas = NINT(86400./pdtphys)
564       IF (MOD(nstep - 1, lmt_pas) == 0) THEN
565          ! Once per day, update the coefficients for ozone chemistry:
566          CALL regr_pr_comb_coefoz(julien, xlat, paprs, pplay)
567       END IF
568       CALL o3_chem(julien, gmtime, t_seri, zmasse, pdtphys, xlat, &
569            xlon, tr_seri(:, :, id_o3))
570    END IF
571
572  END SUBROUTINE traclmdz
573
574
575  SUBROUTINE traclmdz_to_restart(trs_out)
576    ! This SUBROUTINE is called from phyredem.F where the module
577    ! variable trs is written to restart file (restartphy.nc)
578    USE dimphy
579    USE infotrac_phy, ONLY: nbtr
580   
581    REAL,DIMENSION(klon,nbtr), INTENT(OUT) :: trs_out
582    INTEGER :: ierr
583
584    IF ( ALLOCATED(trs) ) THEN
585       trs_out(:,:) = trs(:,:)
586    ELSE
587       ! No previous allocate of trs. This is the case for create_etat0_limit.
588       trs_out(:,:) = 0.0
589    END IF
590   
591  END SUBROUTINE traclmdz_to_restart
592 
593
594END MODULE traclmdz_mod
Note: See TracBrowser for help on using the repository browser.