Context Navigation

phytrac.F @ 343

Last change on this file since 343 was 209, checked in by lmdz, 24 years ago
Manque deux include JQuaas 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 nid_tra
139	SAVE nid_tra
140	c REAL x(klon,klev,nbtr+2) ! traceurs
141	INTEGER ndex(1)
142	REAL zx_tmp_2d(iim,jjm+1), zx_tmp_3d(iim,jjm+1,klev)
143	REAL zx_lon(iim,jjm+1), zx_lat(iim,jjm+1)
144	INTEGER itra
145	SAVE itra ! compteur pour la physique
146	C
147	C Variables liees a l'ecriture de la bande histoire : phytrac.nc
148	c
149	INTEGER ecrit_tra
150	SAVE ecrit_tra
151	logical ok_sync
152	parameter (ok_sync = .true.)
153	C
154	C nature du traceur
155	c
156	logical aerosol(nbtr) ! Nature du traceur
157	c ! aerosol(it) = true => aerosol
158	c ! aerosol(it) = false => gaz
159	c ! nat_trac(it) = 1. aerosolc
160	logical clsol(nbtr) ! clsol(it) = true => CL sol calculee
161	logical radio(nbtr) ! radio(it)=true => decroisssance radioactive
162	save aerosol,clsol,radio
163	C
164	c======================================================================
165	c
166	c Declaration des procedures appelees
167	c
168	c--modif convection tiedtke
169	INTEGER i, k, it,itap
170	save itap
171	REAL delp(klon,klev)
172	c--end modif
173	c
174	c Variables liees a l'ecriture de la bande histoire physique
175	c
176	c Variables locales pour effectuer les appels en serie
177	c----------------------------------------------------
178	c
179	REAL d_tr(klon,klev), d_trs(klon) ! tendances de traceurs
180	REAL d_tr_cl(klon,klev) ! tendance de traceurs couche limite
181	REAL d_tr_cv(klon,klev) ! tendance de traceurs convection
182	REAL d_tr_dec(klon,klev,nbtr) ! la tendance de la decroissance
183	c ! radioactive du rn - > pb
184	REAL d_tr_lessi_impa(klon,klev,nbtr) ! la tendance du lessivage
185	c ! par impaction
186	REAL d_tr_lessi_nucl(klon,klev,nbtr) ! la tendance du lessivage
187	c ! par nucleation
188	REAL fluxrn(klon,klev)
189	REAL fluxpb(klon,klev)
190	REAL pbimpa(klon,klev)
191	REAL pbnucl(klon,klev)
192	REAL rn(klon,klev)
193	REAL pb(klon,klev)
194	REAL flestottr(klon,klev,nbtr) ! flux de lessivage
195	c ! dans chaque couche
196
197	C
198	character*20 modname
199	character*80 abort_message
200	c
201	c Controles
202	c-------------
203	logical first,couchelimite,convection,lessivage,sorties,
204	s rnpb,inirnpb
205	save first,couchelimite,convection,lessivage,sorties,
206	s inirnpb
207	data first,couchelimite,convection,lessivage,sorties
208	s /.true.,.true.,.true.,.true.,.true./
209	c
210	c======================================================================
211	c ecrit_tra = NINT(86400./pdtphys *1.0) ! tous les jours
212	modname='phytrac'
213
214	c print*,'DANS PHYTRAC debutphy=',debutphy
215
216	if (debutphy) then
217
218	print*,'dans phytrac ',pdtphys,ecritphy,ecrit_tra
219	ecrit_tra = NINT(86400./pdtphys/2.) ! tous les 12H
220	c ecrit_tra = NINT(86400./pdtphys) ! tous les 24H
221
222	if(nbtr.lt.nqmax) then
223	print*,'NQMAX=',nqmax
224	print*,'NBTR=',nbtr
225	abort_message='See above'
226	call abort_gcm(modname,abort_message,1)
227	endif
228
229	inirnpb=rnpb
230	PRINT*, 'La frequence de sortie traceurs est ', ecrit_tra
231	itra=0
232	itap=0
233	C
234	CALL ymds2ju(anne_ini, 1, 1, 0.0, zjulian)
235	zjulian = zjulian + day_ini
236	c
237	CALL gr_fi_ecrit(1,klon,iim,jjm+1,xlon,zx_lon)
238	DO i = 1, iim
239	zx_lon(i,1) = xlon(i+1)
240	zx_lon(i,jjm+1) = xlon(i+1)
241	ENDDO
242	c DO ll=1,klev
243	c znivsig(ll)=float(ll)
244	c ENDDO
245	CALL gr_fi_ecrit(1,klon,iim,jjm+1,xlat,zx_lat)
246	CALL histbeg("histrac", iim,zx_lon, jjm+1,zx_lat,
247	. 1,iim,1,jjm+1, 0, zjulian, pdtphys,
248	. nhori, nid_tra)
249	CALL histvert(nid_tra, "presnivs", "Vertical levels", "mb",
250	. klev, presnivs, nvert)
251	zsto = pdtphys
252	zout = pdtphys * FLOAT(ecrit_tra)
253	c
254	CALL histdef(nid_tra, "phis", "Surface geop. height", "-",
255	. iim,jjm+1,nhori, 1,1,1, -99, 32,
256	. "once", zsto,zout)
257	c
258	CALL histdef(nid_tra, "aire", "Grid area", "-",
259	. iim,jjm+1,nhori, 1,1,1, -99, 32,
260	. "once", zsto,zout)
261
262	goto 666
263	CALL histdef(nid_tra, "pyu1", "Vent niv 1", "-",
264	. iim,jjm+1,nhori, 1,1,1, -99, 32,
265	. "inst(X)", zsto,zout)
266
267	CALL histdef(nid_tra, "pyv1", "Vent niv 1", "-",
268	. iim,jjm+1,nhori, 1,1,1, -99, 32,
269	. "inst(X)", zsto,zout)
270	CALL histdef(nid_tra, "psrf1", "nature sol", "-",
271	. iim,jjm+1,nhori, 1,1,1, -99, 32,
272	. "inst(X)", zsto,zout)
273	CALL histdef(nid_tra, "psrf2", "nature sol", "-",
274	. iim,jjm+1,nhori, 1,1,1, -99, 32,
275	. "inst(X)", zsto,zout)
276	CALL histdef(nid_tra, "psrf3", "nature sol", "-",
277	. iim,jjm+1,nhori, 1,1,1, -99, 32,
278	. "inst(X)", zsto,zout)
279	CALL histdef(nid_tra, "psrf4", "nature sol", "-",
280	. iim,jjm+1,nhori, 1,1,1, -99, 32,
281	. "inst(X)", zsto,zout)
282	CALL histdef(nid_tra, "ftsol1", "temper sol", "-",
283	. iim,jjm+1,nhori, 1,1,1, -99, 32,
284	. "inst(X)", zsto,zout)
285	CALL histdef(nid_tra, "ftsol2", "temper sol", "-",
286	. iim,jjm+1,nhori, 1,1,1, -99, 32,
287	. "inst(X)", zsto,zout)
288	CALL histdef(nid_tra, "ftsol3", "temper sol", "-",
289	. iim,jjm+1,nhori, 1,1,1, -99, 32,
290	. "inst", zsto,zout)
291	CALL histdef(nid_tra, "ftsol4", "temper sol", "-",
292	. iim,jjm+1,nhori, 1,1,1, -99, 32,
293	. "inst(X)", zsto,zout)
294	CALL histdef(nid_tra, "pplay", "flux u mont","-",
295	. iim,jjm+1,nhori, klev,1,klev,nvert, 32,
296	. "inst(X)", zsto,zout)
297	CALL histdef(nid_tra, "t", "flux u mont","-",
298	. iim,jjm+1,nhori, klev,1,klev,nvert, 32,
299	. "inst(X)", zsto,zout)
300	CALL histdef(nid_tra, "mfu", "flux u mont","-",
301	. iim,jjm+1,nhori, klev,1,klev,nvert, 32,
302	. "ave(X)", zsto,zout)
303	CALL histdef(nid_tra, "mfd", "flux u decen","-",
304	. iim,jjm+1,nhori, klev,1,klev,nvert, 32,
305	. "ave(X)", zsto,zout)
306	CALL histdef(nid_tra, "en_u", "flux u mont","-",
307	. iim,jjm+1,nhori, klev,1,klev,nvert, 32,
308	. "ave(X)", zsto,zout)
309	CALL histdef(nid_tra, "en_d", "flux u mont","-",
310	. iim,jjm+1,nhori, klev,1,klev,nvert, 32,
311	. "ave(X)", zsto,zout)
312	CALL histdef(nid_tra, "de_u", "flux u mont","-",
313	. iim,jjm+1,nhori, klev,1,klev,nvert, 32,
314	. "ave(X)", zsto,zout)
315	CALL histdef(nid_tra, "de_d", "flux u mont","-",
316	. iim,jjm+1,nhori, klev,1,klev,nvert, 32,
317	. "ave(X)", zsto,zout)
318	CALL histdef(nid_tra, "coefh", "turbulent coef","-",
319	. iim,jjm+1,nhori, klev,1,klev,nvert, 32,
320	. "ave(X)", zsto,zout)
321
322	666 continue
323	c
324	DO it=1,nqmax
325	IF (it.LE.99) THEN
326	WRITE(str2,'(i2.2)') it
327	C champ 3D
328	CALL histdef(nid_tra, "tr"//str2, "Tracer No."//str2, "U/kga",
329	. iim,jjm+1,nhori, klev,1,klev,nvert, 32,
330	. "ave(X)", zsto,zout)
331	IF (lessivage) THEN
332	CALL histdef(nid_tra, "fl"//str2, "Tracer Flux No."//str2,
333	. "U/m2/s",iim,jjm+1,nhori, klev,1,klev,nvert, 32,
334	. "ave(X)", zsto,zout)
335	ENDIF
336	ELSE
337	PRINT*, "Trop de traceurs"
338	CALL abort
339	ENDIF
340	ENDDO
341	CALL histend(nid_tra)
342	ndex(1) = 0
343	c
344	i = NINT(zout/zsto)
345	CALL gr_fi_ecrit(1,klon,iim,jjm+1,pphis,zx_tmp_2d)
346	CALL histwrite(nid_tra,"phis",i,zx_tmp_2d,iim*(jjm+1),ndex)
347	C
348	i = NINT(zout/zsto)
349	CALL gr_fi_ecrit(1,klon,iim,jjm+1,paire,zx_tmp_2d)
350	CALL histwrite(nid_tra,"aire",i,zx_tmp_2d,iim*(jjm+1),ndex)
351
352	c======================================================================
353	c Initialisation de certaines variables pour le Rn et le Pb
354	c======================================================================
355
356	c Initialisation du traceur dans le sol (couche limite radonique)
357	c
358	c print*,'valeur de debut dans phytrac :',debutphy
359	do it=1,nqmax
360	do i=1,klon
361	trs(i,it) = 0.
362	enddo
363	END DO
364
365	open (99,file='starttrac',status='old',
366	. err=999,form='formatted')
367	read(99,*) (trs(i,1),i=1,klon)
368	999 close(99)
369	print*, 'apres starttrac'
370
371	c Initialisation de la fraction d'aerosols lessivee
372	c
373	DO it=1,nqmax
374	DO k = 1, nlev
375	DO i = 1, klon
376	d_tr_lessi_impa(i,k,it) = 0.0
377	d_tr_lessi_nucl(i,k,it) = 0.0
378	ENDDO
379	ENDDO
380	ENDDO
381	c
382	c Initialisation de la nature des traceurs
383	c
384	DO it = 1, nqmax
385	aerosol(it) = .FALSE. ! Tous les traceurs sont des gaz par defaut
386	radio(it) = .FALSE. ! Par defaut pas de passage par radiornpb
387	clsol(it) = .FALSE. ! Par defaut couche limite avec flux prescrit
388	ENDDO
389	c
390	ENDIF ! fin debutphy
391	c Initialisation du traceur dans le sol (couche limite radonique)
392	if(inirnpb) THEN
393	c
394	radio(1)= .true.
395	radio(2)= .true.
396	clsol(1)= .true.
397	clsol(2)= .true.
398	aerosol(2) = .TRUE. ! le Pb est un aerosol
399	c
400	call initrrnpb (ftsol,pctsrf,masktr,fshtr,hsoltr,tautr
401	. ,vdeptr,scavtr)
402	inirnpb=.false.
403	endif
404	if(nqmax.gt.2) aerosol(3)=.true.
405
406
407	c abder
408	goto 777
409	do i=1,nlon
410	pftsol1(i) = ftsol(i,1)
411	pftsol2(i) = ftsol(i,2)
412	pftsol3(i) = ftsol(i,3)
413	pftsol4(i) = ftsol(i,4)
414
415	ppsrf1(i) = pctsrf(i,1)
416	ppsrf2(i) = pctsrf(i,2)
417	ppsrf3(i) = pctsrf(i,3)
418	ppsrf4(i) = pctsrf(i,4)
419
420	enddo
421	ndex(1)=0
422	itap=itap+1
423	CALL gr_fi_ecrit(1,klon,iim,jjm+1,yu1,zx_tmp_2d)
424	CALL histwrite(nid_tra,"pyu1",itap,zx_tmp_2d,
425	s iim*(jjm+1),ndex)
426
427	CALL gr_fi_ecrit(1,klon,iim,jjm+1,yv1,zx_tmp_2d)
428	CALL histwrite(nid_tra,"pyv1",itap,zx_tmp_2d,
429	s iim*(jjm+1),ndex)
430
431	CALL gr_fi_ecrit(1,klon,iim,jjm+1,pftsol1,zx_tmp_2d)
432	CALL histwrite(nid_tra,"ftsol1",itap,zx_tmp_2d,
433	s iim*(jjm+1),ndex)
434
435	CALL gr_fi_ecrit(1,klon,iim,jjm+1,pftsol2,zx_tmp_2d)
436	CALL histwrite(nid_tra,"ftsol2",itap,zx_tmp_2d,
437	s iim*(jjm+1),ndex)
438
439	CALL gr_fi_ecrit(1,klon,iim,jjm+1,pftsol3,zx_tmp_2d)
440	CALL histwrite(nid_tra,"ftsol3",itap,zx_tmp_2d,
441	s iim*(jjm+1),ndex)
442
443	CALL gr_fi_ecrit(1,klon,iim,jjm+1,pftsol4,zx_tmp_2d)
444	CALL histwrite(nid_tra,"ftsol4",itap,zx_tmp_2d,
445	s iim*(jjm+1),ndex)
446
447	CALL gr_fi_ecrit(1,klon,iim,jjm+1,ppsrf1,zx_tmp_2d)
448	CALL histwrite(nid_tra,"psrf1",itap,zx_tmp_2d,
449	s iim*(jjm+1),ndex)
450
451	CALL gr_fi_ecrit(1,klon,iim,jjm+1,ppsrf2,zx_tmp_2d)
452	CALL histwrite(nid_tra,"psrf2",itap,zx_tmp_2d,
453	s iim*(jjm+1),ndex)
454
455	CALL gr_fi_ecrit(1,klon,iim,jjm+1,ppsrf3,zx_tmp_2d)
456	CALL histwrite(nid_tra,"psrf3",itap,zx_tmp_2d,
457	s iim*(jjm+1),ndex)
458
459	CALL gr_fi_ecrit(1,klon,iim,jjm+1,ppsrf4,zx_tmp_2d)
460	CALL histwrite(nid_tra,"psrf4",itap,zx_tmp_2d,
461	s iim*(jjm+1),ndex)
462	777 continue
463	c======================================================================
464	c Calcul de l'effet de la convection
465	c======================================================================
466	print*,'Avant convection'
467	do it=1,nqmax
468	WRITE(itn,'(i1)') it
469	c call diagtracphy(tr_seri(:,:,it),paprs,'Avant conv'//itn)
470	enddo
471
472	if (convection) then
473
474	print*,'Pas de temps dans phytrac : ',pdtphys
475	DO it=1, nqmax
476	CALL nflxtr(pdtphys, pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
477	. pplay, paprs, tr_seri(1,1,it), d_tr_cv)
478	DO k = 1, nlev
479	DO i = 1, klon
480	tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cv(i,k)
481	ENDDO
482	ENDDO
483	c WRITE(itn,'(i1)') it
484	c CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'convection it='//itn)
485	ENDDO
486	c print*,'apres nflxtr'
487
488
489	endif ! convection
490	c print*,'Apres convection'
491	c do it=1,nqmax
492	c WRITE(itn,'(i1)') it
493	c call diagtracphy(tr_seri(:,:,it),paprs,'Avant conv'//itn)
494	c enddo
495
496	c======================================================================
497	c Calcul de l'effet de la couche limite
498	c======================================================================
499	c print *,'Avant couchelimite'
500	c do it=1,nqmax
501	c WRITE(itn,'(i1)') it
502	c call diagtracphy(tr_seri(:,:,it),paprs,'Avant CL '//itn)
503	c enddo
504
505	if (couchelimite) then
506
507	DO k = 1, nlev
508	DO i = 1, klon
509	delp(i,k) = paprs(i,k)-paprs(i,k+1)
510	ENDDO
511	ENDDO
512
513	C maf modif pour tenir compte du cas rnpb + traceur
514	DO it=1, nqmax
515	c print *,'it',it,clsol(it)
516	if (clsol(it)) then ! couche limite avec quantite dans le sol calculee
517	CALL cltracrn(it, pdtphys, yu1, yv1,
518	e coefh,t_seri,ftsol,pctsrf,
519	e tr_seri(1,1,it),trs(1,it),
520	e paprs, pplay, delp,
521	e masktr(1,it),fshtr(1,it),hsoltr(it),
522	e tautr(it),vdeptr(it),
523	e xlat,
524	s d_tr_cl,d_trs)
525	DO k = 1, nlev
526	DO i = 1, klon
527	tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cl(i,k)
528	ENDDO
529	ENDDO
530	c
531	c Traceur ds sol
532	c
533	DO i = 1, klon
534	trs(i,it) = trs(i,it) + d_trs(i)
535	END DO
536	C
537	C maf provisoire suppression des prints
538	C WRITE(itn,'(i1)') it
539	C CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'cltracrn it='//itn)
540	else ! couche limite avec flux prescrit
541	Cmaf provisoire source / traceur a creer
542	DO i=1, klon
543	source(i) = 0.0 ! pas de source, pour l'instant
544	ENDDO
545	C
546	CALL cltrac(pdtphys, coefh,t_seri,
547	s tr_seri(1,1,it), source,
548	e paprs, pplay, delp,
549	s d_tr )
550	DO k = 1, nlev
551	DO i = 1, klon
552	tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr(i,k)
553	ENDDO
554	ENDDO
555	Cmaf provisoire suppression des prints
556	Cmaf WRITE(itn,'(i1)') it
557	cmaf CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'cltracn it='//itn)
558	endif
559	ENDDO
560	c
561	endif ! couche limite
562
563	c print*,'Apres couchelimite'
564	c do it=1,nqmax
565	c WRITE(itn,'(i1)') it
566	c call diagtracphy(tr_seri(:,:,it),paprs,'Avant CL '//itn)
567	c enddo
568
569	c======================================================================
570	c Calcul de l'effet du puits radioactif
571	c======================================================================
572
573	C MAF il faudrait faire une modification pour passer dans radiornpb
574	c si radio=true mais pour l'instant radiornpb propre au cas rnpb
575	if(rnpb) then
576	call radiornpb (tr_seri,pdtphys,tautr,d_tr_dec)
577	C
578	DO it=1,nqmax
579	if(radio(it)) then
580	DO k = 1, nlev
581	DO i = 1, klon
582	tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_dec(i,k,it)
583	ENDDO
584	ENDDO
585	WRITE(itn,'(i1)') it
586	CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'puits rn it='//itn)
587	endif
588	ENDDO
589	c
590	endif ! rnpb decroissance radioactive
591	C
592	c======================================================================
593	c Calcul de l'effet de la precipitation
594	c======================================================================
595
596	print*,'LESSIVAGE =',lessivage
597	IF (lessivage) THEN
598
599	c print*,'avant lessivage'
600
601	DO it = 1, nqmax
602	DO k = 1, nlev
603	DO i = 1, klon
604	d_tr_lessi_nucl(i,k,it) = 0.0
605	d_tr_lessi_impa(i,k,it) = 0.0
606	flestottr(i,k,it) = 0.0
607	ENDDO
608	ENDDO
609	ENDDO
610	c
611	c tendance des aerosols nuclees et impactes
612	c
613	DO it = 1, nqmax
614	IF (aerosol(it)) THEN
615	DO k = 1, nlev
616	DO i = 1, klon
617	d_tr_lessi_nucl(i,k,it) = d_tr_lessi_nucl(i,k,it) +
618	s ( 1 - frac_nucl(i,k) )*tr_seri(i,k,it)
619	d_tr_lessi_impa(i,k,it) = d_tr_lessi_impa(i,k,it) +
620	s ( 1 - frac_impa(i,k) )*tr_seri(i,k,it)
621	ENDDO
622	ENDDO
623	ENDIF
624	ENDDO
625	c
626	c Mises a jour des traceurs + calcul des flux de lessivage
627	c Mise a jour due a l'impaction et a la nucleation
628	c
629	c call dump2d(iim,jjm-1,frac_impa(2:klon-1,10),'FRACIMPA')
630	c call dump2d(iim,jjm-1,frac_nucl(2:klon-1,10),'FRACNUCL')
631	c call dump2d(iim,jjm-1,tr_seri(2:klon-1,10,3),'TRACEUR3')
632	DO it = 1, nqmax
633	c print*,'IT=',it,aerosol(it)
634	IF (aerosol(it)) THEN
635	c print*,'IT=',it,' On lessive'
636	DO k = 1, nlev
637	DO i = 1, klon
638	tr_seri(i,k,it)=tr_seri(i,k,it)
639	s frac_impa(i,k)frac_nucl(i,k)
640	ENDDO
641	ENDDO
642	ENDIF
643	ENDDO
644	c call dump2d(iim,jjm-1,tr_seri(2:klon-1,10,3),'TRACEUR3B')
645	c
646	c Flux lessivage total
647	c
648	DO it = 1, nqmax
649	DO k = 1, nlev
650	DO i = 1, klon
651	flestottr(i,k,it) = flestottr(i,k,it) -
652	s ( d_tr_lessi_nucl(i,k,it) +
653	s d_tr_lessi_impa(i,k,it) ) *
654	s ( paprs(i,k)-paprs(i,k+1) ) /
655	s (RG * pdtphys)
656	ENDDO
657	ENDDO
658	c
659	Cmaf suppression provisoire
660	Cmaf WRITE(itn,'(i1)') it
661	Cmaf CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'tr(lessi) it='//itn)
662	ENDDO
663	c
664	c print*,'apres lessivage'
665	ENDIF
666	Cc
667	DO k = 1, nlev
668	DO i = 1, klon
669	fluxrn(i,k) = flestottr(i,k,1)
670	fluxpb(i,k) = flestottr(i,k,2)
671	rn(i,k) = tr_seri(i,k,1)
672	pb(i,k) = tr_seri(i,k,2)
673	pbnucl(i,k)=d_tr_lessi_nucl(i,k,2)
674	pbimpa(i,k)=d_tr_lessi_impa(i,k,2)
675	ENDDO
676	ENDDO
677	itra=itra+1
678	ndex(1) = 0
679	DO it=1,nqmax
680	IF (it.LE.99) THEN
681	WRITE(str2,'(i2.2)') it
682	C champs 3D
683	CALL gr_fi_ecrit(klev,klon,iim,jjm+1,tr_seri(1,1,it),zx_tmp_3d)
684	CALL histwrite(nid_tra,"tr"//str2,itra,zx_tmp_3d,
685	. iim(jjm+1)klev,ndex)
686	c IF (lessivage) THEN
687	c CALL gr_fi_ecrit(klev,klon,iim,jjm+1,flestottr(1,1,it),zx_tmp_3d)
688	c CALL histwrite(nid_tra,"fl"//str2,itra,zx_tmp_3d,
689	c . iim(jjm+1)klev,ndex)
690	c ENDIF
691	ELSE
692	PRINT*, "Trop de traceurs"
693	CALL abort
694	ENDIF
695	ENDDO
696
697	goto 888
698	CALL gr_fi_ecrit(klev,klon,iim,jjm+1,pplay,zx_tmp_3d)
699	CALL histwrite(nid_tra,"pplay",itra,zx_tmp_3d,
700	. iim(jjm+1)klev,ndex)
701
702	CALL gr_fi_ecrit(klev,klon,iim,jjm+1,t_seri,zx_tmp_3d)
703	CALL histwrite(nid_tra,"t",itra,zx_tmp_3d,
704	. iim(jjm+1)klev,ndex)
705	CALL gr_fi_ecrit(klev,klon,iim,jjm+1,pmfu,zx_tmp_3d)
706	CALL histwrite(nid_tra,"mfu",itra,zx_tmp_3d,
707	. iim(jjm+1)klev,ndex)
708	CALL gr_fi_ecrit(klev,klon,iim,jjm+1,pmfd,zx_tmp_3d)
709	CALL histwrite(nid_tra,"mfd",itra,zx_tmp_3d,
710	. iim(jjm+1)klev,ndex)
711	CALL gr_fi_ecrit(klev,klon,iim,jjm+1,pen_u,zx_tmp_3d)
712	CALL histwrite(nid_tra,"en_u",itra,zx_tmp_3d,
713	. iim(jjm+1)klev,ndex)
714	CALL gr_fi_ecrit(klev,klon,iim,jjm+1,pen_d,zx_tmp_3d)
715	CALL histwrite(nid_tra,"en_d",itra,zx_tmp_3d,
716	. iim(jjm+1)klev,ndex)
717	CALL gr_fi_ecrit(klev,klon,iim,jjm+1,pde_d,zx_tmp_3d)
718	CALL histwrite(nid_tra,"de_d",itra,zx_tmp_3d,
719	. iim(jjm+1)klev,ndex)
720	CALL gr_fi_ecrit(klev,klon,iim,jjm+1,pde_u,zx_tmp_3d)
721	CALL histwrite(nid_tra,"de_u",itra,zx_tmp_3d,
722	. iim(jjm+1)klev,ndex)
723	CALL gr_fi_ecrit(klev,klon,iim,jjm+1,coefh,zx_tmp_3d)
724	CALL histwrite(nid_tra,"coefh",itra,zx_tmp_3d,
725	. iim(jjm+1)klev,ndex)
726
727	888 continue
728
729	c print*,'Sortie phytrac'
730	c do it=1,nqmax
731	c WRITE(itn,'(i1)') it
732	c call diagtracphy(tr_seri(:,:,it),paprs,'Fin Phys '//itn)
733	c enddo
734
735	if (lafin) then
736	print*, 'c est la fin de la physique'
737	open (99,file='restarttrac', form='formatted')
738	do i=1,klon
739	write(99,*) trs(i,1)
740	enddo
741	PRINT*, 'Ecriture du fichier restarttrac'
742	close(99)
743	else
744	print*, 'physique pas fini'
745	endif
746
747
748	RETURN
749	END

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

Download in other formats:

Original Format