Context Navigation

traclmdz_mod.f90 @ 5452

Last change on this file since 5452 was 5285, checked in by abarral, 2 months ago
As discussed internally, remove generic ONLY: ... for new _mod_h modules
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory
File size: 22.2 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format