Context Navigation

phytrac.F @ 740

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

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

Download in other formats:

Original Format