source: LMDZ4/trunk/libf/phylmd/traclmdz_mod.F90 @ 5427

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