Context Navigation

phytrac.F @ 281

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

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

Download in other formats:

Original Format