Context Navigation

phytrac.F @ 4289

Last change on this file since 4289 was 682, checked in by lmdzadmin, 19 years ago
Ajout de bilans pour INCA - AC MAF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 29.0 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	#ifdef INCA_AER
555
556	c Budget calculation for aerosol species
557	CALL tbudget(airephy,pdtphys,nstep,tr_seri,.false.)
558
559	c-- summary info----------------------------------------------------------------
560
561	if (MOD(nstep,nint(86400./pdtphys)) .eq. 0) then
562	print *, "global aerosol optical thickness "
563
564	write (form,'(A,I2,A)') "(A,",trnx-trmx+1,"A10)"
565	print form,"lamba [nm] ", (solsym(it),it=trmx,trnx)
566
567	write (form,'(A,I2,A)') "(I11,",trnx-trmx+1,"F10.4)"
568	do i=1,las
569	print form,int(lambda(i)),(sum(tausum(:,i,it)*airephy)/sum(airephy),it=trmx,trnx)
570	enddo
571
572	print ,"global mean angstroem component ", sum(angstairephy)/sum(airephy)
573	endif
574	#endif
575	#endif
576
577	#ifdef INCAINFO
578	PRINT *, 'OK.'
579	DO it = 1, nbtr
580	PRINT *, solsym(it), MINVAL(tr_seri(:,:,it)),
581	$ MAXVAL(tr_seri(:,:,it))
582	END DO
583	#endif
584	#else
585
586	c Abder
587	ctestmaf if(nqmax.gt.2) aerosol(3)=.true.
588
589	do i=1,nlon
590	pftsol1(i) = ftsol(i,1)
591	pftsol2(i) = ftsol(i,2)
592	pftsol3(i) = ftsol(i,3)
593	pftsol4(i) = ftsol(i,4)
594
595	ppsrf1(i) = pctsrf(i,1)
596	ppsrf2(i) = pctsrf(i,2)
597	ppsrf3(i) = pctsrf(i,3)
598	ppsrf4(i) = pctsrf(i,4)
599
600	enddo
601	c Abder
602	#endif
603	c======================================================================
604	c Calcul de l'effet de la convection
605	c======================================================================
606	c print*,'Avant convection'
607	c do it=1,nqmax
608	c WRITE(itn,'(i2)') it
609	c call diagtracphy(tr_seri(:,:,it),paprs,'Avant conv'//itn)
610	c enddo
611
612	if (convection) then
613
614	c print*,'Pas de temps dans phytrac : ',pdtphys
615	DO it=1, nqmax
616	#ifdef INCA
617	IF ( conv_flg(it) == 0 ) CYCLE
618	#endif
619	if (iflag_con.eq.2) then
620	c tiedke
621	CALL nflxtr(pdtphys, pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
622	. pplay, paprs, tr_seri(1,1,it), d_tr_cv(1,1,it))
623	else if (iflag_con.eq.3) then
624	c KE
625	call cvltr(pdtphys, da, phi, mp, paprs,pplay, tr_seri(1,1,it),
626	. upwd,dnwd,d_tr_cv(1,1,it))
627	endif
628
629	DO k = 1, nlev
630	DO i = 1, klon
631	tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cv(i,k,it)
632	ENDDO
633	ENDDO
634	#ifdef INCA
635	CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'convection it = '
636	. //solsym(it))
637	#else
638	CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'convection it = '//itn)
639	#endif
640	ENDDO
641	c print*,'apres nflxtr'
642
643	endif ! convection
644	c print*,'Apres convection'
645	c do it=1,nqmax
646	c WRITE(itn,'(i1)') it
647	c call diagtracphy(tr_seri(:,:,it),paprs,'Avant conv'//itn)
648	c enddo
649
650
651	c======================================================================
652	c Calcul de l'effet des thermiques
653	c======================================================================
654
655	do k=1,klev
656	do i=1,klon
657	zmasse(i,k)=(paprs(i,k)-paprs(i,k+1))/rg
658	enddo
659	enddo
660
661	c print*,'masse dans ph ',zmasse
662	do it=1,nqmax
663	do k=1,klev
664	do i=1,klon
665	d_tr_th(i,k,it)=0.
666	tr_seri(i,k,it)=max(tr_seri(i,k,it),0.)
667	tr_seri(i,k,it)=min(tr_seri(i,k,it),1.e10)
668	enddo
669	enddo
670	enddo
671
672	if (iflag_thermals.gt.0) then
673	c print*,'calcul de leffet des thermiques'
674	nsplit=10
675	DO it=1, nqmax
676	c WRITE(itn,'(i1)') it
677	c CALL minmaxqfi(tr_seri(1,1,it),1.e10,-1.e33,'conv it='//itn)
678	c print*,'avant dqthermiquesretro'
679	c call dump2d(iim,jjm-1,tr_seri(2,1,1),'TR_SERI ')
680
681	do isplit=1,nsplit
682	c Abderr 25 11 02
683	C Thermiques
684	c print*,'Avant dans phytrac'
685	call dqthermcell(klon,klev,pdtphys/nsplit
686	. ,fm_therm,entr_therm,zmasse
687	. ,tr_seri(1:klon,1:klev,it),d_tr,ztra_th)
688
689	do k=1,klev
690	do i=1,klon
691	d_tr(i,k)=pdtphys*d_tr(i,k)/nsplit
692	d_tr_th(i,k,it)=d_tr_th(i,k,it)+d_tr(i,k)
693	tr_seri(i,k,it)=max(tr_seri(i,k,it)+d_tr(i,k),0.)
694	enddo
695	enddo
696	enddo ! nsplit
697	c print*,'apres thermiques'
698	c call dump2d(iim,jjm-1,d_tr_th(1,1,1),'d_tr_th ')
699	c do k=1,klev
700	c print*,'d_tr_th(',k,')=',tr_seri(280,k,1)
701	c enddo
702
703	c WRITE(itn,'(i1)') it
704	c CALL minmaxqfi(tr_seri(1,1,it),1.e10,-1.e33,'therm it='//itn)
705	ENDDO ! it
706	endif ! Thermiques
707	c print*,'ATTENTION: sdans thermniques'
708
709	c======================================================================
710	c Calcul de l'effet de la couche limite
711	c======================================================================
712	c print *,'Avant couchelimite'
713	c do it=1,nqmax
714	c WRITE(itn,'(i1)') it
715	c call diagtracphy(tr_seri(:,:,it),paprs,'Avant CL '//itn)
716	c enddo
717
718	if (couchelimite) then
719
720	DO k = 1, nlev
721	DO i = 1, klon
722	delp(i,k) = paprs(i,k)-paprs(i,k+1)
723	ENDDO
724	ENDDO
725
726	C maf modif pour tenir compte du cas rnpb + traceur
727	DO it=1, nqmax
728	#ifdef INCA
729	IF ( pbl_flg(it) == 0 ) CYCLE
730	#endif
731	c print *,'it',it,clsol(it)
732	if (clsol(it)) then ! couche limite avec quantite dans le sol calculee
733	CALL cltracrn(it, pdtphys, yu1, yv1,
734	e coefh,t_seri,ftsol,pctsrf,
735	e tr_seri(1,1,it),trs(1,it),
736	e paprs, pplay, delp,
737	e masktr(1,it),fshtr(1,it),hsoltr(it),
738	e tautr(it),vdeptr(it),
739	e xlat,
740	s d_tr_cl(1,1,it),d_trs)
741	DO k = 1, nlev
742	DO i = 1, klon
743	tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cl(i,k,it)
744	ENDDO
745	ENDDO
746	c
747	c Traceur ds sol
748	c
749	DO i = 1, klon
750	trs(i,it) = trs(i,it) + d_trs(i)
751	END DO
752	C
753	C maf provisoire suppression des prints
754	C WRITE(itn,'(i1)') it
755	C CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'cltracrn it='//itn)
756	else ! couche limite avec flux prescrit
757	#ifdef INCA
758	DO k = 1, klon
759	source(k) = eflux(k,it)-dflux(k,it)
760	END DO
761	#else
762
763	Cmaf provisoire source / traceur a creer
764	DO i=1, klon
765	source(i) = 0.0 ! pas de source, pour l'instant
766	ENDDO
767	C
768	#endif
769	CALL cltrac(pdtphys, coefh,t_seri,
770	s tr_seri(1,1,it), source,
771	e paprs, pplay, delp,
772	s d_tr_cl(1,1,it))
773	DO k = 1, nlev
774	DO i = 1, klon
775	tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cl(i,k,it)
776	ENDDO
777	ENDDO
778	Cmaf WRITE(itn,'(i1)') it
779	cmaf CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'cltracn it='//itn)
780	endif
781	ENDDO
782	c
783	endif ! couche limite
784
785	c print*,'Apres couchelimite'
786	c do it=1,nqmax
787	c WRITE(itn,'(i1)') it
788	c call diagtracphy(tr_seri(:,:,it),paprs,'Avant CL '//itn)
789	c enddo
790
791	c======================================================================
792	c Calcul de l'effet du puits radioactif
793	c======================================================================
794
795	C MAF il faudrait faire une modification pour passer dans radiornpb
796	c si radio=true mais pour l'instant radiornpb propre au cas rnpb
797	if(rnpb) then
798	print *, 'decroissance radiactive activee'
799	call radiornpb (tr_seri,pdtphys,tautr,d_tr_dec)
800	C
801	DO it=1,nqmax
802	if(radio(it)) then
803	DO k = 1, nlev
804	DO i = 1, klon
805	tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_dec(i,k,it)
806	ENDDO
807	ENDDO
808	WRITE(itn,'(i1)') it
809	CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'puits rn it='//itn)
810	endif
811	ENDDO
812	c
813	endif ! rnpb decroissance radioactive
814	C
815	c======================================================================
816	c Calcul de l'effet de la precipitation
817	c======================================================================
818
819	c print*,'LESSIVAGE =',lessivage
820	IF (lessivage) THEN
821
822	c print*,'avant lessivage'
823
824	d_tr_lessi_nucl(:,:,:) = 0.
825	d_tr_lessi_impa(:,:,:) = 0.
826	flestottr(:,:,:) = 0.
827	c
828	c tendance des aerosols nuclees et impactes
829	c
830	DO it = 1, nqmax
831	IF (aerosol(it)) THEN
832	DO k = 1, nlev
833	DO i = 1, klon
834	d_tr_lessi_nucl(i,k,it) = d_tr_lessi_nucl(i,k,it) +
835	s ( 1 - frac_nucl(i,k) )*tr_seri(i,k,it)
836	d_tr_lessi_impa(i,k,it) = d_tr_lessi_impa(i,k,it) +
837	s ( 1 - frac_impa(i,k) )*tr_seri(i,k,it)
838	ENDDO
839	ENDDO
840	ENDIF
841	ENDDO
842	c
843	c Mises a jour des traceurs + calcul des flux de lessivage
844	c Mise a jour due a l'impaction et a la nucleation
845	c
846	c call dump2d(iim,jjm-1,frac_impa(2:klon-1,10),'FRACIMPA')
847	c call dump2d(iim,jjm-1,frac_nucl(2:klon-1,10),'FRACNUCL')
848	c call dump2d(iim,jjm-1,tr_seri(2:klon-1,10,3),'TRACEUR3')
849	DO it = 1, nqmax
850	c print*,'IT=',it,aerosol(it)
851	IF (aerosol(it)) THEN
852	c print*,'IT=',it,' On lessive'
853	DO k = 1, nlev
854	DO i = 1, klon
855	tr_seri(i,k,it)=tr_seri(i,k,it)
856	s frac_impa(i,k)frac_nucl(i,k)
857	ENDDO
858	ENDDO
859	ENDIF
860	ENDDO
861	c call dump2d(iim,jjm-1,tr_seri(2:klon-1,10,3),'TRACEUR3B')
862	c
863	c Flux lessivage total
864	c
865	DO it = 1, nqmax
866	DO k = 1, nlev
867	DO i = 1, klon
868	flestottr(i,k,it) = flestottr(i,k,it) -
869	s ( d_tr_lessi_nucl(i,k,it) +
870	s d_tr_lessi_impa(i,k,it) ) *
871	s ( paprs(i,k)-paprs(i,k+1) ) /
872	s (RG * pdtphys)
873	ENDDO
874	ENDDO
875	c
876	Cmaf WRITE(itn,'(i1)') it
877	Cmaf CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'tr(lessi) it='//itn)
878	ENDDO
879	c
880	c print*,'apres lessivage'
881	ENDIF
882	Cc
883	DO k = 1, nlev
884	DO i = 1, klon
885	fluxrn(i,k) = flestottr(i,k,1)
886	fluxpb(i,k) = flestottr(i,k,2)
887	rn(i,k) = tr_seri(i,k,1)
888	pb(i,k) = tr_seri(i,k,2)
889	pbnucl(i,k)=d_tr_lessi_nucl(i,k,2)
890	pbimpa(i,k)=d_tr_lessi_impa(i,k,2)
891	ENDDO
892	ENDDO
893
894	c=============================================================
895	c Ecriture des sorties
896	c=============================================================
897
898	#ifdef CPP_IOIPSL
899	#include "write_histrac.h"
900	#endif
901
902	c=============================================================
903
904	if (lafin) then
905	print*, 'c est la fin de la physique'
906	open (99,file='restarttrac', form='formatted')
907	do i=1,klon
908	write(99,*) trs(i,1)
909	enddo
910	PRINT*, 'Ecriture du fichier restarttrac'
911	close(99)
912	else
913	c print*, 'physique pas fini'
914	endif
915
916
917	RETURN
918	END

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

Download in other formats:

Original Format