Context Navigation

phytrac.F @ 958

Last change on this file since 958 was 956, checked in by lmdzadmin, 17 years ago

Nettoyage du controle des parametres physiques. FH

Les parametres cycle_diurne, soil_model, new_oliq, ok_orodr, ok_orolf, ok_limitvrai, nbapp_rad et iflag_con
sont maintenant geres par la physique uniquement.
ecritphy est elimine.
dimphy.F90 et clesphys.h ne sont plus utilises par le code dynamique.
Le test academique obtenu en compilant avec
makegcm -p nophys gcm
fonctionne. FH
IM

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	logical ok_sync
276	parameter (ok_sync = .true.)
277	C
278	C nature du traceur
279	c
280	logical,save,allocatable :: aerosol(:) ! Nature du traceur
281	c ! aerosol(it) = true => aerosol
282	c ! aerosol(it) = false => gaz
283	c ! nat_trac(it) = 1. aerosol
284	logical,save,allocatable :: clsol(:) ! clsol(it) = true => CL sol calculee
285	logical,save,allocatable :: radio(:) ! radio(it)=true => decroisssance radioactive
286	c$OMP THREADPRIVATE(aerosol,clsol,radio)
287	cym save aerosol,clsol,radio
288	C
289	c======================================================================
290	c
291	c Declaration des procedures appelees
292	c
293	c--modif convection tiedtke
294	INTEGER i, k, it
295	INTEGER iq, iiq
296	REAL delp(klon,klev)
297	c--end modif
298	c
299	c Variables liees a l'ecriture de la bande histoire physique
300	c
301	c Variables locales pour effectuer les appels en serie
302	c----------------------------------------------------
303	c
304	REAL d_tr(klon,klev), d_trs(klon) ! tendances de traceurs
305	REAL d_tr_cl(klon,klev,nbtr) ! tendance de traceurs couche limite
306	REAL d_tr_cv(klon,klev,nbtr) ! tendance de traceurs conv pour chq traceur
307	REAL d_tr_th(klon,klev,nbtr) ! la tendance des thermiques
308	REAL d_tr_dec(klon,klev,nbtr) ! la tendance de la decroissance
309	c ! radioactive du rn - > pb
310	REAL d_tr_lessi_impa(klon,klev,nbtr) ! la tendance du lessivage
311	c ! par impaction
312	REAL d_tr_lessi_nucl(klon,klev,nbtr) ! la tendance du lessivage
313	c ! par nucleation
314	REAL fluxrn(klon,klev)
315	REAL fluxpb(klon,klev)
316	REAL pbimpa(klon,klev)
317	REAL pbnucl(klon,klev)
318	REAL rn(klon,klev)
319	REAL pb(klon,klev)
320	REAL flestottr(klon,klev,nbtr) ! flux de lessivage
321	c ! dans chaque couche
322	real zmasse(klon,klev)
323	real ztra_th(klon,klev)
324
325	C
326	character*20 modname
327	character*80 abort_message
328	c
329	c Controles
330	c-------------
331	logical first,couchelimite,convection,lessivage,sorties,
332	s rnpb,inirnpb
333	save first,couchelimite,convection,lessivage,
334	s sorties,inirnpb
335	c$OMP THREADPRIVATE(first,couchelimite,convection,lessivage,
336	c$OMP+ sorties,inirnpb)
337	c data first,couchelimite,convection,lessivage,sorties
338	c s /.true.,.true.,.false.,.true.,.true./
339	c Olivia
340	data first,couchelimite,convection,lessivage,
341	s sorties
342	s /.true.,.true.,.true.,.true.,.true./
343
344
345	#ifdef INCA
346	INTEGER :: lastgas
347	INTEGER :: ncsec
348	INTEGER :: prt_flag_ts(nbtr)
349
350	REAL, PARAMETER :: dry_mass = 28.966
351	REAL, POINTER :: hbuf(:)
352	REAL, ALLOCATABLE :: obuf(:)
353	REAL :: calday
354	REAL :: pdel(klon,klev)
355	REAL :: dummy(klon,klev)
356	#endif
357	c
358	c======================================================================
359
360	#ifdef INCA
361	prt_flag_ts(:)=(/
362	#ifdef INCA_CH4
363	. 1,1,1,0,0,1,1,1,1,1,
364	. 0,1,0,0,0,0,0,1,0,0,
365	. 0,1,1,1,1,0,1,1,1,0,
366	. 1,1,1,1,1,1,1,1,1,1,
367	. 1,0,0
368	#ifdef INCA_AER
369	. ,1,1,1,1,0,1,1,1,1,0,
370	. 1,1,1,1,1,1,0,1,0,1,
371	. 1,1,1,1,0,1,0,1,1,1
372	#endif
373	#endif
374	#ifdef INCA_NMHC
375	. 1,1,1,1,1,1,1,1,1,1,
376	. 1,1,1,1,1,1,1,1,1,1,
377	. 1,1,1,1,1,1,1,1,1,1,
378	. 1,1,1,1,1,1,1,1,1,1,
379	. 1,1,1,1,1,1,1,1,1,1,
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
384	#ifdef INCA_AER
385	. ,1,1,1,1,0,1,1,1,1,0,
386	. 1,1,1,1,1,1,0,1,0,1,
387	. 1,1,1,1,0,1,0,1,1,1
388	#endif
389	#endif
390	#if defined(INCA_AER) && !defined(INCA_CH4) && !defined(INCA_NMHC)
391	. 1,1,1,1,1,1,1,1,1,1,
392	. 1,1,1,1,1,1,1,1,1,1,
393	. 1,1,1,1,1,1,1,1,1
394	#endif
395	#if defined(INCA) && !defined(INCA_CH4) && !defined(INCA_NMHC) && !defined(INCA_AER)
396	. 1,1,1,1,1,1,1,1,1,1,
397	. 1
398	#endif
399
400	. /)
401	dummy(:,:) = 0.
402
403	#endif
404	modname='phytrac'
405
406	ps(:)=paprs(:,1)
407
408	if (debutphy) then
409	allocate( trs(klon,nbtr) )
410	allocate( masktr(klon,nbtr))
411	allocate( fshtr(klon,nbtr) )
412	allocate( hsoltr(nbtr))
413	allocate( tautr(nbtr))
414	allocate( vdeptr(nbtr))
415	allocate( scavtr(nbtr))
416	allocate( aerosol(nbtr))
417	allocate( clsol(nbtr))
418	allocate( radio(nbtr))
419
420
421	c jg: c'est ca qu'on veut?????
422	! FH 2008/05/09 correction de la frequence d'ecriture des traceurs
423	! ecrit_tra = FLOAT(NINT(86400./pdtphys *ecritphy))
424	print*,'dans phytrac ',pdtphys,ecrit_tra
425
426	if(nbtr.lt.nqmax) then
427	c print*,'NQMAX=',nqmax
428	c print*,'NBTR=',nbtr
429	abort_message='See above'
430	call abort_gcm(modname,abort_message,1)
431	endif
432
433	inirnpb=rnpb
434	PRINT*, 'La frequence de sortie traceurs est ', ecrit_tra
435	C
436	c=============================================================
437	c Initialisation des sorties
438	c=============================================================
439
440	#ifdef CPP_IOIPSL
441	#include "ini_histrac.h"
442	#endif
443
444	c======================================================================
445	c Initialisation de certaines variables pour le Rn et le Pb
446	c======================================================================
447
448	c Initialisation du traceur dans le sol (couche limite radonique)
449	c
450	c print*,'valeur de debut dans phytrac :',debutphy
451	trs(:,:) = 0.
452	c$OMP MASTER
453	if (is_mpi_root) then
454	trs_tmp(:)=0.
455	open (99,file='starttrac',status='old',
456	. err=999,form='formatted')
457	read(99,*) (trs_tmp(i),i=1,klon_glo)
458	999 close(99)
459	endif
460	c$OMP END MASTER
461	call Scatter(trs_tmp,trs(:,1))
462
463	c print*, 'apres starttrac'
464
465	c Initialisation de la fraction d'aerosols lessivee
466	c
467	d_tr_lessi_impa(:,:,:) = 0.
468	d_tr_lessi_nucl(:,:,:) = 0.
469	c
470	c Initialisation de la nature des traceurs
471	c
472	DO it = 1, nqmax
473	aerosol(it) = .FALSE. ! Tous les traceurs sont des gaz par defaut
474	radio(it) = .FALSE. ! Par defaut pas de passage par radiornpb
475	clsol(it) = .FALSE. ! Par defaut couche limite avec flux prescrit
476	ENDDO
477	c
478	ENDIF ! fin debutphy
479	c Initialisation du traceur dans le sol (couche limite radonique)
480	if(inirnpb) THEN
481	c
482	radio(1)= .true.
483	radio(2)= .true.
484	clsol(1)= .true.
485	clsol(2)= .true.
486	aerosol(2) = .TRUE. ! le Pb est un aerosol
487	c
488	call initrrnpb (ftsol,pctsrf,masktr,fshtr,hsoltr,tautr
489	. ,vdeptr,scavtr)
490	inirnpb=.false.
491	endif
492	#ifdef INCA
493	call VTe(VTphysiq)
494	call VTb(VTinca)
495	!======================================================================
496	! Chimie
497	!======================================================================
498
499	calday = FLOAT(julien) + gmtime
500	ncsec = NINT (86400.*gmtime)
501
502	DO k = 1, nlev
503	pdel(:,k) = paprs(:,k) - paprs (:,k+1)
504	END DO
505
506	#ifdef INCAINFO
507	PRINT *, 'CHEMMAIN @ ', calday, ' ... '
508	DO it = 1, nbtr
509	PRINT *, solsym(it), MINVAL(tr_seri(:,:,it)),
510	$ MAXVAL(tr_seri(:,:,it))
511	END DO
512	#endif
513
514
515	#ifdef INCA_AER
516	CALL aerosolmain (aerosol_couple,
517	$ tr_seri,
518	$ pdtphys,
519	$ pplay,
520	$ pdel,
521	$ prfl,
522	$ pmflxr,
523	$ psfl,
524	$ pmflxs,
525	$ pmfu,
526	$ itop_con,
527	$ ibas_con,
528	$ pphi,
529	$ airephy, ! paire,
530	$ nstep,
531	$ rneb, ! for chimiaq
532	$ t_seri, ! for chimiaq
533	$ rh,
534	$ tau_inca,
535	$ piz_inca,
536	$ cg_inca,
537	$ rfname,
538	$ ccm,
539	$ lafin)
540	#endif
541
542	CALL chemmain (tr_seri, !mmr
543	$ nstep, !nstep
544	$ calday, !calday
545	$ julien, !ncdate
546	$ ncsec, !ncsec
547	$ 1, !lat
548	$ pdtphys, !delt
549	$ paprs(1,1), !ps
550	$ pplay, !pmid
551	$ pdel, !pdel
552	$ airephy,
553	$ pctsrf(1,1),!oro
554	$ ftsol, !tsurf
555	$ albsol, !albs
556	$ pphi, !zma
557	$ pphis, !phis
558	$ cldfra, !cldfr
559	$ rneb, !cldfr_st
560	$ diafra, !cldfr_cv
561	$ itop_con, !cldtop
562	$ ibas_con, !cldbot
563	$ cldliq, !cwat
564	$ prfl, !flxrst
565	$ pmflxr, !flxrcv
566	$ psfl, !flxsst
567	$ pmflxs, !flxscv
568	$ pmfu, !flxupd
569	$ flxmass_w, !flxmass_w
570	$ t_seri, !tfld
571	$ sh, !sh
572	$ rh, !rh
573	$ .false., !wrhstts
574	$ hbuf, !hbuf
575	$ obuf, !obuf
576	$ iip1, !nx
577	$ jjp1, !ny
578	$ source,
579	$ solsym)
580	#ifdef INCAINFO
581	#ifdef INCA_AER
582
583	c Budget calculation for aerosol species
584	CALL tbudget(airephy,pdtphys,nstep,tr_seri,.false.)
585
586	c-- summary info----------------------------------------------------------------
587
588	if (MOD(nstep,nint(86400./pdtphys)) .eq. 0) then
589	print *, "global aerosol optical thickness "
590
591	write (form,'(A,I2,A)') "(A,",trnx-trmx+1,"A10)"
592	print form,"lamba [nm] ", (solsym(it),it=trmx,trnx)
593
594	write (form,'(A,I2,A)') "(I11,",trnx-trmx+1,"F10.4)"
595	do i=1,las
596	print form,int(lambda(i)),(sum(tausum(:,i,it)*airephy)/sum(airephy),it=trmx,trnx)
597	enddo
598
599	print ,"global mean angstroem component ", sum(angstairephy)/sum(airephy)
600	endif
601	#endif
602	#endif
603
604	#ifdef INCAINFO
605	PRINT *, 'OK.'
606	DO it = 1, nbtr
607	PRINT *, solsym(it), MINVAL(tr_seri(:,:,it)),
608	$ MAXVAL(tr_seri(:,:,it))
609	END DO
610	#endif
611	call VTe(VTinca)
612	call VTb(VTphysiq)
613	#else
614
615	c Abder
616	ctestmaf if(nqmax.gt.2) aerosol(3)=.true.
617
618	do i=1,nlon
619	pftsol1(i) = ftsol(i,1)
620	pftsol2(i) = ftsol(i,2)
621	pftsol3(i) = ftsol(i,3)
622	pftsol4(i) = ftsol(i,4)
623
624	ppsrf1(i) = pctsrf(i,1)
625	ppsrf2(i) = pctsrf(i,2)
626	ppsrf3(i) = pctsrf(i,3)
627	ppsrf4(i) = pctsrf(i,4)
628
629	enddo
630	c Abder
631	#endif
632	c======================================================================
633	c Calcul de l'effet de la convection
634	c======================================================================
635	c print*,'Avant convection'
636	do it=1,nqmax
637	WRITE(itn,'(i2)') it
638	c call diagtracphy(tr_seri(:,:,it),paprs,'Avant conv'//itn)
639	enddo
640
641	if (convection) then
642
643	c print*,'Pas de temps dans phytrac : ',pdtphys
644	DO it=1, nqmax
645	#ifdef INCA
646	IF ( conv_flg(it) == 0 ) CYCLE
647	#endif
648	if (iflag_con.lt.2) then
649	d_tr_cv=0.
650	else if (iflag_con.eq.2) then
651	c tiedke
652	CALL nflxtr(pdtphys, pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
653	. pplay, paprs, tr_seri(1,1,it), d_tr_cv(1,1,it))
654	else
655	c KE
656	call cvltr(pdtphys, da, phi, mp, paprs,pplay, tr_seri(1,1,it),
657	. upwd,dnwd,d_tr_cv(1,1,it))
658	endif
659
660	DO k = 1, nlev
661	DO i = 1, klon
662	tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cv(i,k,it)
663	ENDDO
664	ENDDO
665	#ifdef INCA
666	CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'convection it = '
667	. //solsym(it))
668	#else
669	CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'convection it = '//itn)
670	#endif
671	ENDDO
672	c print*,'apres nflxtr'
673
674	endif ! convection
675	c print*,'Apres convection'
676	c do it=1,nqmax
677	c WRITE(itn,'(i1)') it
678	c call diagtracphy(tr_seri(:,:,it),paprs,'Avant conv'//itn)
679	c enddo
680
681
682	c======================================================================
683	c Calcul de l'effet des thermiques
684	c======================================================================
685
686	do k=1,klev
687	do i=1,klon
688	zmasse(i,k)=(paprs(i,k)-paprs(i,k+1))/rg
689	enddo
690	enddo
691
692	c print*,'masse dans ph ',zmasse
693	do it=1,nqmax
694	do k=1,klev
695	do i=1,klon
696	d_tr_th(i,k,it)=0.
697	tr_seri(i,k,it)=max(tr_seri(i,k,it),0.)
698	tr_seri(i,k,it)=min(tr_seri(i,k,it),1.e10)
699	enddo
700	enddo
701	enddo
702
703	if (iflag_thermals.gt.0) then
704	c print*,'calcul de leffet des thermiques'
705	nsplit=10
706	DO it=1, nqmax
707	c WRITE(itn,'(i1)') it
708	c CALL minmaxqfi(tr_seri(1,1,it),1.e10,-1.e33,'conv it='//itn)
709	c print*,'avant dqthermiquesretro'
710	c call dump2d(iim,jjm-1,tr_seri(2,1,1),'TR_SERI ')
711
712	do isplit=1,nsplit
713	c Abderr 25 11 02
714	C Thermiques
715	c print*,'Avant dans phytrac'
716	call dqthermcell(klon,klev,pdtphys/nsplit
717	. ,fm_therm,entr_therm,zmasse
718	. ,tr_seri(1:klon,1:klev,it),d_tr,ztra_th)
719
720	do k=1,klev
721	do i=1,klon
722	d_tr(i,k)=pdtphys*d_tr(i,k)/nsplit
723	d_tr_th(i,k,it)=d_tr_th(i,k,it)+d_tr(i,k)
724	tr_seri(i,k,it)=max(tr_seri(i,k,it)+d_tr(i,k),0.)
725	enddo
726	enddo
727	enddo ! nsplit
728	c print*,'apres thermiques'
729	c call dump2d(iim,jjm-1,d_tr_th(1,1,1),'d_tr_th ')
730	c do k=1,klev
731	c print*,'d_tr_th(',k,')=',tr_seri(280,k,1)
732	c enddo
733
734	c WRITE(itn,'(i1)') it
735	c CALL minmaxqfi(tr_seri(1,1,it),1.e10,-1.e33,'therm it='//itn)
736	ENDDO ! it
737	endif ! Thermiques
738	c print*,'ATTENTION: sdans thermniques'
739
740	c======================================================================
741	c Calcul de l'effet de la couche limite
742	c======================================================================
743	c print *,'Avant couchelimite'
744	c do it=1,nqmax
745	c WRITE(itn,'(i1)') it
746	c call diagtracphy(tr_seri(:,:,it),paprs,'Avant CL '//itn)
747	c enddo
748
749	if (couchelimite) then
750
751	DO k = 1, nlev
752	DO i = 1, klon
753	delp(i,k) = paprs(i,k)-paprs(i,k+1)
754	ENDDO
755	ENDDO
756
757	C maf modif pour tenir compte du cas rnpb + traceur
758	DO it=1, nqmax
759	#ifdef INCA
760	IF ( pbl_flg(it) == 0 ) CYCLE
761	#endif
762	c print *,'it',it,clsol(it)
763	if (clsol(it)) then ! couche limite avec quantite dans le sol calculee
764	CALL cltracrn(it, pdtphys, yu1, yv1,
765	e coefh,t_seri,ftsol,pctsrf,
766	e tr_seri(1,1,it),trs(1,it),
767	e paprs, pplay, delp,
768	e masktr(1,it),fshtr(1,it),hsoltr(it),
769	e tautr(it),vdeptr(it),
770	e xlat,
771	s d_tr_cl(1,1,it),d_trs)
772	DO k = 1, nlev
773	DO i = 1, klon
774	tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cl(i,k,it)
775	ENDDO
776	ENDDO
777	c
778	c Traceur ds sol
779	c
780	DO i = 1, klon
781	trs(i,it) = trs(i,it) + d_trs(i)
782	END DO
783	C
784	C maf provisoire suppression des prints
785	C WRITE(itn,'(i1)') it
786	C CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'cltracrn it='//itn)
787	else ! couche limite avec flux prescrit
788	#ifndef INCA
789
790	Cmaf provisoire source / traceur a creer
791	DO i=1, klon
792	source(i,it) = 0.0 ! pas de source, pour l'instant
793	ENDDO
794	C
795	#endif
796	CALL cltrac(pdtphys, coefh,t_seri,
797	s tr_seri(1,1,it), source(:,it),
798	e paprs, pplay, delp,
799	s d_tr_cl(1,1,it))
800	DO k = 1, nlev
801	DO i = 1, klon
802	tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cl(i,k,it)
803	ENDDO
804	ENDDO
805	Cmaf WRITE(itn,'(i1)') it
806	cmaf CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'cltracn it='//itn)
807	endif
808	ENDDO
809	c
810	endif ! couche limite
811
812	c print*,'Apres couchelimite'
813	c do it=1,nqmax
814	c WRITE(itn,'(i1)') it
815	c call diagtracphy(tr_seri(:,:,it),paprs,'Avant CL '//itn)
816	c enddo
817
818	c======================================================================
819	c Calcul de l'effet du puits radioactif
820	c======================================================================
821
822	C MAF il faudrait faire une modification pour passer dans radiornpb
823	c si radio=true mais pour l'instant radiornpb propre au cas rnpb
824	if(rnpb) then
825	c print *, 'decroissance radiactive activee'
826	call radiornpb (tr_seri,pdtphys,tautr,d_tr_dec)
827	C
828	DO it=1,nqmax
829	if(radio(it)) then
830	DO k = 1, nlev
831	DO i = 1, klon
832	tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_dec(i,k,it)
833	ENDDO
834	ENDDO
835	WRITE(itn,'(i1)') it
836	CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'puits rn it='//itn)
837	endif
838	ENDDO
839	c
840	endif ! rnpb decroissance radioactive
841	C
842	c======================================================================
843	c Calcul de l'effet de la precipitation
844	c======================================================================
845
846	c print*,'LESSIVAGE =',lessivage
847	IF (lessivage) THEN
848
849	c print*,'avant lessivage'
850
851	d_tr_lessi_nucl(:,:,:) = 0.
852	d_tr_lessi_impa(:,:,:) = 0.
853	flestottr(:,:,:) = 0.
854	c
855	c tendance des aerosols nuclees et impactes
856	c
857	DO it = 1, nqmax
858	IF (aerosol(it)) THEN
859	DO k = 1, nlev
860	DO i = 1, klon
861	d_tr_lessi_nucl(i,k,it) = d_tr_lessi_nucl(i,k,it) +
862	s ( 1 - frac_nucl(i,k) )*tr_seri(i,k,it)
863	d_tr_lessi_impa(i,k,it) = d_tr_lessi_impa(i,k,it) +
864	s ( 1 - frac_impa(i,k) )*tr_seri(i,k,it)
865	ENDDO
866	ENDDO
867	ENDIF
868	ENDDO
869	c
870	c Mises a jour des traceurs + calcul des flux de lessivage
871	c Mise a jour due a l'impaction et a la nucleation
872	c
873	c call dump2d(iim,jjm-1,frac_impa(2:klon-1,10),'FRACIMPA')
874	c call dump2d(iim,jjm-1,frac_nucl(2:klon-1,10),'FRACNUCL')
875	c call dump2d(iim,jjm-1,tr_seri(2:klon-1,10,3),'TRACEUR3')
876	DO it = 1, nqmax
877	c print*,'IT=',it,aerosol(it)
878	IF (aerosol(it)) THEN
879	c print*,'IT=',it,' On lessive'
880	DO k = 1, nlev
881	DO i = 1, klon
882	tr_seri(i,k,it)=tr_seri(i,k,it)
883	s frac_impa(i,k)frac_nucl(i,k)
884	ENDDO
885	ENDDO
886	ENDIF
887	ENDDO
888	c call dump2d(iim,jjm-1,tr_seri(2:klon-1,10,3),'TRACEUR3B')
889	c
890	c Flux lessivage total
891	c
892	DO it = 1, nqmax
893	DO k = 1, nlev
894	DO i = 1, klon
895	flestottr(i,k,it) = flestottr(i,k,it) -
896	s ( d_tr_lessi_nucl(i,k,it) +
897	s d_tr_lessi_impa(i,k,it) ) *
898	s ( paprs(i,k)-paprs(i,k+1) ) /
899	s (RG * pdtphys)
900	ENDDO
901	ENDDO
902	c
903	Cmaf WRITE(itn,'(i1)') it
904	Cmaf CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'tr(lessi) it='//itn)
905	ENDDO
906	c
907	c print*,'apres lessivage'
908	ENDIF
909	Cc
910	DO k = 1, nlev
911	DO i = 1, klon
912	fluxrn(i,k) = flestottr(i,k,1)
913	fluxpb(i,k) = flestottr(i,k,2)
914	rn(i,k) = tr_seri(i,k,1)
915	pb(i,k) = tr_seri(i,k,2)
916	pbnucl(i,k)=d_tr_lessi_nucl(i,k,2)
917	pbimpa(i,k)=d_tr_lessi_impa(i,k,2)
918	ENDDO
919	ENDDO
920
921	c=============================================================
922	c Ecriture des sorties
923	c=============================================================
924
925	#ifdef CPP_IOIPSL
926	#include "write_histrac.h"
927	#endif
928
929	c=============================================================
930
931	if (lafin) then
932	print*, 'c est la fin de la physique'
933	call Gather(trs(:,1),trs_tmp)
934	c$OMP MASTER
935	if (is_mpi_root) then
936
937	open (99,file='restarttrac', form='formatted')
938	do i=1,klon_glo
939	write(99,*) trs_tmp(i)
940	enddo
941	PRINT*, 'Ecriture du fichier restarttrac'
942	close(99)
943	endif
944	c$OMP END MASTER
945	else
946	c print*, 'physique pas fini'
947	endif
948
949
950	RETURN
951	END

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

Download in other formats:

Original Format