Context Navigation

phytrac.F @ 955

Last change on this file since 955 was 955, checked in by lsce, 16 years ago

ACo+JG

Ajout du flag aerosol_couple :
false=lecture des sulfates dans un fichier(par defaut) - configuration existante
true=calcul des aerosol par INCA

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 29.7 KB

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

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

Download in other formats:

Original Format