Context Navigation

traclmdz_mod.F90 @ 5151

Last change on this file since 5151 was 5144, checked in by abarral, 7 weeks ago
Put YOMCST.h into 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.1 KB

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

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

Download in other formats:

Original Format