Context Navigation

phytrac.F @ 180

Last change on this file since 180 was 177, checked in by lmdzadmin, 23 years ago

Lots of stuff, plus particulierement:

appel a ORCHIDEE en etat de marche (pb de grille subsiste)
modifs de Pascale sur soil dans le cas ou ok_veget=false -

Property svn:eol-style set to native
Property svn:executable set to *
Property svn:keywords set to Author Date Id Revision

File size: 17.6 KB

Line
1	SUBROUTINE phytrac (rnpb,
2	I debutphy,
3	I nqmax,
4	I nlon,nlev,pdtphys,
5	I t_seri,paprs,pplay,
6	I pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
7	I coefh,yu1,yv1,ftsol,pctsrf,xlat,
8	I frac_impa, frac_nucl,
9	I xlon,presnivs,paire,pphis,
10	O tr_seri)
11	USE ioipsl
12	USE histcom
13
14	IMPLICIT none
15	c======================================================================
16	c Auteur(s) FH
17	c Objet: Moniteur general des tendances traceurs
18	c
19	cAA Remarques en vrac:
20	cAA--------------------
21	cAA 1/ le call phytrac se fait avec nqmax-2 donc nous avons bien
22	cAA les vrais traceurs (nbtr) dans phytrac (pas la vapeur ni eau liquide)
23	cAA 2/ Le choix du radon et du pb se fait juste avec un data
24	cAA (peu propre). Peut-etre pourrait-on prevoir dans l'avenir
25	cAA une variable "type de traceur"
26	c======================================================================
27	#include "YOMCST.h"
28	#include "dimensions.h"
29	#include "dimphy.h"
30	#include "indicesol.h"
31	#include "control.h"
32	#include "temps.h"
33	c======================================================================
34
35	c Arguments:
36	c
37	c EN ENTREE:
38	c ==========
39	c
40	c divers:
41	c -------
42	c
43	integer nlon ! nombre de points horizontaux
44	integer nlev ! nombre de couches verticales
45	integer nqmax ! nombre de traceurs auxquels on applique la physique
46	real pdtphys ! pas d'integration pour la physique (seconde)
47	real t_seri(nlon,nlev) ! temperature
48	real tr_seri(nlon,nlev,nbtr) ! traceur
49	real paprs(nlon,nlev+1) ! pression pour chaque inter-couche (en Pa)
50	real pplay(nlon,nlev) ! pression pour le mileu de chaque couche (en Pa)
51	real presnivs(klev) ! pressions approximat. des milieux couches ( en PA)
52	real znivsig(klev) ! niveaux sigma
53	real paire(klon)
54	real pphis(klon)
55	logical debutphy ! le flag de l'initialisation de la physique
56	integer ll
57	c
58	cAA Rem : nbtr : nombre de vrais traceurs est defini dans dimphy.h
59	c
60	c convection:
61	c -----------
62	c
63	REAL pmfu(nlon,nlev) ! flux de masse dans le panache montant
64	REAL pmfd(nlon,nlev) ! flux de masse dans le panache descendant
65	REAL pen_u(nlon,nlev) ! flux entraine dans le panache montant
66	REAL pde_u(nlon,nlev) ! flux detraine dans le panache montant
67	REAL pen_d(nlon,nlev) ! flux entraine dans le panache descendant
68	REAL pde_d(nlon,nlev) ! flux detraine dans le panache descendant
69	c
70	c Couche limite:
71	c --------------
72	c
73	REAL coefh(nlon,nlev) ! coeff melange CL
74	REAL yu1(nlon) ! vents au premier niveau
75	REAL yv1(nlon) ! vents au premier niveau
76	REAL xlat(nlon) ! latitudes pour chaque point
77	REAL xlon(nlon) ! longitudes pour chaque point
78
79	c
80	c Lessivage:
81	c ----------
82	c
83	c pour le ON-LINE
84	c
85	REAL frac_impa(nlon,nlev) ! fraction d'aerosols impactes
86	REAL frac_nucl(nlon,nlev) ! fraction d'aerosols nuclees
87	c
88	cAA
89	cAA Arguments necessaires pour les sources et puits de traceur:
90	cAA ----------------
91	cAA
92	real ftsol(nlon,nbsrf) ! Temperature du sol (surf)(Kelvin)
93	real pctsrf(nlon,nbsrf) ! Pourcentage de sol f(nature du sol)
94
95	cAA ----------------------------
96	cAA VARIABLES LOCALES TRACEURS
97	cAA ----------------------------
98	cAA
99	cAA Sources et puits des traceurs:
100	cAA ------------------------------
101	cAA
102	cAA Pour l'instant seuls les cas du rn et du pb ont ete envisages.
103
104	REAL source(klon) ! a voir lorsque le flux est prescrit
105	REAL topuits(klev,nbtr) ! a voir
106	cAA
107	cAA Pour la source de radon et son reservoir de sol
108	cAA ................................................
109
110	REAL trs(klon,nbtr) ! Conc. radon ds le sol
111	SAVE trs
112
113	REAL masktr(klon,nbtr) ! Masque reservoir de sol traceur
114	c Masque de l'echange avec la surface
115	c (1 = reservoir) ou (possible => 1 )
116	SAVE masktr
117	REAL fshtr(klon,nbtr) ! Flux surfacique dans le reservoir de sol
118	SAVE fshtr
119	REAL hsoltr(nbtr) ! Epaisseur equivalente du reservoir de sol
120	SAVE hsoltr
121	REAL tautr(nbtr) ! Constante de decroissance radioactive
122	SAVE tautr
123	REAL vdeptr(nbtr) ! Vitesse de depot sec dans la couche Brownienne
124	SAVE vdeptr
125	REAL scavtr(nbtr) ! Coefficient de lessivage
126	SAVE scavtr
127	cAA
128	CHARACTER*1 itn
129	C maf ioipsl
130	CHARACTER*2 str2
131	INTEGER nhori, nvert
132	REAL zsto, zout, zjulian
133	INTEGER nid_tra
134	SAVE nid_tra
135	INTEGER ndex2d(iim(jjm+1)),ndex3d(iim(jjm+1)*klev)
136	REAL zx_tmp_2d(iim,jjm+1), zx_tmp_3d(iim,jjm+1,klev)
137	REAL zx_lon(iim,jjm+1), zx_lat(iim,jjm+1)
138	INTEGER itra
139	SAVE itra ! compteur pour la physique
140	C
141	C Variables liees a l'ecriture de la bande histoire : phytrac.nc
142	c
143	INTEGER ecrit_tra
144	SAVE ecrit_tra
145	logical ok_sync
146	parameter (ok_sync = .true.)
147	C
148	C nature du traceur
149	c
150	logical aerosol(nbtr) ! Nature du traceur
151	c ! aerosol(it) = true => aerosol
152	c ! aerosol(it) = false => gaz
153	c ! nat_trac(it) = 1. aerosolc
154	logical clsol(nbtr) ! clsol(it) = true => CL sol calculee
155	logical radio(nbtr) ! radio(it)=true => decroisssance radioactive
156	save aerosol,clsol,radio
157	C
158	c======================================================================
159	c
160	c Declaration des procedures appelees
161	c
162	c--modif convection tiedtke
163	INTEGER i, k, it
164
165	REAL delp(klon,klev)
166	c--end modif
167	c
168	c Variables liees a l'ecriture de la bande histoire physique
169	c
170	c Variables locales pour effectuer les appels en serie
171	c----------------------------------------------------
172	c
173	REAL d_tr(klon,klev), d_trs(klon) ! tendances de traceurs
174	REAL d_tr_cl(klon,klev) ! tendance de traceurs couche limite
175	REAL d_tr_cv(klon,klev) ! tendance de traceurs convection
176	REAL d_tr_dec(klon,klev,nbtr) ! la tendance de la decroissance
177	c ! radioactive du rn - > pb
178	REAL d_tr_lessi_impa(klon,klev,nbtr) ! la tendance du lessivage
179	c ! par impaction
180	REAL d_tr_lessi_nucl(klon,klev,nbtr) ! la tendance du lessivage
181	c ! par nucleation
182	REAL fluxrn(klon,klev)
183	REAL fluxpb(klon,klev)
184	REAL pbimpa(klon,klev)
185	REAL pbnucl(klon,klev)
186	REAL rn(klon,klev)
187	REAL pb(klon,klev)
188	REAL flestottr(klon,klev,nbtr) ! flux de lessivage
189	c ! dans chaque couche
190
191	C
192	character*20 modname
193	character*80 abort_message
194	c
195	c Controles
196	c-------------
197	logical first,couchelimite,convection,lessivage,sorties,
198	s rnpb,inirnpb
199	save first,couchelimite,convection,lessivage,sorties,
200	s inirnpb
201	data first,couchelimite,convection,lessivage,sorties
202	s /.true.,.true.,.true.,.true.,.true./
203	c
204	c======================================================================
205	c ecrit_tra = NINT(86400./pdtphys *1.0) ! tous les jours
206	modname='phytrac'
207
208	c print*,'DANS PHYTRAC debutphy=',debutphy
209
210	ecrit_tra = NINT(86400./pdtphys *ecritphy)
211	zsto = pdtphys
212	zout = pdtphys * FLOAT(ecrit_tra)
213	if (debutphy) then
214
215	if(nbtr.lt.nqmax) then
216	print*,'NQMAX=',nqmax
217	print*,'NBTR=',nbtr
218	abort_message='See above'
219	call abort_gcm(modname,abort_message,1)
220	endif
221
222	inirnpb=rnpb
223	PRINT*, 'La frequence de sortie traceurs est ', ecrit_tra
224	itra=0
225	C
226	CALL ymds2ju(anne_ini, 1, 1, 0.0, zjulian)
227	zjulian = zjulian + day_ini
228	c
229	CALL gr_fi_ecrit(1,klon,iim,jjm+1,xlon,zx_lon)
230	DO i = 1, iim
231	zx_lon(i,1) = xlon(i+1)
232	zx_lon(i,jjm+1) = xlon(i+1)
233	ENDDO
234	DO ll=1,klev
235	znivsig(ll)=float(ll)
236	ENDDO
237	CALL gr_fi_ecrit(1,klon,iim,jjm+1,xlat,zx_lat)
238	CALL histbeg("histrac", iim,zx_lon, jjm+1,zx_lat,
239	. 1,iim,1,jjm+1, 0, zjulian, pdtphys,
240	. nhori, nid_tra)
241	call histvert(nid_tra, 'sig_s', 'Niveaux sigma','-',
242	. klev, znivsig, nvert)
243	C
244	C CALL histvert(nid_tra, "presnivs", "Vertical levels", "mb",
245	C . klev, presnivs, nvert)
246	c
247	CALL histdef(nid_tra, "phis", "Surface geop. height", "-",
248	. iim,jjm+1,nhori, 1,1,1, -99, 32,
249	. "once", zsto,zout)
250	c
251	CALL histdef(nid_tra, "aire", "Grid area", "-",
252	. iim,jjm+1,nhori, 1,1,1, -99, 32,
253	. "once", zsto,zout)
254	c
255	DO it=1,nqmax
256	IF (it.LE.99) THEN
257	WRITE(str2,'(i2.2)') it
258	C champ 3D
259	CALL histdef(nid_tra, "tr"//str2, "Tracer No."//str2, "U/kga",
260	. iim,jjm+1,nhori, klev,1,klev,nvert, 32,
261	. "ave(X)", zsto,zout)
262	IF (lessivage) THEN
263	CALL histdef(nid_tra, "fl"//str2, "Tracer Flux No."//str2,
264	. "U/m2/s",iim,jjm+1,nhori, klev,1,klev,nvert, 32,
265	. "ave(X)", zsto,zout)
266	ENDIF
267	ELSE
268	PRINT*, "Trop de traceurs"
269	CALL abort
270	ENDIF
271	ENDDO
272	CALL histend(nid_tra)
273
274	c======================================================================
275	c Initialisation de certaines variables pour le Rn et le Pb
276	c======================================================================
277
278	c Initialisation du traceur dans le sol (couche limite radonique)
279	c
280	c print*,'valeur de debut dans phytrac :',debutphy
281	do it=1,nqmax
282	do i=1,klon
283	trs(i,it) = 0.
284	enddo
285	END DO
286	c Initialisation de la fraction d'aerosols lessivee
287	c
288	DO it=1,nqmax
289	DO k = 1, nlev
290	DO i = 1, klon
291	d_tr_lessi_impa(i,k,it) = 0.0
292	d_tr_lessi_nucl(i,k,it) = 0.0
293	ENDDO
294	ENDDO
295	ENDDO
296	c
297	c Initialisation de la nature des traceurs
298	c
299	DO it = 1, nqmax
300	aerosol(it) = .FALSE. ! Tous les traceurs sont des gaz par defaut
301	radio(it) = .FALSE. ! Par defaut pas de passage par radiornpb
302	clsol(it) = .FALSE. ! Par defaut couche limite avec flux prescrit
303	ENDDO
304	c
305	ENDIF ! fin debutphy
306	c Initialisation du traceur dans le sol (couche limite radonique)
307	if(inirnpb) THEN
308	c
309	radio(1)= .true.
310	radio(2)= .true.
311	clsol(1)= .true.
312	clsol(2)= .true.
313	aerosol(2) = .TRUE. ! le Pb est un aerosol
314	c
315	call initrrnpb (ftsol,pctsrf,masktr,fshtr,hsoltr,tautr
316	. ,vdeptr,scavtr)
317	inirnpb=.false.
318	endif
319	c======================================================================
320	c Calcul de l'effet de la convection
321	c======================================================================
322
323	if (convection) then
324
325	c print*,'Pas de temps dans phytrac : ',pdtphys
326	DO it=1, nqmax
327	CALL nflxtr(pdtphys, pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
328	. pplay, paprs, tr_seri(1,1,it), d_tr_cv)
329	DO k = 1, nlev
330	DO i = 1, klon
331	tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cv(i,k)
332	ENDDO
333	ENDDO
334	WRITE(itn,'(i1)') it
335	CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'convection it='//itn)
336	ENDDO
337	c print*,'apres nflxtr'
338
339
340	endif ! convection
341
342	c======================================================================
343	c Calcul de l'effet de la couche limite
344	c======================================================================
345
346	c print*,'avant couchelimite'
347	if (couchelimite) then
348
349	DO k = 1, nlev
350	DO i = 1, klon
351	delp(i,k) = paprs(i,k)-paprs(i,k+1)
352	ENDDO
353	ENDDO
354
355	C maf modif pour tenir compte du cas rnpb + traceur
356	DO it=1, nqmax
357	c print *,'it',it,clsol(it)
358	if (clsol(it)) then ! couche limite avec quantite dans le sol calculee
359	CALL cltracrn(it, pdtphys, yu1, yv1,
360	e coefh,t_seri,ftsol,pctsrf,
361	e tr_seri(1,1,it),trs(1,it),
362	e paprs, pplay, delp,
363	e masktr(1,it),fshtr(1,it),hsoltr(it),
364	e tautr(it),vdeptr(it),
365	e xlat,
366	s d_tr_cl,d_trs)
367	DO k = 1, nlev
368	DO i = 1, klon
369	tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cl(i,k)
370	ENDDO
371	ENDDO
372	c
373	c Traceur ds sol
374	c
375	DO i = 1, klon
376	trs(i,it) = trs(i,it) + d_trs(i)
377	END DO
378	C
379	C maf provisoire suppression des prints
380	C WRITE(itn,'(i1)') it
381	C CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'cltracrn it='//itn)
382	else ! couche limite avec flux prescrit
383	Cmaf provisoire source / traceur a creer
384	DO i=1, klon
385	source(i) = 0.0 ! pas de source, pour l'instant
386	ENDDO
387	C
388	CALL cltrac(pdtphys, coefh,t_seri,
389	s tr_seri(1,1,it), source,
390	e paprs, pplay, delp,
391	s d_tr )
392	DO k = 1, nlev
393	DO i = 1, klon
394	tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr(i,k)
395	ENDDO
396	ENDDO
397	Cmaf provisoire suppression des prints
398	Cmaf WRITE(itn,'(i1)') it
399	cmaf CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'cltracn it='//itn)
400	endif
401	ENDDO
402	c
403	endif ! couche limite
404
405	c print*,'apres couchelimite'
406
407	c======================================================================
408	c Calcul de l'effet du puits radioactif
409	c======================================================================
410
411	C MAF il faudrait faire une modification pour passer dans radiornpb
412	c si radio=true mais pour l'instant radiornpb propre au cas rnpb
413	if(rnpb) then
414	call radiornpb (tr_seri,pdtphys,tautr,d_tr_dec)
415	C
416	DO it=1,nqmax
417	if(radio(it)) then
418	DO k = 1, nlev
419	DO i = 1, klon
420	tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_dec(i,k,it)
421	ENDDO
422	ENDDO
423	WRITE(itn,'(i1)') it
424	CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'puits rn it='//itn)
425	endif
426	ENDDO
427	c
428	endif ! rnpb decroissance radioactive
429	C
430	c======================================================================
431	c Calcul de l'effet de la precipitation
432	c======================================================================
433
434	IF (lessivage) THEN
435
436	c print*,'avant lessivage'
437
438	DO it = 1, nqmax
439	DO k = 1, nlev
440	DO i = 1, klon
441	d_tr_lessi_nucl(i,k,it) = 0.0
442	d_tr_lessi_impa(i,k,it) = 0.0
443	flestottr(i,k,it) = 0.0
444	ENDDO
445	ENDDO
446	ENDDO
447	c
448	c tendance des aerosols nuclees et impactes
449	c
450	DO it = 1, nqmax
451	IF (aerosol(it)) THEN
452	DO k = 1, nlev
453	DO i = 1, klon
454	d_tr_lessi_nucl(i,k,it) = d_tr_lessi_nucl(i,k,it) +
455	s ( 1 - frac_nucl(i,k) )*tr_seri(i,k,it)
456	d_tr_lessi_impa(i,k,it) = d_tr_lessi_impa(i,k,it) +
457	s ( 1 - frac_impa(i,k) )*tr_seri(i,k,it)
458	ENDDO
459	ENDDO
460	ENDIF
461	ENDDO
462	c
463	c Mises a jour des traceurs + calcul des flux de lessivage
464	c Mise a jour due a l'impaction et a la nucleation
465	c
466	DO it = 1, nqmax
467	IF (aerosol(it)) THEN
468	DO k = 1, nlev
469	DO i = 1, klon
470	tr_seri(i,k,it) = tr_seri(i,k,it) *
471	s ( frac_impa(i,k) + frac_nucl(i,k) - 1. )
472	ENDDO
473	ENDDO
474	ENDIF
475	ENDDO
476	c
477	c Flux lessivage total
478	c
479	DO it = 1, nqmax
480	DO k = 1, nlev
481	DO i = 1, klon
482	flestottr(i,k,it) = flestottr(i,k,it) -
483	s ( d_tr_lessi_nucl(i,k,it) +
484	s d_tr_lessi_impa(i,k,it) ) *
485	s ( paprs(i,k)-paprs(i,k+1) ) /
486	s (RG * pdtphys)
487	ENDDO
488	ENDDO
489	c
490	Cmaf suppression provisoire
491	Cmaf WRITE(itn,'(i1)') it
492	Cmaf CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'tr(lessi) it='//itn)
493	ENDDO
494	c
495	c print*,'apres lessivage'
496	ENDIF
497	Cc
498	DO k = 1, nlev
499	DO i = 1, klon
500	fluxrn(i,k) = flestottr(i,k,1)
501	fluxpb(i,k) = flestottr(i,k,2)
502	rn(i,k) = tr_seri(i,k,1)
503	pb(i,k) = tr_seri(i,k,2)
504	pbnucl(i,k)=d_tr_lessi_nucl(i,k,2)
505	pbimpa(i,k)=d_tr_lessi_impa(i,k,2)
506	ENDDO
507	ENDDO
508	itra=itra+1
509
510	C
511	C Sorties IOIPSL
512	ndex2d = 0
513	ndex3d = 0
514	c
515	c write(,)'sorties ioipsl phytrac',zsto,zout
516	CALL gr_fi_ecrit(1,klon,iim,jjm+1,pphis,zx_tmp_2d)
517	CALL histwrite(nid_tra,"phis",itra,zx_tmp_2d,iim*(jjm+1),ndex2d)
518	C
519	CALL gr_fi_ecrit(1,klon,iim,jjm+1,paire,zx_tmp_2d)
520	CALL histwrite(nid_tra,"aire",itra,zx_tmp_2d,iim*(jjm+1),ndex2d)
521	DO it=1,nqmax
522	IF (it.LE.99) THEN
523	WRITE(str2,'(i2.2)') it
524	C champs 3D
525	CALL gr_fi_ecrit(klev,klon,iim,jjm+1,tr_seri(1,1,it),zx_tmp_3d)
526	CALL histwrite(nid_tra,"tr"//str2,itra,zx_tmp_3d,
527	. iim(jjm+1)klev,ndex3d)
528	IF (lessivage) THEN
529	CALL gr_fi_ecrit(klev,klon,iim,jjm+1,flestottr(1,1,it),zx_tmp_3d)
530	CALL histwrite(nid_tra,"fl"//str2,itra,zx_tmp_3d,
531	. iim(jjm+1)klev,ndex3d)
532	ENDIF
533	ELSE
534	PRINT*, "Trop de traceurs"
535	CALL abort
536	ENDIF
537	ENDDO
538	if (ok_sync) call histsync(nid_tra)
539
540	RETURN
541	END

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

Download in other formats:

Original Format