Context Navigation

phytrac.F @ 673

Last change on this file since 673 was 673, checked in by lmdzadmin, 19 years ago
suppression de prt_flag_ts (INCA) AC MAF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 26.4 KB

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

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

Download in other formats:

Original Format