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

Last change on this file since 1320 was 1279, checked in by Laurent Fairhead, 15 years ago

Merged LMDZ4-dev branch changes r1241:1278 into the trunk
Running trunk and LMDZ4-dev in LMDZOR configuration on local
machine (sequential) and SX8 (4-proc) yields identical results
(restart and restartphy are identical binarily)
Log history from r1241 to r1278 is available by switching to
source:LMDZ4/branches/LMDZ4-dev-20091210
File size: 11.6 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  LOGICAL,SAVE :: rnpb=.TRUE. ! Presence du couple Rn222, Pb210
38!$OMP THREADPRIVATE(rnpb)
39
40
41CONTAINS
42
43  SUBROUTINE traclmdz_from_restart(trs_in)
44    ! This subroutine initialize the module saved variable trs with values from restart file (startphy.nc).
45    ! This subroutine is called from phyetat0 after the field trs_in has been read.
46   
47    USE dimphy
48    USE infotrac
49    IMPLICIT NONE
50   
51    ! Input argument
52    REAL,DIMENSION(klon,nbtr), INTENT(IN) :: trs_in
53   
54    ! Local variables
55    INTEGER :: ierr
56   
57    ! Allocate restart variables trs
58    ALLOCATE( trs(klon,nbtr), stat=ierr)
59    IF (ierr /= 0) CALL abort_gcm('traclmdz_from_restart', 'pb in allocation 1',1)
60   
61    ! Initialize trs with values read from restart file
62    trs(:,:) = trs_in(:,:)
63   
64  END SUBROUTINE traclmdz_from_restart
65
66
67  SUBROUTINE traclmdz_init(pctsrf, ftsol, tr_seri, aerosol, lessivage)
68    ! This subroutine allocates and initialize module variables and control variables.
69    USE dimphy
70    USE infotrac
71    USE carbon_cycle_mod, ONLY : carbon_cycle_init, carbon_cycle_tr, carbon_cycle_cpl
72
73    IMPLICIT NONE
74
75    INCLUDE "indicesol.h"
76
77! Input variables
78    REAL,DIMENSION(klon,nbsrf),INTENT(IN)     :: pctsrf ! Pourcentage de sol f(nature du sol)
79    REAL,DIMENSION(klon,nbsrf),INTENT(IN)     :: ftsol  ! Temperature du sol (surf)(Kelvin)
80    REAL,DIMENSION(klon,klev,nbtr),INTENT(IN) :: tr_seri! Concentration Traceur [U/KgA] 
81
82! Output variables
83    LOGICAL,DIMENSION(nbtr), INTENT(OUT) :: aerosol
84    LOGICAL,INTENT(OUT)                  :: lessivage
85       
86! Local variables   
87    INTEGER :: ierr, it, iiq
88   
89! --------------------------------------------
90! Allocation
91! --------------------------------------------
92
93    ALLOCATE( scavtr(nbtr), stat=ierr)
94    IF (ierr /= 0) CALL abort_gcm('traclmdz_init', 'pb in allocation 9',1)
95    scavtr(:)=1.
96   
97    ALLOCATE( radio(nbtr), stat=ierr)
98    IF (ierr /= 0) CALL abort_gcm('traclmdz_init', 'pb in allocation 11',1)
99    radio(:) = .false.    ! Par defaut pas decroissance radioactive
100   
101    ALLOCATE( masktr(klon,nbtr), stat=ierr)
102    IF (ierr /= 0) CALL abort_gcm('traclmdz_init', 'pb in allocation 2',1)
103   
104    ALLOCATE( fshtr(klon,nbtr), stat=ierr)
105    IF (ierr /= 0) CALL abort_gcm('traclmdz_init', 'pb in allocation 3',1)
106   
107    ALLOCATE( hsoltr(nbtr), stat=ierr)
108    IF (ierr /= 0) CALL abort_gcm('traclmdz_init', 'pb in allocation 4',1)
109   
110    ALLOCATE( tautr(nbtr), stat=ierr)
111    IF (ierr /= 0) CALL abort_gcm('traclmdz_init', 'pb in allocation 5',1)
112    tautr(:)  = 0.
113   
114    ALLOCATE( vdeptr(nbtr), stat=ierr)
115    IF (ierr /= 0) CALL abort_gcm('traclmdz_init', 'pb in allocation 6',1)
116    vdeptr(:) = 0.
117
118
119    lessivage  = .TRUE.
120    aerosol(:) = .FALSE.  ! Tous les traceurs sont des gaz par defaut
121   
122!
123! Recherche des traceurs connus : Be7, CO2,...
124! --------------------------------------------
125    id_be=0
126    DO it=1,nbtr
127       iiq=niadv(it+2)
128       IF ( tname(iiq) == "BE" .OR. tname(iiq) == "Be" .OR.  &
129            tname(iiq) == "BE7" .OR. tname(iiq) == "Be7" ) THEN 
130          ! Recherche du Beryllium 7
131          id_be=it
132          ALLOCATE( srcbe(klon,klev) )
133          radio(id_be) = .TRUE.
134          aerosol(id_be) = .TRUE. ! le Be est un aerosol
135          CALL init_be(pctsrf,masktr(:,id_be),tautr(id_be),vdeptr(id_be),scavtr(id_be),srcbe)
136          WRITE(*,*) 'Initialisation srcBe: OK'
137       END IF   
138    END DO
139!
140! Valeurs specifiques pour les traceurs Rn222 et Pb210
141! ----------------------------------------------
142    IF (rnpb) THEN
143       
144       radio(1)= .TRUE.
145       radio(2)= .TRUE.
146       pbl_flg(1) = 0 ! au lieu de clsol=true ! CL au sol calcule
147       pbl_flg(2) = 0 ! au lieu de clsol=true
148       
149       aerosol(2) = .TRUE. ! le Pb est un aerosol
150       
151       CALL initrrnpb (ftsol,pctsrf,masktr,fshtr,hsoltr,tautr,vdeptr,scavtr)
152    END IF
153
154!
155! Initialisation de module carbon_cycle_mod
156! ----------------------------------------------
157    IF (carbon_cycle_tr .OR. carbon_cycle_cpl) THEN
158       CALL carbon_cycle_init(tr_seri, aerosol, radio)
159    END IF
160
161  END SUBROUTINE traclmdz_init
162
163  SUBROUTINE traclmdz(                           &
164       nstep,    pdtphys,      t_seri,           &
165       paprs,    pplay,        cdragh,  coefh,   &
166       yu1,      yv1,          ftsol,   pctsrf,  &
167       xlat,     couchelimite,                   &
168       tr_seri,  source,       solsym,  d_tr_cl)
169   
170    USE dimphy
171    USE infotrac
172    USE carbon_cycle_mod, ONLY : carbon_cycle, carbon_cycle_tr, carbon_cycle_cpl
173   
174    IMPLICIT NONE
175   
176    INCLUDE "YOMCST.h"
177    INCLUDE "indicesol.h"
178
179!==========================================================================
180!                   -- DESCRIPTION DES ARGUMENTS --
181!==========================================================================
182
183! Input arguments
184!
185!Configuration grille,temps:
186    INTEGER,INTENT(IN) :: nstep      ! nombre d'appels de la physiq
187    REAL,INTENT(IN)    :: pdtphys    ! Pas d'integration pour la physique (seconde) 
188    REAL,DIMENSION(klon),INTENT(IN) :: xlat    ! latitudes pour chaque point
189
190!
191!Physique:
192!--------
193    REAL,DIMENSION(klon,klev),INTENT(IN)   :: t_seri  ! Temperature
194    REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs   ! pression pour chaque inter-couche (en Pa)
195    REAL,DIMENSION(klon,klev),INTENT(IN)   :: pplay   ! pression pour le mileu de chaque couche (en Pa)
196
197
198!Couche limite:
199!--------------
200!
201    REAL,DIMENSION(klon,klev),INTENT(IN) :: cdragh     ! coeff drag pour T et Q
202    REAL,DIMENSION(klon,klev),INTENT(IN) :: coefh      ! coeff melange CL (m**2/s)
203    REAL,DIMENSION(klon),INTENT(IN)      :: yu1        ! vents au premier niveau
204    REAL,DIMENSION(klon),INTENT(IN)      :: yv1        ! vents au premier niveau
205    LOGICAL,INTENT(IN)                   :: couchelimite
206
207! Arguments necessaires pour les sources et puits de traceur:
208    REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: ftsol  ! Temperature du sol (surf)(Kelvin)
209    REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: pctsrf ! Pourcentage de sol f(nature du sol)
210
211
212! InOutput argument
213    REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT)  :: tr_seri ! Concentration Traceur [U/KgA] 
214
215! Output argument
216    CHARACTER(len=8),DIMENSION(nbtr), INTENT(OUT) :: solsym
217    REAL,DIMENSION(klon,nbtr), INTENT(OUT)        :: source  ! a voir lorsque le flux de surface est prescrit
218    REAL,DIMENSION(klon,klev,nbtr), INTENT(OUT)   :: d_tr_cl ! Td couche limite/traceur
219
220!=======================================================================================
221!                        -- VARIABLES LOCALES TRACEURS --
222!=======================================================================================
223
224    INTEGER :: i, k, it
225
226    REAL,DIMENSION(klon)           :: d_trs    ! Td dans le reservoir
227    REAL,DIMENSION(klon,klev)      :: delp     ! epaisseur de couche (Pa)
228   
229    REAL,DIMENSION(klon,klev,nbtr) :: d_tr_dec ! Td radioactive
230    REAL                           :: zrho      ! Masse Volumique de l'air KgA/m3
231
232!
233!
234!=================================================================
235!  Ajout de la production en  Be7 (Beryllium) srcbe U/s/kgA
236!=================================================================
237!
238    IF ( id_be /= 0 ) THEN
239       DO k = 1, klev
240          DO i = 1, klon
241             tr_seri(i,k,id_be) = tr_seri(i,k,id_be)+srcbe(i,k)*pdtphys
242          END DO
243       END DO
244       WRITE(*,*) 'Ajout srcBe dans tr_seri: OK'
245    END IF
246 
247
248    DO it=1,nbtr
249       WRITE(solsym(it),'(i2)') it
250    END DO
251!======================================================================
252!     -- Calcul de l'effet de la couche limite --
253!======================================================================
254
255    IF (couchelimite) THEN             
256       source(:,:) = 0.0
257
258       IF (id_be /=0) THEN
259          DO i=1, klon
260             zrho = pplay(i,1)/t_seri(i,1)/RD
261             source(i,id_be) = - vdeptr(id_be)*tr_seri(i,1,id_be)*zrho
262          END DO
263       END IF
264
265    END IF
266   
267
268    DO k = 1, klev
269       DO i = 1, klon
270          delp(i,k) = paprs(i,k)-paprs(i,k+1)
271       END DO
272    END DO
273   
274    DO it=1, nbtr
275       IF (couchelimite .AND. pbl_flg(it) == 0 ) THEN ! couche limite avec quantite dans le sol calculee
276          CALL cltracrn(it, pdtphys, yu1, yv1,     &
277               cdragh, coefh,t_seri,ftsol,pctsrf,  &
278               tr_seri(:,:,it),trs(:,it),          &
279               paprs, pplay, delp,                 &
280               masktr(:,it),fshtr(:,it),hsoltr(it),&
281               tautr(it),vdeptr(it),               &
282               xlat,d_tr_cl(:,:,it),d_trs)
283         
284          DO k = 1, klev
285             DO i = 1, klon
286                tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cl(i,k,it)
287             END DO
288          END DO
289       
290          ! Traceur dans le reservoir sol
291          DO i = 1, klon
292             trs(i,it) = trs(i,it) + d_trs(i)
293          END DO
294       END IF
295    END DO
296           
297!======================================================================
298!   Calcul de l'effet du puits radioactif
299!======================================================================
300    CALL radio_decay (radio,rnpb,pdtphys,tautr,tr_seri,d_tr_dec)
301 
302    DO it=1,nbtr
303       IF(radio(it)) then     
304          DO k = 1, klev
305             DO i = 1, klon
306                tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_dec(i,k,it)
307             END DO
308          END DO
309          CALL minmaxqfi(tr_seri(:,:,it),0.,1.e33,'puits rn it='//solsym(it))
310       END IF
311    END DO
312
313!======================================================================
314!   Calcul de cycle de carbon
315!======================================================================
316    IF (carbon_cycle_tr .OR. carbon_cycle_cpl) THEN
317       CALL carbon_cycle(nstep, pdtphys, pctsrf, tr_seri)
318    END IF
319
320  END SUBROUTINE traclmdz
321
322
323  SUBROUTINE traclmdz_to_restart(trs_out)
324    ! This subroutine is called from phyredem.F where the module
325    ! variable trs is written to restart file (restartphy.nc)
326    USE dimphy
327    USE infotrac
328   
329    IMPLICIT NONE
330   
331    REAL,DIMENSION(klon,nbtr), INTENT(OUT) :: trs_out
332    INTEGER :: ierr
333
334    IF ( ALLOCATED(trs) ) THEN
335       trs_out(:,:) = trs(:,:)
336    ELSE
337       ! No previous allocate of trs. This is the case for create_etat0_limit.
338       trs_out(:,:) = 0.0
339    END IF
340   
341  END SUBROUTINE traclmdz_to_restart
342 
343
344END MODULE traclmdz_mod
Note: See TracBrowser for help on using the repository browser.