Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

traclmdz_mod.F90 @ 5452

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

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

Download in other formats:

Original Format