source: LMDZ4/branches/LMDZ4_AR5/libf/phylmd/traclmdz_mod.F90 @ 3400

Last change on this file since 3400 was 1410, checked in by jghattas, 14 years ago

Bug fix : Now tracer Aga can be used in the same run as RN and
PB. Still, in this branche, RN and PB must be tracer number 3 and 4 in tracer.def.

File size: 13.3 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
9  REAL,DIMENSION(:,:),ALLOCATABLE,SAVE :: masktr   ! Masque reservoir de sol traceur
10!$OMP THREADPRIVATE(masktr)                        ! Masque de l'echange avec la surface (1 = reservoir) ou (possible >= 1 )
11  REAL,DIMENSION(:,:),ALLOCATABLE,SAVE :: fshtr    ! Flux surfacique dans le reservoir de sol
12!$OMP THREADPRIVATE(fshtr)
13  REAL,DIMENSION(:),ALLOCATABLE,SAVE   :: hsoltr   ! Epaisseur equivalente du reservoir de sol
14!$OMP THREADPRIVATE(hsoltr)
15!
16!Radioelements:
17!--------------
18!
19  REAL,DIMENSION(:),ALLOCATABLE,SAVE   :: tautr    ! Constante de decroissance radioactive
20!$OMP THREADPRIVATE(tautr)
21  REAL,DIMENSION(:),ALLOCATABLE,SAVE   :: vdeptr   ! Vitesse de depot sec dans la couche Brownienne
22!$OMP THREADPRIVATE(vdeptr)
23  REAL,DIMENSION(:),ALLOCATABLE,SAVE   :: scavtr   ! Coefficient de lessivage
24!$OMP THREADPRIVATE(scavtr)
25  REAL,DIMENSION(:,:),ALLOCATABLE,SAVE :: srcbe    ! Production du beryllium7 dans l atmosphere (U/s/kgA)
26!$OMP THREADPRIVATE(srcbe)
27
28  LOGICAL,DIMENSION(:),ALLOCATABLE,SAVE :: radio    ! radio(it)   = true  => decroisssance radioactive
29!$OMP THREADPRIVATE(radio) 
30
31  REAL,DIMENSION(:,:),ALLOCATABLE,SAVE :: trs     ! Conc. radon ds le sol
32!$OMP THREADPRIVATE(trs)
33
34  INTEGER,SAVE :: id_be       ! Activation et position du traceur Be7 [ id_be=0 -> desactive ]
35!$OMP THREADPRIVATE(id_be)
36
37  INTEGER,SAVE :: id_aga      ! Identification number for tracer : Age of stratospheric air
38!$OMP THREADPRIVATE(id_aga)
39  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.
40!$OMP THREADPRIVATE(lev_1p5km)
41
42  INTEGER,SAVE :: id_rn, id_pb ! Identification number for tracer : radon (Rn222), lead (Pb210)
43  !$OMP THREADPRIVATE(id_rn, id_pb)
44  LOGICAL,SAVE :: rnpb=.FALSE. ! Presence du couple Rn222, Pb210
45!$OMP THREADPRIVATE(rnpb)
46
47 
48
49CONTAINS
50
51  SUBROUTINE traclmdz_from_restart(trs_in)
52    ! This subroutine initialize the module saved variable trs with values from restart file (startphy.nc).
53    ! This subroutine is called from phyetat0 after the field trs_in has been read.
54   
55    USE dimphy
56    USE infotrac
57    IMPLICIT NONE
58   
59    ! Input argument
60    REAL,DIMENSION(klon,nbtr), INTENT(IN) :: trs_in
61   
62    ! Local variables
63    INTEGER :: ierr
64   
65    ! Allocate restart variables trs
66    ALLOCATE( trs(klon,nbtr), stat=ierr)
67    IF (ierr /= 0) CALL abort_gcm('traclmdz_from_restart', 'pb in allocation 1',1)
68   
69    ! Initialize trs with values read from restart file
70    trs(:,:) = trs_in(:,:)
71   
72  END SUBROUTINE traclmdz_from_restart
73
74
75  SUBROUTINE traclmdz_init(pctsrf, ftsol, tr_seri, aerosol, lessivage)
76    ! This subroutine allocates and initialize module variables and control variables.
77    USE dimphy
78    USE infotrac
79    USE carbon_cycle_mod, ONLY : carbon_cycle_init, carbon_cycle_tr, carbon_cycle_cpl
80
81    IMPLICIT NONE
82
83    INCLUDE "indicesol.h"
84
85! Input variables
86    REAL,DIMENSION(klon,nbsrf),INTENT(IN)     :: pctsrf ! Pourcentage de sol f(nature du sol)
87    REAL,DIMENSION(klon,nbsrf),INTENT(IN)     :: ftsol  ! Temperature du sol (surf)(Kelvin)
88    REAL,DIMENSION(klon,klev,nbtr),INTENT(IN) :: tr_seri! Concentration Traceur [U/KgA] 
89
90! Output variables
91    LOGICAL,DIMENSION(nbtr), INTENT(OUT) :: aerosol
92    LOGICAL,INTENT(OUT)                  :: lessivage
93       
94! Local variables   
95    INTEGER :: ierr, it, iiq
96   
97! --------------------------------------------
98! Allocation
99! --------------------------------------------
100
101    ALLOCATE( scavtr(nbtr), stat=ierr)
102    IF (ierr /= 0) CALL abort_gcm('traclmdz_init', 'pb in allocation 9',1)
103    scavtr(:)=1.
104   
105    ALLOCATE( radio(nbtr), stat=ierr)
106    IF (ierr /= 0) CALL abort_gcm('traclmdz_init', 'pb in allocation 11',1)
107    radio(:) = .false.    ! Par defaut pas decroissance radioactive
108   
109    ALLOCATE( masktr(klon,nbtr), stat=ierr)
110    IF (ierr /= 0) CALL abort_gcm('traclmdz_init', 'pb in allocation 2',1)
111   
112    ALLOCATE( fshtr(klon,nbtr), stat=ierr)
113    IF (ierr /= 0) CALL abort_gcm('traclmdz_init', 'pb in allocation 3',1)
114   
115    ALLOCATE( hsoltr(nbtr), stat=ierr)
116    IF (ierr /= 0) CALL abort_gcm('traclmdz_init', 'pb in allocation 4',1)
117   
118    ALLOCATE( tautr(nbtr), stat=ierr)
119    IF (ierr /= 0) CALL abort_gcm('traclmdz_init', 'pb in allocation 5',1)
120    tautr(:)  = 0.
121   
122    ALLOCATE( vdeptr(nbtr), stat=ierr)
123    IF (ierr /= 0) CALL abort_gcm('traclmdz_init', 'pb in allocation 6',1)
124    vdeptr(:) = 0.
125
126
127    lessivage  = .TRUE.
128    aerosol(:) = .FALSE.  ! Tous les traceurs sont des gaz par defaut
129   
130!
131! Recherche des traceurs connus : Be7, CO2,...
132! --------------------------------------------
133    id_rn=0; id_pb=0; id_aga=0; id_be=0
134    DO it=1,nbtr
135       iiq=niadv(it+2)
136       IF ( tname(iiq) == "RN" ) THEN
137          id_rn=it ! radon
138       ELSE IF ( tname(iiq) == "PB") THEN
139          id_pb=it ! plomb
140       ELSE IF ( tname(iiq) == "Aga" .OR. tname(iiq) == "AGA" ) THEN
141          ! Age of stratospheric air
142          id_aga=it
143          radio(id_aga) = .FALSE.
144          aerosol(id_aga) = .FALSE.
145          pbl_flg(id_aga) = 0
146
147          ! Find the first model layer above 1.5km from the surface
148          IF (klev>=30) THEN
149             lev_1p5km=6   ! NB! This value is for klev=39
150          ELSE IF (klev>=10) THEN
151             lev_1p5km=5   ! NB! This value is for klev=19
152          ELSE
153             lev_1p5km=klev/2
154          END IF
155       ELSE IF ( tname(iiq) == "BE" .OR. tname(iiq) == "Be" .OR.  &
156            tname(iiq) == "BE7" .OR. tname(iiq) == "Be7" ) THEN 
157          ! Recherche du Beryllium 7
158          id_be=it
159          ALLOCATE( srcbe(klon,klev) )
160          radio(id_be) = .TRUE.
161          aerosol(id_be) = .TRUE. ! le Be est un aerosol
162          CALL init_be(pctsrf,masktr(:,id_be),tautr(id_be),vdeptr(id_be),scavtr(id_be),srcbe)
163          WRITE(*,*) 'Initialisation srcBe: OK'
164       END IF
165       
166    END DO
167!
168! Valeurs specifiques pour les traceurs Rn222 et Pb210
169! ----------------------------------------------
170    IF ( id_rn/=0 .AND. id_pb/=0 ) THEN
171       rnpb = .TRUE.
172       radio(1)= .TRUE.
173       radio(2)= .TRUE.
174       pbl_flg(1) = 0 ! au lieu de clsol=true ! CL au sol calcule
175       pbl_flg(2) = 0 ! au lieu de clsol=true
176       
177       aerosol(2) = .TRUE. ! le Pb est un aerosol
178       
179       CALL initrrnpb (ftsol,pctsrf,masktr,fshtr,hsoltr,tautr,vdeptr,scavtr)
180    END IF
181
182!
183! Initialisation de module carbon_cycle_mod
184! ----------------------------------------------
185    IF (carbon_cycle_tr .OR. carbon_cycle_cpl) THEN
186       CALL carbon_cycle_init(tr_seri, aerosol, radio)
187    END IF
188
189  END SUBROUTINE traclmdz_init
190
191  SUBROUTINE traclmdz(                           &
192       nstep,    pdtphys,      t_seri,           &
193       paprs,    pplay,        cdragh,  coefh,   &
194       yu1,      yv1,          ftsol,   pctsrf,  &
195       xlat,     couchelimite,                   &
196       tr_seri,  source,       solsym,  d_tr_cl)
197   
198    USE dimphy
199    USE infotrac
200    USE carbon_cycle_mod, ONLY : carbon_cycle, carbon_cycle_tr, carbon_cycle_cpl
201   
202    IMPLICIT NONE
203   
204    INCLUDE "YOMCST.h"
205    INCLUDE "indicesol.h"
206
207!==========================================================================
208!                   -- DESCRIPTION DES ARGUMENTS --
209!==========================================================================
210
211! Input arguments
212!
213!Configuration grille,temps:
214    INTEGER,INTENT(IN) :: nstep      ! nombre d'appels de la physiq
215    REAL,INTENT(IN)    :: pdtphys    ! Pas d'integration pour la physique (seconde) 
216    REAL,DIMENSION(klon),INTENT(IN) :: xlat    ! latitudes pour chaque point
217
218!
219!Physique:
220!--------
221    REAL,DIMENSION(klon,klev),INTENT(IN)   :: t_seri  ! Temperature
222    REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs   ! pression pour chaque inter-couche (en Pa)
223    REAL,DIMENSION(klon,klev),INTENT(IN)   :: pplay   ! pression pour le mileu de chaque couche (en Pa)
224
225
226!Couche limite:
227!--------------
228!
229    REAL,DIMENSION(klon,klev),INTENT(IN) :: cdragh     ! coeff drag pour T et Q
230    REAL,DIMENSION(klon,klev),INTENT(IN) :: coefh      ! coeff melange CL (m**2/s)
231    REAL,DIMENSION(klon),INTENT(IN)      :: yu1        ! vents au premier niveau
232    REAL,DIMENSION(klon),INTENT(IN)      :: yv1        ! vents au premier niveau
233    LOGICAL,INTENT(IN)                   :: couchelimite
234
235! Arguments necessaires pour les sources et puits de traceur:
236    REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: ftsol  ! Temperature du sol (surf)(Kelvin)
237    REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: pctsrf ! Pourcentage de sol f(nature du sol)
238
239
240! InOutput argument
241    REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT)  :: tr_seri ! Concentration Traceur [U/KgA] 
242
243! Output argument
244    CHARACTER(len=8),DIMENSION(nbtr), INTENT(OUT) :: solsym
245    REAL,DIMENSION(klon,nbtr), INTENT(OUT)        :: source  ! a voir lorsque le flux de surface est prescrit
246    REAL,DIMENSION(klon,klev,nbtr), INTENT(OUT)   :: d_tr_cl ! Td couche limite/traceur
247
248!=======================================================================================
249!                        -- VARIABLES LOCALES TRACEURS --
250!=======================================================================================
251
252    INTEGER :: i, k, it
253
254    REAL,DIMENSION(klon)           :: d_trs    ! Td dans le reservoir
255    REAL,DIMENSION(klon,klev)      :: delp     ! epaisseur de couche (Pa)
256   
257    REAL,DIMENSION(klon,klev,nbtr) :: d_tr_dec ! Td radioactive
258    REAL                           :: zrho      ! Masse Volumique de l'air KgA/m3
259
260!
261!
262!=================================================================
263!  Ajout de la production en  Be7 (Beryllium) srcbe U/s/kgA
264!=================================================================
265!
266    IF ( id_be /= 0 ) THEN
267       DO k = 1, klev
268          DO i = 1, klon
269             tr_seri(i,k,id_be) = tr_seri(i,k,id_be)+srcbe(i,k)*pdtphys
270          END DO
271       END DO
272       WRITE(*,*) 'Ajout srcBe dans tr_seri: OK'
273    END IF
274 
275
276    DO it=1,nbtr
277       WRITE(solsym(it),'(i2)') it
278    END DO
279
280!
281! Update tracer : Age of stratospheric air
282!
283    IF (id_aga/=0) THEN
284       
285       ! Bottom layers
286       DO k = 1, lev_1p5km
287          tr_seri(:,k,id_aga) = 0.0
288       END DO
289       
290       ! Layers above 1.5km
291       DO k = lev_1p5km+1,klev-1
292          tr_seri(:,k,id_aga) = tr_seri(:,k,id_aga) + pdtphys
293       END DO
294       
295       ! Top layer
296       tr_seri(:,klev,id_aga) = tr_seri(:,klev-1,id_aga)
297       
298    END IF
299
300!======================================================================
301!     -- Calcul de l'effet de la couche limite --
302!======================================================================
303
304    IF (couchelimite) THEN             
305       source(:,:) = 0.0
306
307       IF (id_be /=0) THEN
308          DO i=1, klon
309             zrho = pplay(i,1)/t_seri(i,1)/RD
310             source(i,id_be) = - vdeptr(id_be)*tr_seri(i,1,id_be)*zrho
311          END DO
312       END IF
313
314    END IF
315   
316
317    DO k = 1, klev
318       DO i = 1, klon
319          delp(i,k) = paprs(i,k)-paprs(i,k+1)
320       END DO
321    END DO
322   
323    DO it=1, nbtr
324       IF (couchelimite .AND. pbl_flg(it) == 0 .AND. (it==id_rn .OR. it==id_pb)) THEN ! couche limite avec quantite dans le sol calculee
325          CALL cltracrn(it, pdtphys, yu1, yv1,     &
326               cdragh, coefh,t_seri,ftsol,pctsrf,  &
327               tr_seri(:,:,it),trs(:,it),          &
328               paprs, pplay, delp,                 &
329               masktr(:,it),fshtr(:,it),hsoltr(it),&
330               tautr(it),vdeptr(it),               &
331               xlat,d_tr_cl(:,:,it),d_trs)
332         
333          DO k = 1, klev
334             DO i = 1, klon
335                tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cl(i,k,it)
336             END DO
337          END DO
338       
339          ! Traceur dans le reservoir sol
340          DO i = 1, klon
341             trs(i,it) = trs(i,it) + d_trs(i)
342          END DO
343       END IF
344    END DO
345           
346!======================================================================
347!   Calcul de l'effet du puits radioactif
348!======================================================================
349    CALL radio_decay (radio,rnpb,pdtphys,tautr,tr_seri,d_tr_dec)
350 
351    DO it=1,nbtr
352       IF(radio(it)) then     
353          DO k = 1, klev
354             DO i = 1, klon
355                tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_dec(i,k,it)
356             END DO
357          END DO
358          CALL minmaxqfi(tr_seri(:,:,it),0.,1.e33,'puits rn it='//solsym(it))
359       END IF
360    END DO
361
362!======================================================================
363!   Calcul de cycle de carbon
364!======================================================================
365    IF (carbon_cycle_tr .OR. carbon_cycle_cpl) THEN
366       CALL carbon_cycle(nstep, pdtphys, pctsrf, tr_seri)
367    END IF
368
369  END SUBROUTINE traclmdz
370
371
372  SUBROUTINE traclmdz_to_restart(trs_out)
373    ! This subroutine is called from phyredem.F where the module
374    ! variable trs is written to restart file (restartphy.nc)
375    USE dimphy
376    USE infotrac
377   
378    IMPLICIT NONE
379   
380    REAL,DIMENSION(klon,nbtr), INTENT(OUT) :: trs_out
381    INTEGER :: ierr
382
383    IF ( ALLOCATED(trs) ) THEN
384       trs_out(:,:) = trs(:,:)
385    ELSE
386       ! No previous allocate of trs. This is the case for create_etat0_limit.
387       trs_out(:,:) = 0.0
388    END IF
389   
390  END SUBROUTINE traclmdz_to_restart
391 
392
393END MODULE traclmdz_mod
Note: See TracBrowser for help on using the repository browser.