Context Navigation

phytrac.F @ 338

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

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

Download in other formats:

Original Format