Context Navigation

traclmdz_mod.F90 @ 5099

Last change on this file since 5099 was 5099, checked in by abarral, 4 months ago
Replace most uses of CPP_DUST by the corresponding logical defined in lmdz_cppkeys_wrapper.F90 Convert several files from .F to .f90 to allow Dust to compile w/o rrtm/ecrad Create lmdz_yoerad.f90 (lint) Remove "!" on otherwise empty line
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
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 mod_grid_phy_lmdz
94	USE mod_phys_lmdz_para
95	USE indice_sol_mod
96	USE print_control_mod, 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
158	lessivage = .TRUE.
159	!!jyg(20130206) : le choix d activation du lessivage est fait dans phytrac avec iflag_lscav
160	!! call getin('lessivage',lessivage)
161	!! if(lessivage) then
162	!! print*,'lessivage lsc ON'
163	!! else
164	!! print*,'lessivage lsc OFF'
165	!! endif
166	aerosol(:) = .FALSE. ! Tous les traceurs sont des gaz par defaut
167
168	! Recherche des traceurs connus : Be7, O3, CO2,...
169	! --------------------------------------------
170	id_rn=0; id_pb=0; id_aga=0; id_be=0; id_o3=0
171	id_pcsat=0; id_pcocsat=0; id_pcq=0; id_pcs0=0; id_pcos0=0; id_pcq0=0
172	it = 0
173	DO iq = 1, nqtot
174	IF(.NOT.(tracers(iq)%isInPhysics)) CYCLE
175	it = it+1
176	SELECT CASE(tracers(iq)%name)
177	CASE("RN"); id_rn = it ! radon
178	CASE("PB"); id_pb = it ! plomb
179	CASE("Aga","AGA"); id_aga = it ! Age of stratospheric air
180	CASE("Be","BE","Be7","BE7"); id_be = it ! Recherche du Beryllium 7
181	CASE("o3","O3"); id_o3 = it ! Recherche de l'ozone
182	CASE("pcsat", "Pcsat"); id_pcsat = it
183	CASE("pcocsat","Pcocsat"); id_pcocsat= it
184	CASE("pcq", "Pcq"); id_pcq = it
185	CASE("pcs0", "Pcs0"); id_pcs0 = it
186	CASE("pcos0", "Pcos0"); id_pcos0 = it
187	CASE("pcq0", "Pcq0"); id_pcq0 = it
188	CASE DEFAULT
189	WRITE(lunout,*) 'This is an unknown tracer in LMDZ : ', trim(tracers(iq)%name)
190	END SELECT
191
192	SELECT CASE(tracers(iq)%name)
193	CASE("PB") !--- RomP >>> profil initial de PB210
194	OPEN(ilesfil2,file='prof.pb210',status='old',iostat=irr2)
195	IF(irr2 == 0) THEN
196	READ(ilesfil2,*) kradio2
197	WRITE(lunout,*)'number of levels for pb210 profile ',kradio2
198	DO k=kradio2,1,-1; READ (ilesfil2,*) plomb(:,k); 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); END DO
228	CLOSE(ilesfil)
229	tr_seri(:,:,id_be)=beryllium(:,:)
230	ELSE
231	WRITE(lunout,*)'Prof. Be7 does not exist: use restart values'
232	END IF
233	CASE("o3","O3") !--- Parametrisation par la chimie de Cariolle
234	CALL alloc_coefoz !--- Allocate ozone coefficients
235	CALL press_coefoz !--- Read input pressure levels
236	CASE("pcs0","Pcs0", "pcos0","Pcos0", "pcq0","Pcq0")
237	conv_flg(it)=0 !--- No transport by convection for this tracer
238	END SELECT
239	END DO
240
241	! Valeurs specifiques pour les traceurs Rn222 et Pb210
242	! ----------------------------------------------
243	IF ( id_rn/=0 .AND. id_pb/=0 ) THEN
244	rnpb = .TRUE.
245	radio(id_rn)= .TRUE.
246	radio(id_pb)= .TRUE.
247	pbl_flg(id_rn) = 0 ! au lieu de clsol=true ! CL au sol calcule
248	pbl_flg(id_pb) = 0 ! au lieu de clsol=true
249	aerosol(id_rn) = .FALSE.
250	aerosol(id_pb) = .TRUE. ! le Pb est un aerosol
251
252	CALL initrrnpb (ftsol,pctsrf,masktr,fshtr,hsoltr,tautr,vdeptr,scavtr)
253	END IF
254
255	! Check if all tracers have restart values
256	! ----------------------------------------------
257	it = 0
258	DO iq = 1, nqtot
259	IF(.NOT.(tracers(iq)%isAdvected .AND. tracers(iq)%isInPhysics)) CYCLE
260	it = it+1
261	! Test if tracer is zero everywhere.
262	! Done by master process MPI and master thread OpenMP
263	CALL gather(tr_seri(:,:,it),varglo)
264	IF (is_mpi_root .AND. is_omp_root) THEN
265	mintmp=MINVAL(varglo,dim=1)
266	maxtmp=MAXVAL(varglo,dim=1)
267	IF (MINVAL(mintmp,dim=1)==0. .AND. MAXVAL(maxtmp,dim=1)==0.) THEN
268	! Tracer is zero everywhere
269	zero=.TRUE.
270	ELSE
271	zero=.FALSE.
272	END IF
273	END IF
274
275	! Distribute variable at all processes
276	CALL bcast(zero)
277
278	! Initalize tracer that was not found in restart file.
279	IF (zero) THEN
280	! The tracer was not found in restart file or it was equal zero everywhere.
281	WRITE(lunout,*) "The tracer ",trim(tracers(iq)%name)," will be initialized"
282	IF (it==id_pcsat .OR. it==id_pcq .OR. &
283	it==id_pcs0 .OR. it==id_pcq0) THEN
284	tr_seri(:,:,it) = 100.
285	ELSE IF (it==id_pcocsat .OR. it==id_pcos0) THEN
286	DO i = 1, klon
287	IF ( pctsrf (i, is_oce) == 0. ) THEN
288	tr_seri(i,:,it) = 0.
289	ELSE
290	tr_seri(i,:,it) = 100.
291	END IF
292	END DO
293	ELSE
294	! No specific initialization exist for this tracer
295	tr_seri(:,:,it) = 0.
296	END IF
297	END IF
298	END DO
299
300	END SUBROUTINE traclmdz_init
301
302	SUBROUTINE traclmdz(nstep, julien, gmtime, pdtphys, t_seri, paprs, pplay, &
303	cdragh, coefh, yu1, yv1, ftsol, pctsrf, xlat, xlon, couchelimite, sh, &
304	rh, pphi, ustar, wstar, ale_bl, ale_wake, zu10m, zv10m, &
305	tr_seri, source, d_tr_cl,d_tr_dec, zmasse) !RomP
306
307	USE dimphy
308	USE infotrac_phy, ONLY: nbtr, pbl_flg
309	USE strings_mod, ONLY: int2str
310	USE regr_pr_comb_coefoz_m, ONLY: regr_pr_comb_coefoz
311	USE o3_chem_m, ONLY: o3_chem
312	USE indice_sol_mod
313
314	INCLUDE "YOMCST.h"
315
316	!==========================================================================
317	! -- DESCRIPTION DES ARGUMENTS --
318	!==========================================================================
319
320	! Input arguments
321
322	!Configuration grille,temps:
323	INTEGER,INTENT(IN) :: nstep ! nombre d'appels de la physiq
324	INTEGER,INTENT(IN) :: julien ! Jour julien
325	REAL,INTENT(IN) :: gmtime
326	REAL,INTENT(IN) :: pdtphys ! Pas d'integration pour la physique (seconde)
327	REAL,DIMENSION(klon),INTENT(IN) :: xlat ! latitudes pour chaque point
328	REAL, INTENT(IN):: xlon(:) ! dim(klon) longitude
329
330	!Physique:
331	!--------
332	REAL,DIMENSION(klon,klev),INTENT(IN) :: t_seri ! Temperature
333	REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs ! pression pour chaque inter-couche (en Pa)
334	REAL,DIMENSION(klon,klev),INTENT(IN) :: pplay ! pression pour le mileu de chaque couche (en Pa)
335	REAL,intent(in):: zmasse (:, :) ! dim(klon,klev) density of air, in kg/m2
336
337
338	!Couche limite:
339	!--------------
340
341	REAL,DIMENSION(klon),INTENT(IN) :: cdragh ! coeff drag pour T et Q
342	REAL,DIMENSION(klon,klev),INTENT(IN) :: coefh ! coeff melange CL (m**2/s)
343	REAL,DIMENSION(klon),INTENT(IN) :: yu1 ! vents au premier niveau
344	REAL,DIMENSION(klon),INTENT(IN) :: yv1 ! vents au premier niveau
345	LOGICAL,INTENT(IN) :: couchelimite
346	REAL,DIMENSION(klon,klev),INTENT(IN) :: sh ! humidite specifique
347	REAL,DIMENSION(klon,klev),INTENT(IN) :: rh ! Humidite relative
348	REAL,DIMENSION(klon,klev),INTENT(IN) :: pphi ! geopotentie
349	REAL,DIMENSION(klon),INTENT(IN) :: ustar ! ustar (m/s)
350	REAL,DIMENSION(klon),INTENT(IN) :: wstar,ale_bl,ale_wake ! wstar (m/s) and Avail. Lifti. Energ.
351	REAL,DIMENSION(klon),INTENT(IN) :: zu10m ! vent zonal 10m (m/s)
352	REAL,DIMENSION(klon),INTENT(IN) :: zv10m ! vent zonal 10m (m/s)
353
354	! Arguments necessaires pour les sources et puits de traceur:
355	REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: ftsol ! Temperature du sol (surf)(Kelvin)
356	REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: pctsrf ! Pourcentage de sol f(nature du sol)
357
358	! InOutput argument
359	REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT) :: tr_seri ! Concentration Traceur [U/KgA]
360
361	! Output argument
362	REAL,DIMENSION(klon,nbtr), INTENT(OUT) :: source ! a voir lorsque le flux de surface est prescrit
363	REAL,DIMENSION(klon,klev,nbtr), INTENT(OUT) :: d_tr_cl ! Td couche limite/traceur
364
365	!=======================================================================================
366	! -- VARIABLES LOCALES TRACEURS --
367	!=======================================================================================
368
369	INTEGER :: i, k, it
370	INTEGER :: lmt_pas ! number of time steps of "physics" per day
371
372	REAL,DIMENSION(klon) :: d_trs ! Td dans le reservoir
373	REAL,DIMENSION(klon,klev) :: qsat ! pression de la vapeur a saturation
374	REAL,DIMENSION(klon,klev,nbtr) :: d_tr_dec ! Td radioactive
375	REAL :: zrho ! Masse Volumique de l'air KgA/m3
376	REAL :: amn, amx
377
378	!=================================================================
379	! Ajout de la production en Be7 (Beryllium) srcbe U/s/kgA
380	!=================================================================
381
382	IF ( id_be /= 0 ) THEN
383	DO k = 1, klev
384	DO i = 1, klon
385	tr_seri(i,k,id_be) = tr_seri(i,k,id_be)+srcbe(i,k)*pdtphys
386	END DO
387	END DO
388	WRITE(,) 'Ajout srcBe dans tr_seri: OK'
389	END IF
390
391
392	!=================================================================
393	! Update pseudo-vapor tracers
394	!=================================================================
395
396	CALL q_sat(klon*klev,t_seri,pplay,qsat)
397
398	IF ( id_pcsat /= 0 ) THEN
399	DO k = 1, klev
400	DO i = 1, klon
401	IF ( pplay(i,k)>=85000.) THEN
402	tr_seri(i,k,id_pcsat) = qsat(i,k)
403	ELSE
404	tr_seri(i,k,id_pcsat) = MIN (qsat(i,k), tr_seri(i,k,id_pcsat))
405	END IF
406	END DO
407	END DO
408	END IF
409
410	IF ( id_pcocsat /= 0 ) THEN
411	DO k = 1, klev
412	DO i = 1, klon
413	IF ( pplay(i,k)>=85000.) THEN
414	IF ( pctsrf (i, is_oce) > 0. ) THEN
415	tr_seri(i,k,id_pcocsat) = qsat(i,k)
416	ELSE
417	tr_seri(i,k,id_pcocsat) = 0.
418	END IF
419	ELSE
420	tr_seri(i,k,id_pcocsat) = MIN (qsat(i,k), tr_seri(i,k,id_pcocsat))
421	END IF
422	END DO
423	END DO
424	END IF
425
426	IF ( id_pcq /= 0 ) THEN
427	DO k = 1, klev
428	DO i = 1, klon
429	IF ( pplay(i,k)>=85000.) THEN
430	tr_seri(i,k,id_pcq) = sh(i,k)
431	ELSE
432	tr_seri(i,k,id_pcq) = MIN (qsat(i,k), tr_seri(i,k,id_pcq))
433	END IF
434	END DO
435	END DO
436	END IF
437
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
452	IF ( id_pcos0 /= 0 ) THEN
453	DO k = 1, klev
454	DO i = 1, klon
455	IF ( pplay(i,k)>=85000.) THEN
456	IF ( pctsrf (i, is_oce) > 0. ) THEN
457	tr_seri(i,k,id_pcos0) = qsat(i,k)
458	ELSE
459	tr_seri(i,k,id_pcos0) = 0.
460	END IF
461	ELSE
462	tr_seri(i,k,id_pcos0) = MIN (qsat(i,k), tr_seri(i,k,id_pcos0))
463	END IF
464	END DO
465	END DO
466	END IF
467
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