Context Navigation

phytrac.F @ 619

Last change on this file since 619 was 619, checked in by lmdzadmin, 20 years ago
Rajout convection Kerry Emanuel pour traceurs- MAF+JYG
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 27.4 KB

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

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

Download in other formats:

Original Format