My Project
 All Classes Files Functions Variables Macros
traclmdz_mod.F90
Go to the documentation of this file.
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  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 
98  include "indicesol.h"
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, zu10m, zv10m, &
340 !! tr_seri, source, solsym, d_tr_cl, zmasse) !RomP
341  tr_seri, source, solsym, d_tr_cl,d_tr_dec, zmasse) !RomP
342 
343  USE dimphy
344  USE infotrac
345  USE regr_pr_comb_coefoz_m, ONLY: regr_pr_comb_coefoz
346  USE o3_chem_m, ONLY: o3_chem
347  USE carbon_cycle_mod, ONLY : carbon_cycle, carbon_cycle_tr, carbon_cycle_cpl
348  include "YOMCST.h"
349  include "indicesol.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) :: zu10m ! vent zonal 10m (m/s)
387  REAL,DIMENSION(klon),INTENT(IN) :: zv10m ! vent zonal 10m (m/s)
388 
389 ! Arguments necessaires pour les sources et puits de traceur:
390  REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: ftsol ! Temperature du sol (surf)(Kelvin)
391  REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: pctsrf ! Pourcentage de sol f(nature du sol)
392 
393 ! InOutput argument
394  REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT) :: tr_seri ! Concentration Traceur [U/KgA]
395 
396 ! Output argument
397  CHARACTER(len=8),DIMENSION(nbtr), INTENT(OUT) :: solsym
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 
641 END MODULE traclmdz_mod