Context Navigation

phytrac.F @ 620

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

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

Download in other formats:

Original Format