source: LMDZ5/trunk/libf/phylmd/traclmdz_mod.F90 @ 1651

Last change on this file since 1651 was 1579, checked in by jghattas, 13 years ago

Added latitudes and longitudes in traclmdz_init to simplify initialization of new tracers.

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