Context Navigation

phytrac.F @ 658

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

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

Download in other formats:

Original Format