Context Navigation

phytrac.F @ 663

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

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

Download in other formats:

Original Format