Context Navigation

phytrac.F @ 679

Last change on this file since 679 was 679, checked in by lmdzadmin, 18 years ago
Rajout variable prt_flag_ts pour limiter sorties INCA - AC MAF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 28.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 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	c 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	INTEGER :: prt_flag_ts(nbtr)=(/
349	#ifdef INCA_CH4
350	. 1,1,1,0,0,1,1,1,1,1,
351	. 0,1,0,0,0,0,0,1,0,0,
352	. 0,1,1,1,1,0,1,1,1,0,
353	. 1,1,1,1,1,1,1,1,1,1,
354	. 1,0,0
355	#ifdef INCA_AER
356	. ,1,1,1,1,0,1,1,1,1,0,
357	. 1,1,1,1,1,1,0,1,0,1,
358	. 1,1,1,1,0,1,0,1,1,1
359	#endif
360	#endif
361	#ifdef INCA_NMHC
362	. 1,1,1,1,1,1,1,1,1,1,
363	. 1,1,1,1,1,1,1,1,1,1,
364	. 1,1,1,1,1,1,1,1,1,1,
365	. 1,1,1,1,1,1,1,1,1,1,
366	. 1,1,1,1,1,1,1,1,1,1,
367	. 1,1,1,1,1,1,1,1,1,1,
368	. 1,1,1,1,1,1,1,1,1,1,
369	. 1,1,1,1,1,1,1,1,1,1,
370	. 1,1,1,1,1,1,1
371	#ifdef INCA_AER
372	. ,1,1,1,1,0,1,1,1,1,0,
373	. 1,1,1,1,1,1,0,1,0,1,
374	. 1,1,1,1,0,1,0,1,1,1
375	#endif
376	#endif
377	#if defined(INCA_AER) && !defined(INCA_CH4) && !defined(INCA_NMHC)
378	. 1,1,1,1,1,1,1,1,1,1,
379	. 1,1,1,1,1,1,1,1,1,1,
380	. 1,1,1,1,1,1,1,1,1
381	#endif
382	#if defined(INCA) && !defined(INCA_CH4) && !defined(INCA_NMHC) && !defined(INCA_AER)
383	. 1,1,1,1,1,1,1,1,1,1,
384	. 1
385	#endif
386
387	. /)
388
389
390	REAL, PARAMETER :: dry_mass = 28.966
391	REAL, POINTER :: hbuf(:)
392	REAL, ALLOCATABLE :: obuf(:)
393	REAL :: calday
394	REAL :: pdel(klon,klev)
395	REAL :: dummy(klon,klev) = 0.
396	#endif
397	#ifdef INCA_AER
398	integer la
399	#endif
400	c
401	c======================================================================
402	modname='phytrac'
403
404	ps(:)=paprs(:,1)
405
406	if (debutphy) then
407
408	c ecrit_tra = NINT(86400./pdtphys *ecritphy)
409	print*,'dans phytrac ',pdtphys,ecritphy,ecrit_tra
410
411	if(nbtr.lt.nqmax) then
412	c print*,'NQMAX=',nqmax
413	c print*,'NBTR=',nbtr
414	abort_message='See above'
415	call abort_gcm(modname,abort_message,1)
416	endif
417
418	inirnpb=rnpb
419	PRINT*, 'La frequence de sortie traceurs est ', ecrit_tra
420	C
421	c=============================================================
422	c Initialisation des sorties
423	c=============================================================
424
425	#ifdef CPP_IOIPSL
426	#include "ini_histrac.h"
427	#endif
428
429	c======================================================================
430	c Initialisation de certaines variables pour le Rn et le Pb
431	c======================================================================
432
433	c Initialisation du traceur dans le sol (couche limite radonique)
434	c
435	c print*,'valeur de debut dans phytrac :',debutphy
436	trs(:,:) = 0.
437
438	open (99,file='starttrac',status='old',
439	. err=999,form='formatted')
440	read(99,*) (trs(i,1),i=1,klon)
441	999 close(99)
442	c print*, 'apres starttrac'
443
444	c Initialisation de la fraction d'aerosols lessivee
445	c
446	d_tr_lessi_impa(:,:,:) = 0.
447	d_tr_lessi_nucl(:,:,:) = 0.
448	c
449	c Initialisation de la nature des traceurs
450	c
451	DO it = 1, nqmax
452	aerosol(it) = .FALSE. ! Tous les traceurs sont des gaz par defaut
453	radio(it) = .FALSE. ! Par defaut pas de passage par radiornpb
454	clsol(it) = .FALSE. ! Par defaut couche limite avec flux prescrit
455	ENDDO
456	c
457	ENDIF ! fin debutphy
458	c Initialisation du traceur dans le sol (couche limite radonique)
459	if(inirnpb) THEN
460	c
461	radio(1)= .true.
462	radio(2)= .true.
463	clsol(1)= .true.
464	clsol(2)= .true.
465	aerosol(2) = .TRUE. ! le Pb est un aerosol
466	c
467	call initrrnpb (ftsol,pctsrf,masktr,fshtr,hsoltr,tautr
468	. ,vdeptr,scavtr)
469	inirnpb=.false.
470	endif
471	#ifdef INCA
472	!======================================================================
473	! Chimie
474	!======================================================================
475
476	calday = FLOAT(julien) + gmtime
477	ncsec = NINT (86400.*gmtime)
478
479	DO k = 1, nlev
480	pdel(:,k) = paprs(:,k) - paprs (:,k+1)
481	END DO
482
483	#ifdef INCAINFO
484	PRINT *, 'CHEMMAIN @ ', calday, ' ... '
485	DO it = 1, nbtr
486	PRINT *, solsym(it), MINVAL(tr_seri(:,:,it)),
487	$ MAXVAL(tr_seri(:,:,it))
488	END DO
489	#endif
490
491
492	#ifdef INCA_AER
493
494	! Changement Anne 01/04/2005
495	CALL aerosolmain (tr_seri,
496	$ pdtphys,
497	$ pplay,
498	$ pdel,
499	$ prfl,
500	$ pmflxr,
501	$ psfl,
502	$ pmflxs,
503	$ pmfu,
504	$ itop_con,
505	$ ibas_con,
506	$ pphi,
507	$ airephy, ! paire,
508	$ nstep,
509	$ rneb, ! for chimiaq
510	$ t_seri, ! for chimiaq
511	$ rh)
512	! fin changement anne
513
514	#endif
515
516	CALL chemmain (tr_seri, !mmr
517	$ nas, !nas
518	$ nstep, !nstep
519	$ calday, !calday
520	$ julien, !ncdate
521	$ ncsec, !ncsec
522	$ 1, !lat
523	$ pdtphys, !delt
524	$ paprs(1,1), !ps
525	$ pplay, !pmid
526	$ pdel, !pdel
527	$ airephy,
528	$ pctsrf(1,1),!oro
529	$ ftsol, !tsurf
530	$ albsol, !albs
531	$ pphi, !zma
532	$ pphis, !phis
533	$ cldfra, !cldfr
534	$ rneb, !cldfr_st
535	$ diafra, !cldfr_cv
536	$ itop_con, !cldtop
537	$ ibas_con, !cldbot
538	$ cldliq, !cwat
539	$ prfl, !flxrst
540	$ pmflxr, !flxrcv
541	$ psfl, !flxsst
542	$ pmflxs, !flxscv
543	$ pmfu, !flxupd
544	$ flxmass_w, !flxmass_w
545	$ t_seri, !tfld
546	$ sh, !sh
547	$ rh, !rh
548	$ .false., !wrhstts
549	$ hbuf, !hbuf
550	$ obuf, !obuf
551	$ iip1, !nx
552	$ jjp1) !ny
553	#ifdef INCAINFO
554	PRINT *, 'OK.'
555	DO it = 1, nbtr
556	PRINT *, solsym(it), MINVAL(tr_seri(:,:,it)),
557	$ MAXVAL(tr_seri(:,:,it))
558	END DO
559	#endif
560	#else
561
562	c Abder
563	ctestmaf if(nqmax.gt.2) aerosol(3)=.true.
564
565	do i=1,nlon
566	pftsol1(i) = ftsol(i,1)
567	pftsol2(i) = ftsol(i,2)
568	pftsol3(i) = ftsol(i,3)
569	pftsol4(i) = ftsol(i,4)
570
571	ppsrf1(i) = pctsrf(i,1)
572	ppsrf2(i) = pctsrf(i,2)
573	ppsrf3(i) = pctsrf(i,3)
574	ppsrf4(i) = pctsrf(i,4)
575
576	enddo
577	c Abder
578	#endif
579	c======================================================================
580	c Calcul de l'effet de la convection
581	c======================================================================
582	c print*,'Avant convection'
583	c do it=1,nqmax
584	c WRITE(itn,'(i2)') it
585	c call diagtracphy(tr_seri(:,:,it),paprs,'Avant conv'//itn)
586	c enddo
587
588	if (convection) then
589
590	c print*,'Pas de temps dans phytrac : ',pdtphys
591	DO it=1, nqmax
592	#ifdef INCA
593	IF ( conv_flg(it) == 0 ) CYCLE
594	#endif
595	if (iflag_con.eq.2) then
596	c tiedke
597	CALL nflxtr(pdtphys, pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
598	. pplay, paprs, tr_seri(1,1,it), d_tr_cv(1,1,it))
599	else if (iflag_con.eq.3) then
600	c KE
601	call cvltr(pdtphys, da, phi, mp, paprs,pplay, tr_seri(1,1,it),
602	. upwd,dnwd,d_tr_cv(1,1,it))
603	endif
604
605	DO k = 1, nlev
606	DO i = 1, klon
607	tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cv(i,k,it)
608	ENDDO
609	ENDDO
610	#ifdef INCA
611	CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'convection it = '
612	. //solsym(it))
613	#else
614	CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'convection it = '//itn)
615	#endif
616	ENDDO
617	c print*,'apres nflxtr'
618
619	endif ! convection
620	c print*,'Apres convection'
621	c do it=1,nqmax
622	c WRITE(itn,'(i1)') it
623	c call diagtracphy(tr_seri(:,:,it),paprs,'Avant conv'//itn)
624	c enddo
625
626
627	c======================================================================
628	c Calcul de l'effet des thermiques
629	c======================================================================
630
631	do k=1,klev
632	do i=1,klon
633	zmasse(i,k)=(paprs(i,k)-paprs(i,k+1))/rg
634	enddo
635	enddo
636
637	c print*,'masse dans ph ',zmasse
638	do it=1,nqmax
639	do k=1,klev
640	do i=1,klon
641	d_tr_th(i,k,it)=0.
642	tr_seri(i,k,it)=max(tr_seri(i,k,it),0.)
643	tr_seri(i,k,it)=min(tr_seri(i,k,it),1.e10)
644	enddo
645	enddo
646	enddo
647
648	if (iflag_thermals.gt.0) then
649	c print*,'calcul de leffet des thermiques'
650	nsplit=10
651	DO it=1, nqmax
652	c WRITE(itn,'(i1)') it
653	c CALL minmaxqfi(tr_seri(1,1,it),1.e10,-1.e33,'conv it='//itn)
654	c print*,'avant dqthermiquesretro'
655	c call dump2d(iim,jjm-1,tr_seri(2,1,1),'TR_SERI ')
656
657	do isplit=1,nsplit
658	c Abderr 25 11 02
659	C Thermiques
660	c print*,'Avant dans phytrac'
661	call dqthermcell(klon,klev,pdtphys/nsplit
662	. ,fm_therm,entr_therm,zmasse
663	. ,tr_seri(1:klon,1:klev,it),d_tr,ztra_th)
664
665	do k=1,klev
666	do i=1,klon
667	d_tr(i,k)=pdtphys*d_tr(i,k)/nsplit
668	d_tr_th(i,k,it)=d_tr_th(i,k,it)+d_tr(i,k)
669	tr_seri(i,k,it)=max(tr_seri(i,k,it)+d_tr(i,k),0.)
670	enddo
671	enddo
672	enddo ! nsplit
673	c print*,'apres thermiques'
674	c call dump2d(iim,jjm-1,d_tr_th(1,1,1),'d_tr_th ')
675	c do k=1,klev
676	c print*,'d_tr_th(',k,')=',tr_seri(280,k,1)
677	c enddo
678
679	c WRITE(itn,'(i1)') it
680	c CALL minmaxqfi(tr_seri(1,1,it),1.e10,-1.e33,'therm it='//itn)
681	ENDDO ! it
682	endif ! Thermiques
683	c print*,'ATTENTION: sdans thermniques'
684
685	c======================================================================
686	c Calcul de l'effet de la couche limite
687	c======================================================================
688	c print *,'Avant couchelimite'
689	c do it=1,nqmax
690	c WRITE(itn,'(i1)') it
691	c call diagtracphy(tr_seri(:,:,it),paprs,'Avant CL '//itn)
692	c enddo
693
694	if (couchelimite) then
695
696	DO k = 1, nlev
697	DO i = 1, klon
698	delp(i,k) = paprs(i,k)-paprs(i,k+1)
699	ENDDO
700	ENDDO
701
702	C maf modif pour tenir compte du cas rnpb + traceur
703	DO it=1, nqmax
704	#ifdef INCA
705	IF ( pbl_flg(it) == 0 ) CYCLE
706	#endif
707	c print *,'it',it,clsol(it)
708	if (clsol(it)) then ! couche limite avec quantite dans le sol calculee
709	CALL cltracrn(it, pdtphys, yu1, yv1,
710	e coefh,t_seri,ftsol,pctsrf,
711	e tr_seri(1,1,it),trs(1,it),
712	e paprs, pplay, delp,
713	e masktr(1,it),fshtr(1,it),hsoltr(it),
714	e tautr(it),vdeptr(it),
715	e xlat,
716	s d_tr_cl(1,1,it),d_trs)
717	DO k = 1, nlev
718	DO i = 1, klon
719	tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cl(i,k,it)
720	ENDDO
721	ENDDO
722	c
723	c Traceur ds sol
724	c
725	DO i = 1, klon
726	trs(i,it) = trs(i,it) + d_trs(i)
727	END DO
728	C
729	C maf provisoire suppression des prints
730	C WRITE(itn,'(i1)') it
731	C CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'cltracrn it='//itn)
732	else ! couche limite avec flux prescrit
733	#ifdef INCA
734	DO k = 1, klon
735	source(k) = eflux(k,it)-dflux(k,it)
736	END DO
737	#else
738
739	Cmaf provisoire source / traceur a creer
740	DO i=1, klon
741	source(i) = 0.0 ! pas de source, pour l'instant
742	ENDDO
743	C
744	#endif
745	CALL cltrac(pdtphys, coefh,t_seri,
746	s tr_seri(1,1,it), source,
747	e paprs, pplay, delp,
748	s d_tr_cl(1,1,it))
749	DO k = 1, nlev
750	DO i = 1, klon
751	tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cl(i,k,it)
752	ENDDO
753	ENDDO
754	Cmaf WRITE(itn,'(i1)') it
755	cmaf CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'cltracn it='//itn)
756	endif
757	ENDDO
758	c
759	endif ! couche limite
760
761	c print*,'Apres couchelimite'
762	c do it=1,nqmax
763	c WRITE(itn,'(i1)') it
764	c call diagtracphy(tr_seri(:,:,it),paprs,'Avant CL '//itn)
765	c enddo
766
767	c======================================================================
768	c Calcul de l'effet du puits radioactif
769	c======================================================================
770
771	C MAF il faudrait faire une modification pour passer dans radiornpb
772	c si radio=true mais pour l'instant radiornpb propre au cas rnpb
773	if(rnpb) then
774	print *, 'decroissance radiactive activee'
775	call radiornpb (tr_seri,pdtphys,tautr,d_tr_dec)
776	C
777	DO it=1,nqmax
778	if(radio(it)) then
779	DO k = 1, nlev
780	DO i = 1, klon
781	tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_dec(i,k,it)
782	ENDDO
783	ENDDO
784	WRITE(itn,'(i1)') it
785	CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'puits rn it='//itn)
786	endif
787	ENDDO
788	c
789	endif ! rnpb decroissance radioactive
790	C
791	c======================================================================
792	c Calcul de l'effet de la precipitation
793	c======================================================================
794
795	c print*,'LESSIVAGE =',lessivage
796	IF (lessivage) THEN
797
798	c print*,'avant lessivage'
799
800	d_tr_lessi_nucl(:,:,:) = 0.
801	d_tr_lessi_impa(:,:,:) = 0.
802	flestottr(:,:,:) = 0.
803	c
804	c tendance des aerosols nuclees et impactes
805	c
806	DO it = 1, nqmax
807	IF (aerosol(it)) THEN
808	DO k = 1, nlev
809	DO i = 1, klon
810	d_tr_lessi_nucl(i,k,it) = d_tr_lessi_nucl(i,k,it) +
811	s ( 1 - frac_nucl(i,k) )*tr_seri(i,k,it)
812	d_tr_lessi_impa(i,k,it) = d_tr_lessi_impa(i,k,it) +
813	s ( 1 - frac_impa(i,k) )*tr_seri(i,k,it)
814	ENDDO
815	ENDDO
816	ENDIF
817	ENDDO
818	c
819	c Mises a jour des traceurs + calcul des flux de lessivage
820	c Mise a jour due a l'impaction et a la nucleation
821	c
822	c call dump2d(iim,jjm-1,frac_impa(2:klon-1,10),'FRACIMPA')
823	c call dump2d(iim,jjm-1,frac_nucl(2:klon-1,10),'FRACNUCL')
824	c call dump2d(iim,jjm-1,tr_seri(2:klon-1,10,3),'TRACEUR3')
825	DO it = 1, nqmax
826	c print*,'IT=',it,aerosol(it)
827	IF (aerosol(it)) THEN
828	c print*,'IT=',it,' On lessive'
829	DO k = 1, nlev
830	DO i = 1, klon
831	tr_seri(i,k,it)=tr_seri(i,k,it)
832	s frac_impa(i,k)frac_nucl(i,k)
833	ENDDO
834	ENDDO
835	ENDIF
836	ENDDO
837	c call dump2d(iim,jjm-1,tr_seri(2:klon-1,10,3),'TRACEUR3B')
838	c
839	c Flux lessivage total
840	c
841	DO it = 1, nqmax
842	DO k = 1, nlev
843	DO i = 1, klon
844	flestottr(i,k,it) = flestottr(i,k,it) -
845	s ( d_tr_lessi_nucl(i,k,it) +
846	s d_tr_lessi_impa(i,k,it) ) *
847	s ( paprs(i,k)-paprs(i,k+1) ) /
848	s (RG * pdtphys)
849	ENDDO
850	ENDDO
851	c
852	Cmaf WRITE(itn,'(i1)') it
853	Cmaf CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'tr(lessi) it='//itn)
854	ENDDO
855	c
856	c print*,'apres lessivage'
857	ENDIF
858	Cc
859	DO k = 1, nlev
860	DO i = 1, klon
861	fluxrn(i,k) = flestottr(i,k,1)
862	fluxpb(i,k) = flestottr(i,k,2)
863	rn(i,k) = tr_seri(i,k,1)
864	pb(i,k) = tr_seri(i,k,2)
865	pbnucl(i,k)=d_tr_lessi_nucl(i,k,2)
866	pbimpa(i,k)=d_tr_lessi_impa(i,k,2)
867	ENDDO
868	ENDDO
869
870	c=============================================================
871	c Ecriture des sorties
872	c=============================================================
873
874	#ifdef CPP_IOIPSL
875	#include "write_histrac.h"
876	#endif
877
878	c=============================================================
879
880	if (lafin) then
881	print*, 'c est la fin de la physique'
882	open (99,file='restarttrac', form='formatted')
883	do i=1,klon
884	write(99,*) trs(i,1)
885	enddo
886	PRINT*, 'Ecriture du fichier restarttrac'
887	close(99)
888	else
889	c print*, 'physique pas fini'
890	endif
891
892
893	RETURN
894	END

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

Download in other formats:

Original Format