Context Navigation

phytrac.F @ 4918

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

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

Download in other formats:

Original Format