Context Navigation

traclmdz_mod.F90 @ 5112

Last change on this file since 5112 was 5112, checked in by abarral, 6 months ago
Rename modules in phy_common from _mod > lmdz_
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.0 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 strings_mod, 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
315	INCLUDE "YOMCST.h"
316
317	!==========================================================================
318	! -- DESCRIPTION DES ARGUMENTS --
319	!==========================================================================
320
321	! Input arguments
322
323	!Configuration grille,temps:
324	INTEGER, INTENT(IN) :: nstep ! nombre d'appels de la physiq
325	INTEGER, INTENT(IN) :: julien ! Jour julien
326	REAL, INTENT(IN) :: gmtime
327	REAL, INTENT(IN) :: pdtphys ! Pas d'integration pour la physique (seconde)
328	REAL, DIMENSION(klon), INTENT(IN) :: xlat ! latitudes pour chaque point
329	REAL, INTENT(IN) :: xlon(:) ! dim(klon) longitude
330
331	!Physique:
332	!--------
333	REAL, DIMENSION(klon, klev), INTENT(IN) :: t_seri ! Temperature
334	REAL, DIMENSION(klon, klev + 1), INTENT(IN) :: paprs ! pression pour chaque inter-couche (en Pa)
335	REAL, DIMENSION(klon, klev), INTENT(IN) :: pplay ! pression pour le mileu de chaque couche (en Pa)
336	REAL, intent(in) :: zmasse (:, :) ! dim(klon,klev) density of air, in kg/m2
337
338
339	!Couche limite:
340	!--------------
341
342	REAL, DIMENSION(klon), INTENT(IN) :: cdragh ! coeff drag pour T et Q
343	REAL, DIMENSION(klon, klev), INTENT(IN) :: coefh ! coeff melange CL (m**2/s)
344	REAL, DIMENSION(klon), INTENT(IN) :: yu1 ! vents au premier niveau
345	REAL, DIMENSION(klon), INTENT(IN) :: yv1 ! vents au premier niveau
346	LOGICAL, INTENT(IN) :: couchelimite
347	REAL, DIMENSION(klon, klev), INTENT(IN) :: sh ! humidite specifique
348	REAL, DIMENSION(klon, klev), INTENT(IN) :: rh ! Humidite relative
349	REAL, DIMENSION(klon, klev), INTENT(IN) :: pphi ! geopotentie
350	REAL, DIMENSION(klon), INTENT(IN) :: ustar ! ustar (m/s)
351	REAL, DIMENSION(klon), INTENT(IN) :: wstar, ale_bl, ale_wake ! wstar (m/s) and Avail. Lifti. Energ.
352	REAL, DIMENSION(klon), INTENT(IN) :: zu10m ! vent zonal 10m (m/s)
353	REAL, DIMENSION(klon), INTENT(IN) :: zv10m ! vent zonal 10m (m/s)
354
355	! Arguments necessaires pour les sources et puits de traceur:
356	REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: ftsol ! Temperature du sol (surf)(Kelvin)
357	REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: pctsrf ! Pourcentage de sol f(nature du sol)
358
359	! InOutput argument
360	REAL, DIMENSION(klon, klev, nbtr), INTENT(INOUT) :: tr_seri ! Concentration Traceur [U/KgA]
361
362	! Output argument
363	REAL, DIMENSION(klon, nbtr), INTENT(OUT) :: source ! a voir lorsque le flux de surface est prescrit
364	REAL, DIMENSION(klon, klev, nbtr), INTENT(OUT) :: d_tr_cl ! Td couche limite/traceur
365
366	!=======================================================================================
367	! -- VARIABLES LOCALES TRACEURS --
368	!=======================================================================================
369
370	INTEGER :: i, k, it
371	INTEGER :: lmt_pas ! number of time steps of "physics" per day
372
373	REAL, DIMENSION(klon) :: d_trs ! Td dans le reservoir
374	REAL, DIMENSION(klon, klev) :: qsat ! pression de la vapeur a saturation
375	REAL, DIMENSION(klon, klev, nbtr) :: d_tr_dec ! Td radioactive
376	REAL :: zrho ! Masse Volumique de l'air KgA/m3
377	REAL :: amn, amx
378
379	!=================================================================
380	! Ajout de la production en Be7 (Beryllium) srcbe U/s/kgA
381	!=================================================================
382
383	IF (id_be /= 0) THEN
384	DO k = 1, klev
385	DO i = 1, klon
386	tr_seri(i, k, id_be) = tr_seri(i, k, id_be) + srcbe(i, k) * pdtphys
387	END DO
388	END DO
389	WRITE(, ) 'Ajout srcBe dans tr_seri: OK'
390	END IF
391
392
393	!=================================================================
394	! Update pseudo-vapor tracers
395	!=================================================================
396
397	CALL q_sat(klon * klev, t_seri, pplay, qsat)
398
399	IF (id_pcsat /= 0) THEN
400	DO k = 1, klev
401	DO i = 1, klon
402	IF (pplay(i, k)>=85000.) THEN
403	tr_seri(i, k, id_pcsat) = qsat(i, k)
404	ELSE
405	tr_seri(i, k, id_pcsat) = MIN (qsat(i, k), tr_seri(i, k, id_pcsat))
406	END IF
407	END DO
408	END DO
409	END IF
410
411	IF (id_pcocsat /= 0) THEN
412	DO k = 1, klev
413	DO i = 1, klon
414	IF (pplay(i, k)>=85000.) THEN
415	IF (pctsrf (i, is_oce) > 0.) THEN
416	tr_seri(i, k, id_pcocsat) = qsat(i, k)
417	ELSE
418	tr_seri(i, k, id_pcocsat) = 0.
419	END IF
420	ELSE
421	tr_seri(i, k, id_pcocsat) = MIN (qsat(i, k), tr_seri(i, k, id_pcocsat))
422	END IF
423	END DO
424	END DO
425	END IF
426
427	IF (id_pcq /= 0) THEN
428	DO k = 1, klev
429	DO i = 1, klon
430	IF (pplay(i, k)>=85000.) THEN
431	tr_seri(i, k, id_pcq) = sh(i, k)
432	ELSE
433	tr_seri(i, k, id_pcq) = MIN (qsat(i, k), tr_seri(i, k, id_pcq))
434	END IF
435	END DO
436	END DO
437	END IF
438
439	IF (id_pcs0 /= 0) THEN
440	DO k = 1, klev
441	DO i = 1, klon
442	IF (pplay(i, k)>=85000.) THEN
443	tr_seri(i, k, id_pcs0) = qsat(i, k)
444	ELSE
445	tr_seri(i, k, id_pcs0) = MIN (qsat(i, k), tr_seri(i, k, id_pcs0))
446	END IF
447	END DO
448	END DO
449	END IF
450
451	IF (id_pcos0 /= 0) THEN
452	DO k = 1, klev
453	DO i = 1, klon
454	IF (pplay(i, k)>=85000.) THEN
455	IF (pctsrf (i, is_oce) > 0.) THEN
456	tr_seri(i, k, id_pcos0) = qsat(i, k)
457	ELSE
458	tr_seri(i, k, id_pcos0) = 0.
459	END IF
460	ELSE
461	tr_seri(i, k, id_pcos0) = MIN (qsat(i, k), tr_seri(i, k, id_pcos0))
462	END IF
463	END DO
464	END DO
465	END IF
466
467	IF (id_pcq0 /= 0) THEN
468	DO k = 1, klev
469	DO i = 1, klon
470	IF (pplay(i, k)>=85000.) THEN
471	tr_seri(i, k, id_pcq0) = sh(i, k)
472	ELSE
473	tr_seri(i, k, id_pcq0) = MIN (qsat(i, k), tr_seri(i, k, id_pcq0))
474	END IF
475	END DO
476	END DO
477	END IF
478
479	!=================================================================
480	! Update tracer : Age of stratospheric air
481	!=================================================================
482	IF (id_aga/=0) THEN
483
484	! Bottom layers
485	DO k = 1, lev_1p5km
486	tr_seri(:, k, id_aga) = 0.0
487	END DO
488
489	! Layers above 1.5km
490	DO k = lev_1p5km + 1, klev - 1
491	tr_seri(:, k, id_aga) = tr_seri(:, k, id_aga) + pdtphys
492	END DO
493
494	! Top layer
495	tr_seri(:, klev, id_aga) = tr_seri(:, klev - 1, id_aga)
496
497	END IF
498
499	!======================================================================
500	! -- Calcul de l'effet de la couche limite --
501	!======================================================================
502
503	IF (couchelimite) THEN
504	source(:, :) = 0.0
505
506	IF (id_be /=0) THEN
507	DO i = 1, klon
508	zrho = pplay(i, 1) / t_seri(i, 1) / RD
509	source(i, id_be) = - vdeptr(id_be) * tr_seri(i, 1, id_be) * zrho
510	END DO
511	END IF
512
513	END IF
514
515	DO it = 1, nbtr
516	IF (couchelimite .AND. pbl_flg(it) == 0 .AND. (it==id_rn .OR. it==id_pb)) THEN
517	! couche limite avec quantite dans le sol calculee
518	CALL cltracrn(it, pdtphys, yu1, yv1, &
519	cdragh, coefh, t_seri, ftsol, pctsrf, &
520	tr_seri(:, :, it), trs(:, it), &
521	paprs, pplay, zmasse * rg, &
522	masktr(:, it), fshtr(:, it), hsoltr(it), &
523	tautr(it), vdeptr(it), &
524	xlat, d_tr_cl(:, :, it), d_trs)
525
526	DO k = 1, klev
527	DO i = 1, klon
528	tr_seri(i, k, it) = tr_seri(i, k, it) + d_tr_cl(i, k, it)
529	END DO
530	END DO
531
532	! Traceur dans le reservoir sol
533	DO i = 1, klon
534	trs(i, it) = trs(i, it) + d_trs(i)
535	END DO
536	END IF
537	END DO
538
539
540	!======================================================================
541	! Calcul de l'effet du puits radioactif
542	!======================================================================
543	CALL radio_decay (radio, rnpb, pdtphys, tautr, tr_seri, d_tr_dec)
544
545	DO it = 1, nbtr
546	IF(radio(it)) then
547	DO k = 1, klev
548	DO i = 1, klon
549	tr_seri(i, k, it) = tr_seri(i, k, it) + d_tr_dec(i, k, it)
550	END DO
551	END DO
552	CALL minmaxqfi(tr_seri(:, :, it), 0., 1.e33, 'puits rn it=' // TRIM(int2str(it)))
553	END IF
554	END DO
555
556	!======================================================================
557	! Parameterization of ozone chemistry
558	!======================================================================
559
560	IF (id_o3 /= 0) then
561	lmt_pas = NINT(86400. / pdtphys)
562	IF (MOD(nstep - 1, lmt_pas) == 0) THEN
563	! Once per day, update the coefficients for ozone chemistry:
564	CALL regr_pr_comb_coefoz(julien, xlat, paprs, pplay)
565	END IF
566	CALL o3_chem(julien, gmtime, t_seri, zmasse, pdtphys, xlat, &
567	xlon, tr_seri(:, :, id_o3))
568	END IF
569
570	END SUBROUTINE traclmdz
571
572
573	SUBROUTINE traclmdz_to_restart(trs_out)
574	! This SUBROUTINE is called from phyredem.F where the module
575	! variable trs is written to restart file (restartphy.nc)
576	USE dimphy
577	USE infotrac_phy, ONLY: nbtr
578
579	REAL, DIMENSION(klon, nbtr), INTENT(OUT) :: trs_out
580	INTEGER :: ierr
581
582	IF (ALLOCATED(trs)) THEN
583	trs_out(:, :) = trs(:, :)
584	ELSE
585	! No previous allocate of trs. This is the case for create_etat0_limit.
586	trs_out(:, :) = 0.0
587	END IF
588
589	END SUBROUTINE traclmdz_to_restart
590
591
592	END MODULE traclmdz_mod

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

Download in other formats:

Original Format