Context Navigation

phytrac.F @ 616

Last change on this file since 616 was 616, checked in by lmdzadmin, 19 years ago
Mise a jour pour INCA.2.0 Anne C MAFi+LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 27.3 KB

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

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

Download in other formats:

Original Format