Context Navigation

traclmdz_mod.F90 @ 3816

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

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

Download in other formats:

Original Format