Context Navigation

phytrac.F @ 633

Last change on this file since 633 was 633, checked in by (none), 19 years ago
This commit was manufactured by cvs2svn to create branch 'LMDZ4_par_0'.
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 (iflag_con,
7	I rnpb,
8	I nstep,
9	I julien,
10	I gmtime,
11	I debutphy,
12	I lafin,
13	I nqmax,
14	I nlon,
15	I nlev,
16	I pdtphys,
17	I u,
18	I v,
19	I t_seri,
20	I paprs,
21	I pplay,
22	I pmfu,
23	I pmfd,
24	I pen_u,
25	I pde_u,
26	I pen_d,
27	I pde_d,
28	I coefh,
29	I fm_therm,
30	I entr_therm,
31	I yu1,
32	I yv1,
33	I ftsol,
34	I pctsrf,
35	I xlat,
36	I frac_impa,
37	I frac_nucl,
38	I xlon,
39	I presnivs,
40	I pphis,
41	I pphi,
42	I albsol,
43	I sh,
44	I rh,
45	I cldfra,
46	I rneb,
47	I diafra,
48	I cldliq,
49	I itop_con,
50	I ibas_con,
51	I pmflxr,
52	I pmflxs,
53	I prfl,
54	I psfl,
55	I da,
56	I phi,
57	I mp,
58	I upwd,
59	I dnwd,
60	#ifdef INCA
61	I flxmass_w,
62	#endif
63	O tr_seri)
64
65	USE ioipsl
66
67	#ifdef INCA
68	USE sflx
69	USE chem_tracnm
70	USE species_names
71	USE chem_mods
72	USE pht_tables, ONLY : jrates
73	USE transport_controls, ONLY : conv_flg, pbl_flg
74	USE airplane_src, ONLY : ptrop
75	USE lightning, ONLY : prod_light
76	#ifdef INCA_AER
77	USE AEROSOL_MOD, only : ntr,trmx,trnx
78	USE AEROSOL_DIAG,only : cla,las,tausum,angst,aload,cload,totaerh2o,tau,
79	$ emiss20,sconc,scavcoef_st,scavcoef_cv
80	$ ,cload05ss ,cload05bc ,cload05pom ,cload05dust ,cload05so4
81	$ ,cload125ss ,cload125bc ,cload125pom ,cload125dust ,cload125so4
82	USE AEROSOL_PROGNOS, ONLY : md,mdw
83	USE AEROSOL_METEO, only : airm
84	#endif
85	#ifdef INCA_NMHC
86	USE RESISTANCE_DIAGNOSE, ONLY : surf_alb, sol_irrad, surf_temp, surf_wind,
87	$ aero_resist, lamin_resist, surf_resist
88	#endif
89	#endif
90	IMPLICIT none
91	c======================================================================
92	c Auteur(s) FH
93	c Objet: Moniteur general des tendances traceurs
94	c
95	cAA Remarques en vrac:
96	cAA--------------------
97	cAA 1/ le call phytrac se fait avec nqmax-2 donc nous avons bien
98	cAA les vrais traceurs (nbtr) dans phytrac (pas la vapeur ni eau liquide)
99	cAA 2/ Le choix du radon et du pb se fait juste avec un data
100	cAA (peu propre). Peut-etre pourrait-on prevoir dans l'avenir
101	cAA une variable "type de traceur"
102	c======================================================================
103	#include "YOMCST.h"
104	#include "dimensions.h"
105	#include "dimphy.h"
106	#include "indicesol.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	INTEGER ecrit_tra
269	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	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	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	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	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