Context Navigation

traclmdz_mod.F90 @ 3993

Last change on this file since 3993 was 3852, checked in by dcugnet, 4 years ago

Extension of the tracers management.

The tracers files can be:

1) "traceur.def": old format, with:

the number of tracers on the first line
one line for each tracer: <tracer name> <hadv> <vadv> [<parent name>]

2) "tracer.def": new format with one section each model component.
3) "tracer_<name>.def": new format with a single section.

The formats 2 and 3 reading is driven by the "type_trac" key, which can be a

coma-separated list of components.

Format 2: read the sections from the "tracer.def" file.
format 3: read one section each "tracer_<section name>.def" file.
the first line of a section is "&<section name>
the other lines start with a tracer name followed by <key>=<val> pairs.
the "default" tracer name is reserved ; the other tracers of the section inherit its <key>=<val>, except for the keys that are redefined locally.

This format helps keeping the tracers files compact, thanks to the "default"
special tracer and the three levels of factorization:

on the tracers names: a tracer name can be a coma-separated list of tracers => all the tracers of the list have the same <key>=<val> properties
on the parents names: the value of the "parent" property can be a coma-separated list of tracers => only possible for geographic tagging tracers
on the phases: the property "phases" is [g](l][s] (gas/liquid/solid)

Read information is stored in the vector "tracers(:)", of derived type "tra".

"isotopes_params.def" is a similar file, with one section each isotopes family.
It contains a database of isotopes properties ; if there are second generation
tracers (isotopes), the corresponding sections are read.

Read information is stored in the vector "isotopes(:)", of derived type "iso".

The "getKey" function helps to get the values of the parameters stored in
"tracers" or "isotopes".

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

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

Download in other formats:

Original Format