Context Navigation

phytrac.F @ 5357

Last change on this file since 5357 was 349, checked in by lmdz, 23 years ago
Regle le probleme de decalage de 1 jour au debut de chaque simulation LF
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 26.0 KB

Line
1	c
2	c $Header$
3	c
4	SUBROUTINE phytrac (rnpb,
5	I debutphy,lafin,
6	I nqmax,
7	I nlon,nlev,pdtphys,
8	I t_seri,paprs,pplay,
9	I pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
10	I coefh,yu1,yv1,ftsol,pctsrf,xlat,
11	I frac_impa, frac_nucl,
12	I xlon,presnivs,paire,pphis,
13	O tr_seri)
14	USE ioipsl
15	IMPLICIT none
16	c======================================================================
17	c Auteur(s) FH
18	c Objet: Moniteur general des tendances traceurs
19	c
20	cAA Remarques en vrac:
21	cAA--------------------
22	cAA 1/ le call phytrac se fait avec nqmax-2 donc nous avons bien
23	cAA les vrais traceurs (nbtr) dans phytrac (pas la vapeur ni eau liquide)
24	cAA 2/ Le choix du radon et du pb se fait juste avec un data
25	cAA (peu propre). Peut-etre pourrait-on prevoir dans l'avenir
26	cAA une variable "type de traceur"
27	c======================================================================
28	#include "YOMCST.h"
29	#include "dimensions.h"
30	#include "dimphy.h"
31	#include "indicesol.h"
32	#include "temps.h"
33	#include "control.h"
34	c======================================================================
35
36	c Arguments:
37	c
38	c EN ENTREE:
39	c ==========
40	c
41	c divers:
42	c -------
43	c
44	integer nlon ! nombre de points horizontaux
45	integer nlev ! nombre de couches verticales
46	integer nqmax ! nombre de traceurs auxquels on applique la physique
47	real pdtphys ! pas d'integration pour la physique (seconde)
48	real t_seri(nlon,nlev) ! temperature
49	real tr_seri(nlon,nlev,nbtr) ! traceur
50	real paprs(nlon,nlev+1) ! pression pour chaque inter-couche (en Pa)
51	real pplay(nlon,nlev) ! pression pour le mileu de chaque couche (en Pa)
52	real presnivs(klev) ! pressions approximat. des milieux couches ( en PA)
53	real paire(klon)
54	real pphis(klon)
55	logical debutphy ! le flag de l'initialisation de la physique
56	logical lafin ! le flag de la fin de la physique
57
58	integer ll
59	c
60	cAA Rem : nbtr : nombre de vrais traceurs est defini dans dimphy.h
61	c
62	c convection:
63	c -----------
64	c
65	REAL pmfu(nlon,nlev) ! flux de masse dans le panache montant
66	REAL pmfd(nlon,nlev) ! flux de masse dans le panache descendant
67	REAL pen_u(nlon,nlev) ! flux entraine dans le panache montant
68	REAL pde_u(nlon,nlev) ! flux detraine dans le panache montant
69	REAL pen_d(nlon,nlev) ! flux entraine dans le panache descendant
70	REAL pde_d(nlon,nlev) ! flux detraine dans le panache descendant
71	c
72	c Couche limite:
73	c --------------
74	c
75	REAL coefh(nlon,nlev) ! coeff melange CL
76	REAL yu1(nlon) ! vents au premier niveau
77	REAL yv1(nlon) ! vents au premier niveau
78	REAL xlat(nlon) ! latitudes pour chaque point
79	REAL xlon(nlon) ! longitudes pour chaque point
80
81	c
82	c Lessivage:
83	c ----------
84	c
85	c pour le ON-LINE
86	c
87	REAL frac_impa(nlon,nlev) ! fraction d'aerosols impactes
88	REAL frac_nucl(nlon,nlev) ! fraction d'aerosols nuclees
89	c
90	cAA
91	cAA Arguments necessaires pour les sources et puits de traceur:
92	cAA ----------------
93	cAA
94	real ftsol(nlon,nbsrf) ! Temperature du sol (surf)(Kelvin)
95	real pctsrf(nlon,nbsrf) ! Pourcentage de sol f(nature du sol)
96	c abder
97	real pftsol1(nlon),pftsol2(nlon),pftsol3(nlon),pftsol4(nlon)
98	real ppsrf1(nlon),ppsrf2(nlon),ppsrf3(nlon),ppsrf4(nlon)
99	c fin
100	cAA ----------------------------
101	cAA VARIABLES LOCALES TRACEURS
102	cAA ----------------------------
103	cAA
104	cAA Sources et puits des traceurs:
105	cAA ------------------------------
106	cAA
107	cAA Pour l'instant seuls les cas du rn et du pb ont ete envisages.
108
109	REAL source(klon) ! a voir lorsque le flux est prescrit
110	REAL topuits(klev,nbtr) ! a voir
111	cAA
112	cAA Pour la source de radon et son reservoir de sol
113	cAA ................................................
114
115	REAL trs(klon,nbtr) ! Conc. radon ds le sol
116	SAVE trs
117
118	REAL masktr(klon,nbtr) ! Masque reservoir de sol traceur
119	c Masque de l'echange avec la surface
120	c (1 = reservoir) ou (possible => 1 )
121	SAVE masktr
122	REAL fshtr(klon,nbtr) ! Flux surfacique dans le reservoir de sol
123	SAVE fshtr
124	REAL hsoltr(nbtr) ! Epaisseur equivalente du reservoir de sol
125	SAVE hsoltr
126	REAL tautr(nbtr) ! Constante de decroissance radioactive
127	SAVE tautr
128	REAL vdeptr(nbtr) ! Vitesse de depot sec dans la couche Brownienne
129	SAVE vdeptr
130	REAL scavtr(nbtr) ! Coefficient de lessivage
131	SAVE scavtr
132	cAA
133	CHARACTER*1 itn
134	C maf ioipsl
135	CHARACTER*2 str2
136	INTEGER nhori, nvert
137	REAL zsto, zout, zjulian
138	integer idayref
139	INTEGER nid_tra
140	SAVE nid_tra
141	c REAL x(klon,klev,nbtr+2) ! traceurs
142	INTEGER ndex(1)
143	REAL zx_tmp_2d(iim,jjm+1), zx_tmp_3d(iim,jjm+1,klev)
144	REAL zx_lon(iim,jjm+1), zx_lat(iim,jjm+1)
145	INTEGER itra
146	SAVE itra ! compteur pour la physique
147	C
148	C Variables liees a l'ecriture de la bande histoire : phytrac.nc
149	c
150	INTEGER ecrit_tra
151	SAVE ecrit_tra
152	logical ok_sync
153	parameter (ok_sync = .true.)
154	C
155	C nature du traceur
156	c
157	logical aerosol(nbtr) ! Nature du traceur
158	c ! aerosol(it) = true => aerosol
159	c ! aerosol(it) = false => gaz
160	c ! nat_trac(it) = 1. aerosolc
161	logical clsol(nbtr) ! clsol(it) = true => CL sol calculee
162	logical radio(nbtr) ! radio(it)=true => decroisssance radioactive
163	save aerosol,clsol,radio
164	C
165	c======================================================================
166	c
167	c Declaration des procedures appelees
168	c
169	c--modif convection tiedtke
170	INTEGER i, k, it,itap
171	save itap
172	REAL delp(klon,klev)
173	c--end modif
174	c
175	c Variables liees a l'ecriture de la bande histoire physique
176	c
177	c Variables locales pour effectuer les appels en serie
178	c----------------------------------------------------
179	c
180	REAL d_tr(klon,klev), d_trs(klon) ! tendances de traceurs
181	REAL d_tr_cl(klon,klev) ! tendance de traceurs couche limite
182	REAL d_tr_cv(klon,klev) ! tendance de traceurs convection
183	REAL d_tr_dec(klon,klev,nbtr) ! la tendance de la decroissance
184	c ! radioactive du rn - > pb
185	REAL d_tr_lessi_impa(klon,klev,nbtr) ! la tendance du lessivage
186	c ! par impaction
187	REAL d_tr_lessi_nucl(klon,klev,nbtr) ! la tendance du lessivage
188	c ! par nucleation
189	REAL fluxrn(klon,klev)
190	REAL fluxpb(klon,klev)
191	REAL pbimpa(klon,klev)
192	REAL pbnucl(klon,klev)
193	REAL rn(klon,klev)
194	REAL pb(klon,klev)
195	REAL flestottr(klon,klev,nbtr) ! flux de lessivage
196	c ! dans chaque couche
197
198	C
199	character*20 modname
200	character*80 abort_message
201	c
202	c Controles
203	c-------------
204	logical first,couchelimite,convection,lessivage,sorties,
205	s rnpb,inirnpb
206	save first,couchelimite,convection,lessivage,sorties,
207	s inirnpb
208	data first,couchelimite,convection,lessivage,sorties
209	s /.true.,.true.,.true.,.true.,.true./
210	c
211	c======================================================================
212	c ecrit_tra = NINT(86400./pdtphys *1.0) ! tous les jours
213	modname='phytrac'
214
215	c print*,'DANS PHYTRAC debutphy=',debutphy
216
217	if (debutphy) then
218
219	print*,'dans phytrac ',pdtphys,ecritphy,ecrit_tra
220	ecrit_tra = NINT(86400./pdtphys/2.) ! tous les 12H
221	c ecrit_tra = NINT(86400./pdtphys) ! tous les 24H
222
223	if(nbtr.lt.nqmax) then
224	print*,'NQMAX=',nqmax
225	print*,'NBTR=',nbtr
226	abort_message='See above'
227	call abort_gcm(modname,abort_message,1)
228	endif
229
230	inirnpb=rnpb
231	PRINT*, 'La frequence de sortie traceurs est ', ecrit_tra
232	itra=0
233	itap=0
234	C
235	idayref = day_ini
236	CALL ymds2ju(anne_ini, 1, idayref, 0.0, zjulian)
237	c
238	CALL gr_fi_ecrit(1,klon,iim,jjm+1,xlon,zx_lon)
239	DO i = 1, iim
240	zx_lon(i,1) = xlon(i+1)
241	zx_lon(i,jjm+1) = xlon(i+1)
242	ENDDO
243	c DO ll=1,klev
244	c znivsig(ll)=float(ll)
245	c ENDDO
246	CALL gr_fi_ecrit(1,klon,iim,jjm+1,xlat,zx_lat)
247	CALL histbeg("histrac", iim,zx_lon, jjm+1,zx_lat,
248	. 1,iim,1,jjm+1, 0, zjulian, pdtphys,
249	. nhori, nid_tra)
250	CALL histvert(nid_tra, "presnivs", "Vertical levels", "mb",
251	. klev, presnivs, nvert)
252	zsto = pdtphys
253	zout = pdtphys * FLOAT(ecrit_tra)
254	c
255	CALL histdef(nid_tra, "phis", "Surface geop. height", "-",
256	. iim,jjm+1,nhori, 1,1,1, -99, 32,
257	. "once", zsto,zout)
258	c
259	CALL histdef(nid_tra, "aire", "Grid area", "-",
260	. iim,jjm+1,nhori, 1,1,1, -99, 32,
261	. "once", zsto,zout)
262
263	goto 666
264	CALL histdef(nid_tra, "pyu1", "Vent niv 1", "-",
265	. iim,jjm+1,nhori, 1,1,1, -99, 32,
266	. "inst(X)", zsto,zout)
267
268	CALL histdef(nid_tra, "pyv1", "Vent niv 1", "-",
269	. iim,jjm+1,nhori, 1,1,1, -99, 32,
270	. "inst(X)", zsto,zout)
271	CALL histdef(nid_tra, "psrf1", "nature sol", "-",
272	. iim,jjm+1,nhori, 1,1,1, -99, 32,
273	. "inst(X)", zsto,zout)
274	CALL histdef(nid_tra, "psrf2", "nature sol", "-",
275	. iim,jjm+1,nhori, 1,1,1, -99, 32,
276	. "inst(X)", zsto,zout)
277	CALL histdef(nid_tra, "psrf3", "nature sol", "-",
278	. iim,jjm+1,nhori, 1,1,1, -99, 32,
279	. "inst(X)", zsto,zout)
280	CALL histdef(nid_tra, "psrf4", "nature sol", "-",
281	. iim,jjm+1,nhori, 1,1,1, -99, 32,
282	. "inst(X)", zsto,zout)
283	CALL histdef(nid_tra, "ftsol1", "temper sol", "-",
284	. iim,jjm+1,nhori, 1,1,1, -99, 32,
285	. "inst(X)", zsto,zout)
286	CALL histdef(nid_tra, "ftsol2", "temper sol", "-",
287	. iim,jjm+1,nhori, 1,1,1, -99, 32,
288	. "inst(X)", zsto,zout)
289	CALL histdef(nid_tra, "ftsol3", "temper sol", "-",
290	. iim,jjm+1,nhori, 1,1,1, -99, 32,
291	. "inst", zsto,zout)
292	CALL histdef(nid_tra, "ftsol4", "temper sol", "-",
293	. iim,jjm+1,nhori, 1,1,1, -99, 32,
294	. "inst(X)", zsto,zout)
295	CALL histdef(nid_tra, "pplay", "flux u mont","-",
296	. iim,jjm+1,nhori, klev,1,klev,nvert, 32,
297	. "inst(X)", zsto,zout)
298	CALL histdef(nid_tra, "t", "flux u mont","-",
299	. iim,jjm+1,nhori, klev,1,klev,nvert, 32,
300	. "inst(X)", zsto,zout)
301	CALL histdef(nid_tra, "mfu", "flux u mont","-",
302	. iim,jjm+1,nhori, klev,1,klev,nvert, 32,
303	. "ave(X)", zsto,zout)
304	CALL histdef(nid_tra, "mfd", "flux u decen","-",
305	. iim,jjm+1,nhori, klev,1,klev,nvert, 32,
306	. "ave(X)", zsto,zout)
307	CALL histdef(nid_tra, "en_u", "flux u mont","-",
308	. iim,jjm+1,nhori, klev,1,klev,nvert, 32,
309	. "ave(X)", zsto,zout)
310	CALL histdef(nid_tra, "en_d", "flux u mont","-",
311	. iim,jjm+1,nhori, klev,1,klev,nvert, 32,
312	. "ave(X)", zsto,zout)
313	CALL histdef(nid_tra, "de_u", "flux u mont","-",
314	. iim,jjm+1,nhori, klev,1,klev,nvert, 32,
315	. "ave(X)", zsto,zout)
316	CALL histdef(nid_tra, "de_d", "flux u mont","-",
317	. iim,jjm+1,nhori, klev,1,klev,nvert, 32,
318	. "ave(X)", zsto,zout)
319	CALL histdef(nid_tra, "coefh", "turbulent coef","-",
320	. iim,jjm+1,nhori, klev,1,klev,nvert, 32,
321	. "ave(X)", zsto,zout)
322
323	666 continue
324	c
325	DO it=1,nqmax
326	IF (it.LE.99) THEN
327	WRITE(str2,'(i2.2)') it
328	C champ 3D
329	CALL histdef(nid_tra, "tr"//str2, "Tracer No."//str2, "U/kga",
330	. iim,jjm+1,nhori, klev,1,klev,nvert, 32,
331	. "ave(X)", zsto,zout)
332	IF (lessivage) THEN
333	CALL histdef(nid_tra, "fl"//str2, "Tracer Flux No."//str2,
334	. "U/m2/s",iim,jjm+1,nhori, klev,1,klev,nvert, 32,
335	. "ave(X)", zsto,zout)
336	ENDIF
337	ELSE
338	PRINT*, "Trop de traceurs"
339	CALL abort
340	ENDIF
341	ENDDO
342	CALL histend(nid_tra)
343	ndex(1) = 0
344	c
345	i = NINT(zout/zsto)
346	CALL gr_fi_ecrit(1,klon,iim,jjm+1,pphis,zx_tmp_2d)
347	CALL histwrite(nid_tra,"phis",i,zx_tmp_2d,iim*(jjm+1),ndex)
348	C
349	i = NINT(zout/zsto)
350	CALL gr_fi_ecrit(1,klon,iim,jjm+1,paire,zx_tmp_2d)
351	CALL histwrite(nid_tra,"aire",i,zx_tmp_2d,iim*(jjm+1),ndex)
352
353	c======================================================================
354	c Initialisation de certaines variables pour le Rn et le Pb
355	c======================================================================
356
357	c Initialisation du traceur dans le sol (couche limite radonique)
358	c
359	c print*,'valeur de debut dans phytrac :',debutphy
360	do it=1,nqmax
361	do i=1,klon
362	trs(i,it) = 0.
363	enddo
364	END DO
365
366	open (99,file='starttrac',status='old',
367	. err=999,form='formatted')
368	read(99,*) (trs(i,1),i=1,klon)
369	999 close(99)
370	print*, 'apres starttrac'
371
372	c Initialisation de la fraction d'aerosols lessivee
373	c
374	DO it=1,nqmax
375	DO k = 1, nlev
376	DO i = 1, klon
377	d_tr_lessi_impa(i,k,it) = 0.0
378	d_tr_lessi_nucl(i,k,it) = 0.0
379	ENDDO
380	ENDDO
381	ENDDO
382	c
383	c Initialisation de la nature des traceurs
384	c
385	DO it = 1, nqmax
386	aerosol(it) = .FALSE. ! Tous les traceurs sont des gaz par defaut
387	radio(it) = .FALSE. ! Par defaut pas de passage par radiornpb
388	clsol(it) = .FALSE. ! Par defaut couche limite avec flux prescrit
389	ENDDO
390	c
391	ENDIF ! fin debutphy
392	c Initialisation du traceur dans le sol (couche limite radonique)
393	if(inirnpb) THEN
394	c
395	radio(1)= .true.
396	radio(2)= .true.
397	clsol(1)= .true.
398	clsol(2)= .true.
399	aerosol(2) = .TRUE. ! le Pb est un aerosol
400	c
401	call initrrnpb (ftsol,pctsrf,masktr,fshtr,hsoltr,tautr
402	. ,vdeptr,scavtr)
403	inirnpb=.false.
404	endif
405	if(nqmax.gt.2) aerosol(3)=.true.
406
407
408	c abder
409	goto 777
410	do i=1,nlon
411	pftsol1(i) = ftsol(i,1)
412	pftsol2(i) = ftsol(i,2)
413	pftsol3(i) = ftsol(i,3)
414	pftsol4(i) = ftsol(i,4)
415
416	ppsrf1(i) = pctsrf(i,1)
417	ppsrf2(i) = pctsrf(i,2)
418	ppsrf3(i) = pctsrf(i,3)
419	ppsrf4(i) = pctsrf(i,4)
420
421	enddo
422	ndex(1)=0
423	itap=itap+1
424	CALL gr_fi_ecrit(1,klon,iim,jjm+1,yu1,zx_tmp_2d)
425	CALL histwrite(nid_tra,"pyu1",itap,zx_tmp_2d,
426	s iim*(jjm+1),ndex)
427
428	CALL gr_fi_ecrit(1,klon,iim,jjm+1,yv1,zx_tmp_2d)
429	CALL histwrite(nid_tra,"pyv1",itap,zx_tmp_2d,
430	s iim*(jjm+1),ndex)
431
432	CALL gr_fi_ecrit(1,klon,iim,jjm+1,pftsol1,zx_tmp_2d)
433	CALL histwrite(nid_tra,"ftsol1",itap,zx_tmp_2d,
434	s iim*(jjm+1),ndex)
435
436	CALL gr_fi_ecrit(1,klon,iim,jjm+1,pftsol2,zx_tmp_2d)
437	CALL histwrite(nid_tra,"ftsol2",itap,zx_tmp_2d,
438	s iim*(jjm+1),ndex)
439
440	CALL gr_fi_ecrit(1,klon,iim,jjm+1,pftsol3,zx_tmp_2d)
441	CALL histwrite(nid_tra,"ftsol3",itap,zx_tmp_2d,
442	s iim*(jjm+1),ndex)
443
444	CALL gr_fi_ecrit(1,klon,iim,jjm+1,pftsol4,zx_tmp_2d)
445	CALL histwrite(nid_tra,"ftsol4",itap,zx_tmp_2d,
446	s iim*(jjm+1),ndex)
447
448	CALL gr_fi_ecrit(1,klon,iim,jjm+1,ppsrf1,zx_tmp_2d)
449	CALL histwrite(nid_tra,"psrf1",itap,zx_tmp_2d,
450	s iim*(jjm+1),ndex)
451
452	CALL gr_fi_ecrit(1,klon,iim,jjm+1,ppsrf2,zx_tmp_2d)
453	CALL histwrite(nid_tra,"psrf2",itap,zx_tmp_2d,
454	s iim*(jjm+1),ndex)
455
456	CALL gr_fi_ecrit(1,klon,iim,jjm+1,ppsrf3,zx_tmp_2d)
457	CALL histwrite(nid_tra,"psrf3",itap,zx_tmp_2d,
458	s iim*(jjm+1),ndex)
459
460	CALL gr_fi_ecrit(1,klon,iim,jjm+1,ppsrf4,zx_tmp_2d)
461	CALL histwrite(nid_tra,"psrf4",itap,zx_tmp_2d,
462	s iim*(jjm+1),ndex)
463	777 continue
464	c======================================================================
465	c Calcul de l'effet de la convection
466	c======================================================================
467	print*,'Avant convection'
468	do it=1,nqmax
469	WRITE(itn,'(i1)') it
470	c call diagtracphy(tr_seri(:,:,it),paprs,'Avant conv'//itn)
471	enddo
472
473	if (convection) then
474
475	print*,'Pas de temps dans phytrac : ',pdtphys
476	DO it=1, nqmax
477	CALL nflxtr(pdtphys, pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
478	. pplay, paprs, tr_seri(1,1,it), d_tr_cv)
479	DO k = 1, nlev
480	DO i = 1, klon
481	tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cv(i,k)
482	ENDDO
483	ENDDO
484	c WRITE(itn,'(i1)') it
485	c CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'convection it='//itn)
486	ENDDO
487	c print*,'apres nflxtr'
488
489
490	endif ! convection
491	c print*,'Apres convection'
492	c do it=1,nqmax
493	c WRITE(itn,'(i1)') it
494	c call diagtracphy(tr_seri(:,:,it),paprs,'Avant conv'//itn)
495	c enddo
496
497	c======================================================================
498	c Calcul de l'effet de la couche limite
499	c======================================================================
500	c print *,'Avant couchelimite'
501	c do it=1,nqmax
502	c WRITE(itn,'(i1)') it
503	c call diagtracphy(tr_seri(:,:,it),paprs,'Avant CL '//itn)
504	c enddo
505
506	if (couchelimite) then
507
508	DO k = 1, nlev
509	DO i = 1, klon
510	delp(i,k) = paprs(i,k)-paprs(i,k+1)
511	ENDDO
512	ENDDO
513
514	C maf modif pour tenir compte du cas rnpb + traceur
515	DO it=1, nqmax
516	c print *,'it',it,clsol(it)
517	if (clsol(it)) then ! couche limite avec quantite dans le sol calculee
518	CALL cltracrn(it, pdtphys, yu1, yv1,
519	e coefh,t_seri,ftsol,pctsrf,
520	e tr_seri(1,1,it),trs(1,it),
521	e paprs, pplay, delp,
522	e masktr(1,it),fshtr(1,it),hsoltr(it),
523	e tautr(it),vdeptr(it),
524	e xlat,
525	s d_tr_cl,d_trs)
526	DO k = 1, nlev
527	DO i = 1, klon
528	tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cl(i,k)
529	ENDDO
530	ENDDO
531	c
532	c Traceur ds sol
533	c
534	DO i = 1, klon
535	trs(i,it) = trs(i,it) + d_trs(i)
536	END DO
537	C
538	C maf provisoire suppression des prints
539	C WRITE(itn,'(i1)') it
540	C CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'cltracrn it='//itn)
541	else ! couche limite avec flux prescrit
542	Cmaf provisoire source / traceur a creer
543	DO i=1, klon
544	source(i) = 0.0 ! pas de source, pour l'instant
545	ENDDO
546	C
547	CALL cltrac(pdtphys, coefh,t_seri,
548	s tr_seri(1,1,it), source,
549	e paprs, pplay, delp,
550	s d_tr )
551	DO k = 1, nlev
552	DO i = 1, klon
553	tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr(i,k)
554	ENDDO
555	ENDDO
556	Cmaf provisoire suppression des prints
557	Cmaf WRITE(itn,'(i1)') it
558	cmaf CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'cltracn it='//itn)
559	endif
560	ENDDO
561	c
562	endif ! couche limite
563
564	c print*,'Apres couchelimite'
565	c do it=1,nqmax
566	c WRITE(itn,'(i1)') it
567	c call diagtracphy(tr_seri(:,:,it),paprs,'Avant CL '//itn)
568	c enddo
569
570	c======================================================================
571	c Calcul de l'effet du puits radioactif
572	c======================================================================
573
574	C MAF il faudrait faire une modification pour passer dans radiornpb
575	c si radio=true mais pour l'instant radiornpb propre au cas rnpb
576	if(rnpb) then
577	call radiornpb (tr_seri,pdtphys,tautr,d_tr_dec)
578	C
579	DO it=1,nqmax
580	if(radio(it)) then
581	DO k = 1, nlev
582	DO i = 1, klon
583	tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_dec(i,k,it)
584	ENDDO
585	ENDDO
586	WRITE(itn,'(i1)') it
587	CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'puits rn it='//itn)
588	endif
589	ENDDO
590	c
591	endif ! rnpb decroissance radioactive
592	C
593	c======================================================================
594	c Calcul de l'effet de la precipitation
595	c======================================================================
596
597	print*,'LESSIVAGE =',lessivage
598	IF (lessivage) THEN
599
600	c print*,'avant lessivage'
601
602	DO it = 1, nqmax
603	DO k = 1, nlev
604	DO i = 1, klon
605	d_tr_lessi_nucl(i,k,it) = 0.0
606	d_tr_lessi_impa(i,k,it) = 0.0
607	flestottr(i,k,it) = 0.0
608	ENDDO
609	ENDDO
610	ENDDO
611	c
612	c tendance des aerosols nuclees et impactes
613	c
614	DO it = 1, nqmax
615	IF (aerosol(it)) THEN
616	DO k = 1, nlev
617	DO i = 1, klon
618	d_tr_lessi_nucl(i,k,it) = d_tr_lessi_nucl(i,k,it) +
619	s ( 1 - frac_nucl(i,k) )*tr_seri(i,k,it)
620	d_tr_lessi_impa(i,k,it) = d_tr_lessi_impa(i,k,it) +
621	s ( 1 - frac_impa(i,k) )*tr_seri(i,k,it)
622	ENDDO
623	ENDDO
624	ENDIF
625	ENDDO
626	c
627	c Mises a jour des traceurs + calcul des flux de lessivage
628	c Mise a jour due a l'impaction et a la nucleation
629	c
630	c call dump2d(iim,jjm-1,frac_impa(2:klon-1,10),'FRACIMPA')
631	c call dump2d(iim,jjm-1,frac_nucl(2:klon-1,10),'FRACNUCL')
632	c call dump2d(iim,jjm-1,tr_seri(2:klon-1,10,3),'TRACEUR3')
633	DO it = 1, nqmax
634	c print*,'IT=',it,aerosol(it)
635	IF (aerosol(it)) THEN
636	c print*,'IT=',it,' On lessive'
637	DO k = 1, nlev
638	DO i = 1, klon
639	tr_seri(i,k,it)=tr_seri(i,k,it)
640	s frac_impa(i,k)frac_nucl(i,k)
641	ENDDO
642	ENDDO
643	ENDIF
644	ENDDO
645	c call dump2d(iim,jjm-1,tr_seri(2:klon-1,10,3),'TRACEUR3B')
646	c
647	c Flux lessivage total
648	c
649	DO it = 1, nqmax
650	DO k = 1, nlev
651	DO i = 1, klon
652	flestottr(i,k,it) = flestottr(i,k,it) -
653	s ( d_tr_lessi_nucl(i,k,it) +
654	s d_tr_lessi_impa(i,k,it) ) *
655	s ( paprs(i,k)-paprs(i,k+1) ) /
656	s (RG * pdtphys)
657	ENDDO
658	ENDDO
659	c
660	Cmaf suppression provisoire
661	Cmaf WRITE(itn,'(i1)') it
662	Cmaf CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'tr(lessi) it='//itn)
663	ENDDO
664	c
665	c print*,'apres lessivage'
666	ENDIF
667	Cc
668	DO k = 1, nlev
669	DO i = 1, klon
670	fluxrn(i,k) = flestottr(i,k,1)
671	fluxpb(i,k) = flestottr(i,k,2)
672	rn(i,k) = tr_seri(i,k,1)
673	pb(i,k) = tr_seri(i,k,2)
674	pbnucl(i,k)=d_tr_lessi_nucl(i,k,2)
675	pbimpa(i,k)=d_tr_lessi_impa(i,k,2)
676	ENDDO
677	ENDDO
678	itra=itra+1
679	ndex(1) = 0
680	DO it=1,nqmax
681	IF (it.LE.99) THEN
682	WRITE(str2,'(i2.2)') it
683	C champs 3D
684	CALL gr_fi_ecrit(klev,klon,iim,jjm+1,tr_seri(1,1,it),zx_tmp_3d)
685	CALL histwrite(nid_tra,"tr"//str2,itra,zx_tmp_3d,
686	. iim(jjm+1)klev,ndex)
687	c IF (lessivage) THEN
688	c CALL gr_fi_ecrit(klev,klon,iim,jjm+1,flestottr(1,1,it),zx_tmp_3d)
689	c CALL histwrite(nid_tra,"fl"//str2,itra,zx_tmp_3d,
690	c . iim(jjm+1)klev,ndex)
691	c ENDIF
692	ELSE
693	PRINT*, "Trop de traceurs"
694	CALL abort
695	ENDIF
696	ENDDO
697
698	goto 888
699	CALL gr_fi_ecrit(klev,klon,iim,jjm+1,pplay,zx_tmp_3d)
700	CALL histwrite(nid_tra,"pplay",itra,zx_tmp_3d,
701	. iim(jjm+1)klev,ndex)
702
703	CALL gr_fi_ecrit(klev,klon,iim,jjm+1,t_seri,zx_tmp_3d)
704	CALL histwrite(nid_tra,"t",itra,zx_tmp_3d,
705	. iim(jjm+1)klev,ndex)
706	CALL gr_fi_ecrit(klev,klon,iim,jjm+1,pmfu,zx_tmp_3d)
707	CALL histwrite(nid_tra,"mfu",itra,zx_tmp_3d,
708	. iim(jjm+1)klev,ndex)
709	CALL gr_fi_ecrit(klev,klon,iim,jjm+1,pmfd,zx_tmp_3d)
710	CALL histwrite(nid_tra,"mfd",itra,zx_tmp_3d,
711	. iim(jjm+1)klev,ndex)
712	CALL gr_fi_ecrit(klev,klon,iim,jjm+1,pen_u,zx_tmp_3d)
713	CALL histwrite(nid_tra,"en_u",itra,zx_tmp_3d,
714	. iim(jjm+1)klev,ndex)
715	CALL gr_fi_ecrit(klev,klon,iim,jjm+1,pen_d,zx_tmp_3d)
716	CALL histwrite(nid_tra,"en_d",itra,zx_tmp_3d,
717	. iim(jjm+1)klev,ndex)
718	CALL gr_fi_ecrit(klev,klon,iim,jjm+1,pde_d,zx_tmp_3d)
719	CALL histwrite(nid_tra,"de_d",itra,zx_tmp_3d,
720	. iim(jjm+1)klev,ndex)
721	CALL gr_fi_ecrit(klev,klon,iim,jjm+1,pde_u,zx_tmp_3d)
722	CALL histwrite(nid_tra,"de_u",itra,zx_tmp_3d,
723	. iim(jjm+1)klev,ndex)
724	CALL gr_fi_ecrit(klev,klon,iim,jjm+1,coefh,zx_tmp_3d)
725	CALL histwrite(nid_tra,"coefh",itra,zx_tmp_3d,
726	. iim(jjm+1)klev,ndex)
727
728	888 continue
729
730	c print*,'Sortie phytrac'
731	c do it=1,nqmax
732	c WRITE(itn,'(i1)') it
733	c call diagtracphy(tr_seri(:,:,it),paprs,'Fin Phys '//itn)
734	c enddo
735
736	if (lafin) then
737	print*, 'c est la fin de la physique'
738	open (99,file='restarttrac', form='formatted')
739	do i=1,klon
740	write(99,*) trs(i,1)
741	enddo
742	PRINT*, 'Ecriture du fichier restarttrac'
743	close(99)
744	else
745	print*, 'physique pas fini'
746	endif
747
748
749	RETURN
750	END

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

Download in other formats:

Original Format