Context Navigation

phytrac.F @ 3536

Last change on this file since 3536 was 524, checked in by lmdzadmin, 21 years ago
Initial revision
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 23.3 KB

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

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

Download in other formats:

Original Format