Context Navigation

phytrac.F @ 932

Last change on this file since 932 was 928, checked in by lmdzadmin, 18 years ago
Bug variable d_tr_cv non-definie, cas iflag_con different de 2,3 FH IM
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 29.7 KB

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

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

Download in other formats:

Original Format