Context Navigation

traclmdz_mod.F90 @ 1939

Last change on this file since 1939 was 186, checked in by lfita, 10 years ago
Removing checking printings
File size: 23.0 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
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_gcm('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
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
99	INCLUDE "iniprint.h"
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_gcm('traclmdz_init', 'pb in allocation 9',1)
136	scavtr(:)=1.
137
138	ALLOCATE( radio(nbtr), stat=ierr)
139	IF (ierr /= 0) CALL abort_gcm('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_gcm('traclmdz_init', 'pb in allocation 2',1)
144
145	ALLOCATE( fshtr(klon,nbtr), stat=ierr)
146	IF (ierr /= 0) CALL abort_gcm('traclmdz_init', 'pb in allocation 3',1)
147
148	ALLOCATE( hsoltr(nbtr), stat=ierr)
149	IF (ierr /= 0) CALL abort_gcm('traclmdz_init', 'pb in allocation 4',1)
150
151	ALLOCATE( tautr(nbtr), stat=ierr)
152	IF (ierr /= 0) CALL abort_gcm('traclmdz_init', 'pb in allocation 5',1)
153	tautr(:) = 0.
154
155	ALLOCATE( vdeptr(nbtr), stat=ierr)
156	IF (ierr /= 0) CALL abort_gcm('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)
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 (id_rn, id_pb, ftsol,pctsrf,masktr,fshtr,hsoltr,tautr,vdeptr,scavtr)
283	END IF
284
285	!
286	! Initialisation de module carbon_cycle_mod
287	! ----------------------------------------------
288	IF (carbon_cycle_tr .OR. carbon_cycle_cpl) THEN
289	CALL carbon_cycle_init(tr_seri, pdtphys, aerosol, radio)
290	END IF
291
292	! Check if all tracers have restart values
293	! ----------------------------------------------
294	DO it=1,nbtr
295	iiq=niadv(it+2)
296	! Test if tracer is zero everywhere.
297	! Done by master process MPI and master thread OpenMP
298	CALL gather(tr_seri(:,:,it),varglo)
299	IF (is_mpi_root .AND. is_omp_root) THEN
300	mintmp=MINVAL(varglo,dim=1)
301	maxtmp=MAXVAL(varglo,dim=1)
302	IF (MINVAL(mintmp,dim=1)==0. .AND. MAXVAL(maxtmp,dim=1)==0.) THEN
303	! Tracer is zero everywhere
304	zero=.TRUE.
305	ELSE
306	zero=.FALSE.
307	END IF
308	END IF
309
310	! Distribute variable at all processes
311	CALL bcast(zero)
312
313	! Initalize tracer that was not found in restart file.
314	IF (zero) THEN
315	! The tracer was not found in restart file or it was equal zero everywhere.
316	WRITE(lunout,*) "The tracer ",trim(tname(iiq))," will be initialized"
317	IF (it==id_pcsat .OR. it==id_pcq .OR. &
318	it==id_pcs0 .OR. it==id_pcq0) THEN
319	tr_seri(:,:,it) = 100.
320	ELSE IF (it==id_pcocsat .OR. it==id_pcos0) THEN
321	DO i = 1, klon
322	IF ( pctsrf (i, is_oce) == 0. ) THEN
323	tr_seri(i,:,it) = 0.
324	ELSE
325	tr_seri(i,:,it) = 100.
326	END IF
327	END DO
328	ELSE
329	! No specific initialization exist for this tracer
330	tr_seri(:,:,it) = 0.
331	END IF
332	END IF
333	END DO
334
335	END SUBROUTINE traclmdz_init
336
337	SUBROUTINE traclmdz(nstep, julien, gmtime, pdtphys, t_seri, paprs, pplay, &
338	cdragh, coefh, yu1, yv1, ftsol, pctsrf, xlat, xlon, couchelimite, sh, &
339	rh, pphi, ustar, zu10m, zv10m, &
340	!! tr_seri, source, solsym, d_tr_cl, zmasse) !RomP
341	tr_seri, source, solsym, d_tr_cl,d_tr_dec, zmasse) !RomP
342
343	USE dimphy
344	USE infotrac
345	USE regr_pr_comb_coefoz_m, ONLY: regr_pr_comb_coefoz
346	USE o3_chem_m, ONLY: o3_chem
347	USE carbon_cycle_mod, ONLY : carbon_cycle, carbon_cycle_tr, carbon_cycle_cpl
348	USE indice_sol_mod
349
350	INCLUDE "YOMCST.h"
351
352	!==========================================================================
353	! -- DESCRIPTION DES ARGUMENTS --
354	!==========================================================================
355
356	! Input arguments
357	!
358	!Configuration grille,temps:
359	INTEGER,INTENT(IN) :: nstep ! nombre d'appels de la physiq
360	INTEGER,INTENT(IN) :: julien ! Jour julien
361	REAL,INTENT(IN) :: gmtime
362	REAL,INTENT(IN) :: pdtphys ! Pas d'integration pour la physique (seconde)
363	REAL,DIMENSION(klon),INTENT(IN) :: xlat ! latitudes pour chaque point
364	REAL, INTENT(IN):: xlon(:) ! dim(klon) longitude
365
366	!
367	!Physique:
368	!--------
369	REAL,DIMENSION(klon,klev),INTENT(IN) :: t_seri ! Temperature
370	REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs ! pression pour chaque inter-couche (en Pa)
371	REAL,DIMENSION(klon,klev),INTENT(IN) :: pplay ! pression pour le mileu de chaque couche (en Pa)
372	REAL,intent(in):: zmasse (:, :) ! dim(klon,klev) density of air, in kg/m2
373
374
375	!Couche limite:
376	!--------------
377	!
378	REAL,DIMENSION(klon),INTENT(IN) :: cdragh ! coeff drag pour T et Q
379	REAL,DIMENSION(klon,klev),INTENT(IN) :: coefh ! coeff melange CL (m**2/s)
380	REAL,DIMENSION(klon),INTENT(IN) :: yu1 ! vents au premier niveau
381	REAL,DIMENSION(klon),INTENT(IN) :: yv1 ! vents au premier niveau
382	LOGICAL,INTENT(IN) :: couchelimite
383	REAL,DIMENSION(klon,klev),INTENT(IN) :: sh ! humidite specifique
384	REAL,DIMENSION(klon,klev),INTENT(IN) :: rh ! Humidite relative
385	REAL,DIMENSION(klon,klev),INTENT(IN) :: pphi ! geopotentie
386	REAL,DIMENSION(klon),INTENT(IN) :: ustar ! ustar (m/s)
387	REAL,DIMENSION(klon),INTENT(IN) :: zu10m ! vent zonal 10m (m/s)
388	REAL,DIMENSION(klon),INTENT(IN) :: zv10m ! vent zonal 10m (m/s)
389
390	! Arguments necessaires pour les sources et puits de traceur:
391	REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: ftsol ! Temperature du sol (surf)(Kelvin)
392	REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: pctsrf ! Pourcentage de sol f(nature du sol)
393
394	! InOutput argument
395	REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT) :: tr_seri ! Concentration Traceur [U/KgA]
396
397	! Output argument
398	CHARACTER(len=8),DIMENSION(nbtr), INTENT(OUT) :: solsym
399	REAL,DIMENSION(klon,nbtr), INTENT(OUT) :: source ! a voir lorsque le flux de surface est prescrit
400	REAL,DIMENSION(klon,klev,nbtr), INTENT(OUT) :: d_tr_cl ! Td couche limite/traceur
401
402	!=======================================================================================
403	! -- VARIABLES LOCALES TRACEURS --
404	!=======================================================================================
405
406	INTEGER :: i, k, it
407	INTEGER :: lmt_pas ! number of time steps of "physics" per day
408
409	REAL,DIMENSION(klon) :: d_trs ! Td dans le reservoir
410	REAL,DIMENSION(klon,klev) :: qsat ! pression de la vapeur a saturation
411	REAL,DIMENSION(klon,klev,nbtr) :: d_tr_dec ! Td radioactive
412	REAL :: zrho ! Masse Volumique de l'air KgA/m3
413	REAL :: amn, amx
414	!
415	!=================================================================
416	! Ajout de la production en Be7 (Beryllium) srcbe U/s/kgA
417	!=================================================================
418	!
419	IF ( id_be /= 0 ) THEN
420	DO k = 1, klev
421	DO i = 1, klon
422	tr_seri(i,k,id_be) = tr_seri(i,k,id_be)+srcbe(i,k)*pdtphys
423	END DO
424	END DO
425	WRITE(,) 'Ajout srcBe dans tr_seri: OK'
426	END IF
427
428
429	!=================================================================
430	! Update pseudo-vapor tracers
431	!=================================================================
432
433	CALL q_sat(klon*klev,t_seri,pplay,qsat)
434
435	IF ( id_pcsat /= 0 ) THEN
436	DO k = 1, klev
437	DO i = 1, klon
438	IF ( pplay(i,k).GE.85000.) THEN
439	tr_seri(i,k,id_pcsat) = qsat(i,k)
440	ELSE
441	tr_seri(i,k,id_pcsat) = MIN (qsat(i,k), tr_seri(i,k,id_pcsat))
442	END IF
443	END DO
444	END DO
445	END IF
446
447	IF ( id_pcocsat /= 0 ) THEN
448	DO k = 1, klev
449	DO i = 1, klon
450	IF ( pplay(i,k).GE.85000.) THEN
451	IF ( pctsrf (i, is_oce) > 0. ) THEN
452	tr_seri(i,k,id_pcocsat) = qsat(i,k)
453	ELSE
454	tr_seri(i,k,id_pcocsat) = 0.
455	END IF
456	ELSE
457	tr_seri(i,k,id_pcocsat) = MIN (qsat(i,k), tr_seri(i,k,id_pcocsat))
458	END IF
459	END DO
460	END DO
461	END IF
462
463	IF ( id_pcq /= 0 ) THEN
464	DO k = 1, klev
465	DO i = 1, klon
466	IF ( pplay(i,k).GE.85000.) THEN
467	tr_seri(i,k,id_pcq) = sh(i,k)
468	ELSE
469	tr_seri(i,k,id_pcq) = MIN (qsat(i,k), tr_seri(i,k,id_pcq))
470	END IF
471	END DO
472	END DO
473	END IF
474
475
476	IF ( id_pcs0 /= 0 ) THEN
477	DO k = 1, klev
478	DO i = 1, klon
479	IF ( pplay(i,k).GE.85000.) THEN
480	tr_seri(i,k,id_pcs0) = qsat(i,k)
481	ELSE
482	tr_seri(i,k,id_pcs0) = MIN (qsat(i,k), tr_seri(i,k,id_pcs0))
483	END IF
484	END DO
485	END DO
486	END IF
487
488
489	IF ( id_pcos0 /= 0 ) THEN
490	DO k = 1, klev
491	DO i = 1, klon
492	IF ( pplay(i,k).GE.85000.) THEN
493	IF ( pctsrf (i, is_oce) > 0. ) THEN
494	tr_seri(i,k,id_pcos0) = qsat(i,k)
495	ELSE
496	tr_seri(i,k,id_pcos0) = 0.
497	END IF
498	ELSE
499	tr_seri(i,k,id_pcos0) = MIN (qsat(i,k), tr_seri(i,k,id_pcos0))
500	END IF
501	END DO
502	END DO
503	END IF
504
505
506	IF ( id_pcq0 /= 0 ) THEN
507	DO k = 1, klev
508	DO i = 1, klon
509	IF ( pplay(i,k).GE.85000.) THEN
510	tr_seri(i,k,id_pcq0) = sh(i,k)
511	ELSE
512	tr_seri(i,k,id_pcq0) = MIN (qsat(i,k), tr_seri(i,k,id_pcq0))
513	END IF
514	END DO
515	END DO
516	END IF
517
518	!=================================================================
519	! Update tracer : Age of stratospheric air
520	!=================================================================
521	IF (id_aga/=0) THEN
522
523	! Bottom layers
524	DO k = 1, lev_1p5km
525	tr_seri(:,k,id_aga) = 0.0
526	END DO
527
528	! Layers above 1.5km
529	DO k = lev_1p5km+1,klev-1
530	tr_seri(:,k,id_aga) = tr_seri(:,k,id_aga) + pdtphys
531	END DO
532
533	! Top layer
534	tr_seri(:,klev,id_aga) = tr_seri(:,klev-1,id_aga)
535
536	END IF
537
538	!======================================================================
539	! -- Calcul de l'effet de la couche limite --
540	!======================================================================
541
542	IF (couchelimite) THEN
543	source(:,:) = 0.0
544
545	IF (id_be /=0) THEN
546	DO i=1, klon
547	zrho = pplay(i,1)/t_seri(i,1)/RD
548	source(i,id_be) = - vdeptr(id_be)tr_seri(i,1,id_be)zrho
549	END DO
550	END IF
551
552	END IF
553
554	DO it=1, nbtr
555	IF (couchelimite .AND. pbl_flg(it) == 0 .AND. (it==id_rn .OR. it==id_pb)) THEN
556	! couche limite avec quantite dans le sol calculee
557	CALL cltracrn(it, id_rn, id_pb, pdtphys, yu1, yv1, &
558	cdragh, coefh,t_seri,ftsol,pctsrf, &
559	tr_seri(:,:,it),trs(:,it), &
560	paprs, pplay, zmasse * rg, &
561	masktr(:,it),fshtr(:,it),hsoltr(it),&
562	tautr(it),vdeptr(it), &
563	xlat,d_tr_cl(:,:,it),d_trs)
564
565	DO k = 1, klev
566	DO i = 1, klon
567	tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cl(i,k,it)
568	END DO
569	END DO
570
571	! Traceur dans le reservoir sol
572	DO i = 1, klon
573	trs(i,it) = trs(i,it) + d_trs(i)
574	END DO
575	END IF
576	END DO
577
578
579	!======================================================================
580	! Calcul de l'effet du puits radioactif
581	!======================================================================
582	CALL radio_decay (id_rn, id_pb, radio,rnpb,pdtphys,tautr,tr_seri,d_tr_dec)
583
584	DO it=1,nbtr
585	WRITE(solsym(it),'(i2)') it
586	END DO
587
588	DO it=1,nbtr
589	IF(radio(it)) then
590	DO k = 1, klev
591	DO i = 1, klon
592	tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_dec(i,k,it)
593	END DO
594	END DO
595	CALL minmaxqfi(tr_seri(:,:,it),0.,1.e33,'puits rn it='//solsym(it))
596	END IF
597	END DO
598
599	!======================================================================
600	! Parameterization of ozone chemistry
601	!======================================================================
602
603	IF (id_o3 /= 0) then
604	lmt_pas = NINT(86400./pdtphys)
605	IF (MOD(nstep - 1, lmt_pas) == 0) THEN
606	! Once per day, update the coefficients for ozone chemistry:
607	CALL regr_pr_comb_coefoz(julien, xlat, paprs, pplay)
608	END IF
609	CALL o3_chem(julien, gmtime, t_seri, zmasse, pdtphys, xlat, &
610	xlon, tr_seri(:, :, id_o3))
611	END IF
612
613	!======================================================================
614	! Calcul de cycle de carbon
615	!======================================================================
616	IF (carbon_cycle_tr .OR. carbon_cycle_cpl) THEN
617	CALL carbon_cycle(nstep, pdtphys, pctsrf, tr_seri, source)
618	END IF
619
620	END SUBROUTINE traclmdz
621
622
623	SUBROUTINE traclmdz_to_restart(trs_out)
624	! This subroutine is called from phyredem.F where the module
625	! variable trs is written to restart file (restartphy.nc)
626	USE dimphy
627	USE infotrac
628
629	REAL,DIMENSION(klon,nbtr), INTENT(OUT) :: trs_out
630	INTEGER :: ierr
631
632	IF ( ALLOCATED(trs) ) THEN
633	trs_out(:,:) = trs(:,:)
634	ELSE
635	! No previous allocate of trs. This is the case for create_etat0_limit.
636	trs_out(:,:) = 0.0
637	END IF
638
639	END SUBROUTINE traclmdz_to_restart
640
641
642	END MODULE traclmdz_mod

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

Download in other formats:

Original Format