Context Navigation

phytrac.F @ 5434

Last change on this file since 5434 was 954, checked in by lsce, 17 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.6 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	c#ifdef CPP_IOIPSL
445	c#include "ini_histrac.h"
446	c#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	c do it=1,nqmax
641	c WRITE(itn,'(i2)') it
642	c call diagtracphy(tr_seri(:,:,it),paprs,'Avant conv'//itn)
643	c 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.eq.2) then
653	c tiedke
654	CALL nflxtr(pdtphys, pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
655	. pplay, paprs, tr_seri(1,1,it), d_tr_cv(1,1,it))
656	else if (iflag_con.eq.3) then
657	c KE
658	call cvltr(pdtphys, da, phi, mp, paprs,pplay, tr_seri(1,1,it),
659	. upwd,dnwd,d_tr_cv(1,1,it))
660	endif
661
662	DO k = 1, nlev
663	DO i = 1, klon
664	tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cv(i,k,it)
665	ENDDO
666	ENDDO
667	#ifdef INCA
668	CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'convection it = '
669	. //solsym(it))
670	#else
671	CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'convection it = '//itn)
672	#endif
673	ENDDO
674	c print*,'apres nflxtr'
675
676	endif ! convection
677	c print*,'Apres convection'
678	c do it=1,nqmax
679	c WRITE(itn,'(i1)') it
680	c call diagtracphy(tr_seri(:,:,it),paprs,'Avant conv'//itn)
681	c enddo
682
683
684	c======================================================================
685	c Calcul de l'effet des thermiques
686	c======================================================================
687
688	do k=1,klev
689	do i=1,klon
690	zmasse(i,k)=(paprs(i,k)-paprs(i,k+1))/rg
691	enddo
692	enddo
693
694	c print*,'masse dans ph ',zmasse
695	do it=1,nqmax
696	do k=1,klev
697	do i=1,klon
698	d_tr_th(i,k,it)=0.
699	tr_seri(i,k,it)=max(tr_seri(i,k,it),0.)
700	tr_seri(i,k,it)=min(tr_seri(i,k,it),1.e10)
701	enddo
702	enddo
703	enddo
704
705	if (iflag_thermals.gt.0) then
706	c print*,'calcul de leffet des thermiques'
707	nsplit=10
708	DO it=1, nqmax
709	c WRITE(itn,'(i1)') it
710	c CALL minmaxqfi(tr_seri(1,1,it),1.e10,-1.e33,'conv it='//itn)
711	c print*,'avant dqthermiquesretro'
712	c call dump2d(iim,jjm-1,tr_seri(2,1,1),'TR_SERI ')
713
714	do isplit=1,nsplit
715	c Abderr 25 11 02
716	C Thermiques
717	c print*,'Avant dans phytrac'
718	call dqthermcell(klon,klev,pdtphys/nsplit
719	. ,fm_therm,entr_therm,zmasse
720	. ,tr_seri(1:klon,1:klev,it),d_tr,ztra_th)
721
722	do k=1,klev
723	do i=1,klon
724	d_tr(i,k)=pdtphys*d_tr(i,k)/nsplit
725	d_tr_th(i,k,it)=d_tr_th(i,k,it)+d_tr(i,k)
726	tr_seri(i,k,it)=max(tr_seri(i,k,it)+d_tr(i,k),0.)
727	enddo
728	enddo
729	enddo ! nsplit
730	c print*,'apres thermiques'
731	c call dump2d(iim,jjm-1,d_tr_th(1,1,1),'d_tr_th ')
732	c do k=1,klev
733	c print*,'d_tr_th(',k,')=',tr_seri(280,k,1)
734	c enddo
735
736	c WRITE(itn,'(i1)') it
737	c CALL minmaxqfi(tr_seri(1,1,it),1.e10,-1.e33,'therm it='//itn)
738	ENDDO ! it
739	endif ! Thermiques
740	c print*,'ATTENTION: sdans thermniques'
741
742	c======================================================================
743	c Calcul de l'effet de la couche limite
744	c======================================================================
745	c print *,'Avant couchelimite'
746	c do it=1,nqmax
747	c WRITE(itn,'(i1)') it
748	c call diagtracphy(tr_seri(:,:,it),paprs,'Avant CL '//itn)
749	c enddo
750
751	if (couchelimite) then
752
753	DO k = 1, nlev
754	DO i = 1, klon
755	delp(i,k) = paprs(i,k)-paprs(i,k+1)
756	ENDDO
757	ENDDO
758
759	C maf modif pour tenir compte du cas rnpb + traceur
760	DO it=1, nqmax
761	#ifdef INCA
762	IF ( pbl_flg(it) == 0 ) CYCLE
763	#endif
764	c print *,'it',it,clsol(it)
765	if (clsol(it)) then ! couche limite avec quantite dans le sol calculee
766	CALL cltracrn(it, pdtphys, yu1, yv1,
767	e coefh,t_seri,ftsol,pctsrf,
768	e tr_seri(1,1,it),trs(1,it),
769	e paprs, pplay, delp,
770	e masktr(1,it),fshtr(1,it),hsoltr(it),
771	e tautr(it),vdeptr(it),
772	e xlat,
773	s d_tr_cl(1,1,it),d_trs)
774	DO k = 1, nlev
775	DO i = 1, klon
776	tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cl(i,k,it)
777	ENDDO
778	ENDDO
779	c
780	c Traceur ds sol
781	c
782	DO i = 1, klon
783	trs(i,it) = trs(i,it) + d_trs(i)
784	END DO
785	C
786	C maf provisoire suppression des prints
787	C WRITE(itn,'(i1)') it
788	C CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'cltracrn it='//itn)
789	else ! couche limite avec flux prescrit
790	#ifndef INCA
791
792	Cmaf provisoire source / traceur a creer
793	DO i=1, klon
794	source(i,it) = 0.0 ! pas de source, pour l'instant
795	ENDDO
796	C
797	#endif
798	CALL cltrac(pdtphys, coefh,t_seri,
799	s tr_seri(1,1,it), source(:,it),
800	e paprs, pplay, delp,
801	s d_tr_cl(1,1,it))
802	DO k = 1, nlev
803	DO i = 1, klon
804	tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cl(i,k,it)
805	ENDDO
806	ENDDO
807	Cmaf WRITE(itn,'(i1)') it
808	cmaf CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'cltracn it='//itn)
809	endif
810	ENDDO
811	c
812	endif ! couche limite
813
814	c print*,'Apres couchelimite'
815	c do it=1,nqmax
816	c WRITE(itn,'(i1)') it
817	c call diagtracphy(tr_seri(:,:,it),paprs,'Avant CL '//itn)
818	c enddo
819
820	c======================================================================
821	c Calcul de l'effet du puits radioactif
822	c======================================================================
823
824	C MAF il faudrait faire une modification pour passer dans radiornpb
825	c si radio=true mais pour l'instant radiornpb propre au cas rnpb
826	if(rnpb) then
827	c print *, 'decroissance radiactive activee'
828	call radiornpb (tr_seri,pdtphys,tautr,d_tr_dec)
829	C
830	DO it=1,nqmax
831	if(radio(it)) then
832	DO k = 1, nlev
833	DO i = 1, klon
834	tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_dec(i,k,it)
835	ENDDO
836	ENDDO
837	WRITE(itn,'(i1)') it
838	CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'puits rn it='//itn)
839	endif
840	ENDDO
841	c
842	endif ! rnpb decroissance radioactive
843	C
844	c======================================================================
845	c Calcul de l'effet de la precipitation
846	c======================================================================
847
848	c print*,'LESSIVAGE =',lessivage
849	IF (lessivage) THEN
850
851	c print*,'avant lessivage'
852
853	d_tr_lessi_nucl(:,:,:) = 0.
854	d_tr_lessi_impa(:,:,:) = 0.
855	flestottr(:,:,:) = 0.
856	c
857	c tendance des aerosols nuclees et impactes
858	c
859	DO it = 1, nqmax
860	IF (aerosol(it)) THEN
861	DO k = 1, nlev
862	DO i = 1, klon
863	d_tr_lessi_nucl(i,k,it) = d_tr_lessi_nucl(i,k,it) +
864	s ( 1 - frac_nucl(i,k) )*tr_seri(i,k,it)
865	d_tr_lessi_impa(i,k,it) = d_tr_lessi_impa(i,k,it) +
866	s ( 1 - frac_impa(i,k) )*tr_seri(i,k,it)
867	ENDDO
868	ENDDO
869	ENDIF
870	ENDDO
871	c
872	c Mises a jour des traceurs + calcul des flux de lessivage
873	c Mise a jour due a l'impaction et a la nucleation
874	c
875	c call dump2d(iim,jjm-1,frac_impa(2:klon-1,10),'FRACIMPA')
876	c call dump2d(iim,jjm-1,frac_nucl(2:klon-1,10),'FRACNUCL')
877	c call dump2d(iim,jjm-1,tr_seri(2:klon-1,10,3),'TRACEUR3')
878	DO it = 1, nqmax
879	c print*,'IT=',it,aerosol(it)
880	IF (aerosol(it)) THEN
881	c print*,'IT=',it,' On lessive'
882	DO k = 1, nlev
883	DO i = 1, klon
884	tr_seri(i,k,it)=tr_seri(i,k,it)
885	s frac_impa(i,k)frac_nucl(i,k)
886	ENDDO
887	ENDDO
888	ENDIF
889	ENDDO
890	c call dump2d(iim,jjm-1,tr_seri(2:klon-1,10,3),'TRACEUR3B')
891	c
892	c Flux lessivage total
893	c
894	DO it = 1, nqmax
895	DO k = 1, nlev
896	DO i = 1, klon
897	flestottr(i,k,it) = flestottr(i,k,it) -
898	s ( d_tr_lessi_nucl(i,k,it) +
899	s d_tr_lessi_impa(i,k,it) ) *
900	s ( paprs(i,k)-paprs(i,k+1) ) /
901	s (RG * pdtphys)
902	ENDDO
903	ENDDO
904	c
905	Cmaf WRITE(itn,'(i1)') it
906	Cmaf CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'tr(lessi) it='//itn)
907	ENDDO
908	c
909	c print*,'apres lessivage'
910	ENDIF
911	Cc
912	DO k = 1, nlev
913	DO i = 1, klon
914	fluxrn(i,k) = flestottr(i,k,1)
915	fluxpb(i,k) = flestottr(i,k,2)
916	rn(i,k) = tr_seri(i,k,1)
917	pb(i,k) = tr_seri(i,k,2)
918	pbnucl(i,k)=d_tr_lessi_nucl(i,k,2)
919	pbimpa(i,k)=d_tr_lessi_impa(i,k,2)
920	ENDDO
921	ENDDO
922
923	c=============================================================
924	c Ecriture des sorties
925	c=============================================================
926
927	c#ifdef CPP_IOIPSL
928	c#include "write_histrac.h"
929	c#endif
930
931	c=============================================================
932
933	if (lafin) then
934	print*, 'c est la fin de la physique'
935	call Gather(trs(:,1),trs_tmp)
936	c$OMP MASTER
937	if (is_mpi_root) then
938
939	open (99,file='restarttrac', form='formatted')
940	do i=1,klon_glo
941	write(99,*) trs_tmp(i)
942	enddo
943	PRINT*, 'Ecriture du fichier restarttrac'
944	close(99)
945	endif
946	c$OMP END MASTER
947	else
948	c print*, 'physique pas fini'
949	endif
950
951
952	RETURN
953	END

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

Download in other formats:

Original Format