source: LMDZ5/branches/testing/libf/phylmd/traclmdz_mod.F90 @ 2254

Last change on this file since 2254 was 2187, checked in by Laurent Fairhead, 10 years ago

Merged trunk changes -r2158:2186 into testing branch.

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