Context Navigation

phytrac.F @ 207

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

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

Download in other formats:

Original Format