Context Navigation

phytrac.F @ 1000

Last change on this file since 1000 was 888, checked in by slebonnois, 12 years ago
SL: small modifications to the tools, to Venus default .def files and to outputs (including forgotten modifications linked to the 1D); + bug corrections in phytitan
File size: 31.8 KB

Rev	Line
[3]	1	SUBROUTINE phytrac (firstcall,lastcall,
	2	. nqmax,nmicro,ptimestep,appkim,dtkim,
	3	. pplev,pplay,delp,ptemp,pmu0,pfract,pdecli,
[175]	4	. lonsol,
	5	. pu,pv,pzlev,pzlay,ftsol,
	6	. tr_seri,qaer,d_tr_mph,d_tr_kim,
	7	. fclat,reservoir)
[3]	8
	9	c======================================================================
	10	c S. Lebonnois, mai 2008
	11	c
	12	c Arguments:
	13	c
	14	c firstcall----input-L-variable logique indiquant le premier passage
	15	c lastcall-----input-L-variable logique indiquant le dernier passage
	16	c nqmax--------input-I-nombre de traceurs (total)
	17	c nmicro-------input-I-nombre de traceurs microphysiques !! doivent etre toujours en premiers!!
	18	c ptimestep----input-R-pas d'integration pour la physique (seconde)
	19	c appkim-------input-I-appel a la chimie
	20	c dtkim--------input-R-pas de temps chimique (seconde)
	21	c pplev--------input-R-pression pour chaque inter-couche (en Pa)
	22	c pplay--------input-R-pression pour chaque couche (en Pa)
	23	c delp---------input-R-epaisseur d'une couche (en Pa)
	24	c ptemp--------input-R-temperature (K)
	25	c pmu0---------input-R-cos angle zenithal
	26	c pfract-------input-R-fractional day
	27	c pdecli-------input-R-declinaison en radian
	28	c lonsol-------input-R-longitude solaire en radian
[175]	29	c pu-----------input-R-vitesse dans la direction X (de O a E) en m/s (1ere couche)
	30	c pv-----------input-R-vitesse Y (de S a N) en m/s (1ere couche)
	31	c pzlev--------input-R-altitude pour chaque inter-couche (en m)
	32	c pzlay--------input-R-altitude pour chaque couche (en m)
	33	c ftsol--------input-R-temperature au sol (en K)
[3]	34	c tr_seri------input-R-mass mixing ratio traceurs (kg/kg)
	35	c d_tr_mph----output-R-tendance microphysique de "qx" (kg/kg/s)
	36	c d_tr_kim----output-R-tendance chimique de "qx" (kg/kg/s)
[175]	37	c fclat--------output-R-flux de chaleur latente d'evaporation du reservoir CH4 (J/m2/s)
	38	c reservoir----outpur-R-un reservoir de surface !!! (m)
[3]	39	c======================================================================
[102]	40	USE infotrac
	41	use dimphy
	42	IMPLICIT none
[3]	43	#include "dimensions.h"
	44	#include "clesphys.h"
	45	#include "YOMCST.h"
[175]	46	#include "microtab.h"
	47	#include "varmuphy.h"
	48	#include "diagmuphy.h"
	49	#include "itemps.h"
[3]	50
	51	c======================================================================
	52	c Variables argument:
	53	c
	54	LOGICAL firstcall,lastcall
	55	INTEGER nqmax,nmicro,nlat,appkim
	56	REAL ptimestep,dtkim
	57	REAL pplev(klon,klev+1),pplay(klon,klev+1),delp(klon,klev)
	58	REAL ptemp(klon,klev)
	59	REAL pmu0(klon), pfract(klon), pdecli, lonsol
[175]	60	REAL pu(klon),pv(klon)
	61	REAL pzlev(klon,klev+1),pzlay(klon,klev)
	62	REAL ftsol(klon)
[3]	63	REAL tr_seri(klon,klev,nqmax)
[104]	64	REAL qaer(klon,klev,nqmax)
[3]	65	REAL d_tr_mph(klon,klev,nqmax),d_tr_kim(klon,klev,nqmax)
[175]	66	REAL fclat(klon)
	67	REAL reservoir(klon)
[3]	68
	69	c======================================================================
	70	c Local variables
[175]	71	REAL qaer0(klon,klev,nqmax)
[474]	72	REAL prec(klon,5)
[3]	73
[175]	74	c ASTUCE POUR EVITER klon... EN ATTENDANT MIEUX
	75	INTEGER ngrid,NLAYER
	76	PARAMETER (ngrid=(jjm-1)*iim+2) ! = klon
	77	PARAMETER (NLAYER=llm) ! = klev
	78	* common relatifs au nuages
	79	real rmcbar(ngrid,NLAYER),xfbar(ngrid,NLAYER,4)
	80	integer ncount(ngrid,NLAYER)
	81	COMMON/rnuabar/ncount,rmcbar,xfbar
	82
	83	REAL rcloud(klon,klev,nrad),xfrac(klon,klev,4)
	84
	85	REAL vcl,nuc,xgsn,xmsn,xesn,xasn
	86
	87
	88	ReAL gaz1(klon,klev),gaz2(klon,klev),gaz3(klon,klev)
	89
	90	REAL socccld
	91
[3]	92	c grandeurs en moyennes zonales
[175]	93	REAL zplev(jjm+1,klev+1),zplay(jjm+1,klev),ztsol(jjm+1)
	94	REAL zzlev(jjm+1,klev+1),zzlay(jjm+1,klev)
[3]	95	REAL ztemp(jjm+1,klev),zmu0(jjm+1),zfract(jjm+1)
	96	real temp_eq(klev),press_eq(klev)
	97	REAL zqaer(jjm+1,klev,nqmax) ! et non nmicro... Permet nmicro=0.
	98	REAL zqaer0(jjm+1,klev,nqmax)
[175]	99	REAL zdqmufi(jjm+1,klev,nqmax)
[3]	100	REAL ychim(jjm+1,klev,nqmax-nmicro)
[175]	101	REAL zgaz1(jjm+1,klev),zgaz2(jjm+1,klev),zgaz3(jjm+1,klev)
	102	REAL zgaz10(jjm+1,klev),zgaz20(jjm+1,klev),zgaz30(jjm+1,klev)
[3]	103	c La saturation n est calculee qu une seule fois: sauvegarde qysat
	104	c La chimie n est pas calculee tous les pas, il faut donc
	105	c sauvegarder les sorties de la chimie
[104]	106	REAL,save,allocatable :: qysat(:,:),pdyfi(:,:,:)
[3]	107
[104]	108	character*10 nomqy(nqmax-nmicro+1)
[175]	109	integer i,j,k,l,iq,ig0
[104]	110
[474]	111	REAL zprec(jjm+1,5),zsolesp(jjm+1,klev,3),
[175]	112	& zflxesp_i(jjm+1,klev,3)
	113	REAL ztau_drop(jjm+1,klev),ztau_aer(jjm+1,klev,nrad)
	114	c
	115	c indice des esp chimiques utilisees dans la microfi
	116	integer icldch4,icldc2h6,icldc2h2
	117	save icldch4,icldc2h6,icldc2h2
	118
	119	real fte,ftm,Lvch4
	120
	121	REAL tmp,ex,kmin,kmax,dqsq
	122	REAL dqch4
	123	c REAL ch4(jjm+1),ch4b(jjm+1),dch4(jjm+1),ch4c(jjm+1,llm)
	124	c integer ich4
	125	c common/ch4ind/ich4
	126
[3]	127	c======================================================================
	128	c======================================================================
	129
[104]	130	if (firstcall) then
	131	allocate(qysat(klev,nqmax-nmicro),pdyfi(jjm+1,klev,nqmax-nmicro))
[175]	132
	133	c -------- Quelques verifications au demarrage sur les tailles des tableaux.
	134	IF (microfi.ge.1) then
	135	c Faire de la microphysique sans traceurs... bon courage !
	136	if (nmicro.le.0) then
	137	print*,"aLeRtE cRiTiQuE !!!"
	138	print*,"Vous faites de la microphysique sans traceurs"
	139	print*,"microphysique..."
	140	print*,"Je m'arrete et vous laisse reflechir !"
	141	stop
	142	endif
	143	c Nombre de traceurs incompatibles avec la microphysique.
[176]	144	if ((nmicro.ne.ntype*nrad).and.(clouds.eq.1)) then
[175]	145	print*,"aLeRtE cRiTiQuE !!!"
	146	print*,"Nb trac imcompatible avec la microphysique."
	147	print,nmicro,ntypenrad
	148	stop
	149	endif
[176]	150	if ((nmicro.ne.nrad).and.(clouds.eq.0)) then
	151	print*,"aLeRtE cRiTiQuE !!!"
	152	print*,"Nb trac imcompatible avec la microphysique."
	153	print*,nmicro,nrad
	154	stop
	155	endif
[175]	156	ENDIF
	157
[104]	158	endif
	159
[175]	160	c RAZ des sorties : les moyennes se font directement dans IOIPSL :
	161	c
	162	flxesp_i(:,:,:) = 0.
	163	tau_drop(:,:) = 0.
	164	tau_aer(:,:,:) = 0.
	165	solesp(:,:,:) = 0.
[474]	166	precip(:,:) = 0. ! c'est uniquement une sortie en um/s
	167	c
	168	prec(:,:) = 0. ! c'est la variable temporaire des precipitions de la microfi
	169	! prec est en m (metre precipitable)
[175]	170
[3]	171	c-----------------------------------------------------------------------
	172	c convertion moyennes zonales et changement d unites pour microphy
	173	c ---------------------------------
	174
	175	c print*,'CONVERSION 2D ET CHANGEMENT UNITES (PHYTRAC)'
	176
[175]	177	c -------------------
	178	c Gestion de la temperature et de la pression :
	179	c soit la chimie est active, soit la microphysique se fait en 2D.
	180	IF (chimi.or.microfi.eq.1) THEN
	181
	182	zplev = 0.0
	183	zplay = 0.0
	184	zzlev = 0.0
	185	zzlay = 0.0
	186	ztemp = 0.0
	187	zqaer = 0.0
	188	ychim = 0.0
	189	zmu0 = 0.0
	190	zfract= 0.0
	191	zgaz1 = 0.0
	192	zgaz2 = 0.0
	193	zgaz3 = 0.0
	194	zprec = 0.0
	195	zflxesp_i = 0.0
	196	ztau_drop = 0.0
	197	ztau_aer = 0.0
	198	zsolesp = 0.0
	199
[3]	200	do l=1,llm+1
	201	zplev(1,l) = pplev(1,l)
[175]	202	zzlev(1,l) = pzlev(1,l)
[3]	203	do j=2,jjm
	204	ig0=1+(j-2)*iim
	205	do i=1,iim
	206	zplev(j,l) = zplev(j,l) + pplev(ig0+i,l)/iim
[175]	207	zzlev(j,l) = zzlev(j,l) + pzlev(ig0+i,l)/iim
[3]	208	enddo
	209	enddo
	210	zplev(jjm+1,l) = pplev(klon,l)
[175]	211	zzlev(jjm+1,l) = pzlev(klon,l)
[3]	212	enddo
	213
	214	do l=1,llm
	215	ztemp(1,l) = ptemp(1,l)
	216	zplay(1,l) = pplay(1,l)
[175]	217	zzlay(1,l) = pzlay(1,l)
[3]	218	do j=2,jjm
	219	ig0=1+(j-2)*iim
	220	do i=1,iim
	221	ztemp(j,l) = ztemp(j,l) + ptemp(ig0+i,l)/iim
	222	zplay(j,l) = zplay(j,l) + pplay(ig0+i,l)/iim
[175]	223	zzlay(j,l) = zzlay(j,l) + pzlay(ig0+i,l)/iim
[3]	224	enddo
	225	enddo
	226	ztemp(jjm+1,l) = ptemp(klon,l)
	227	zplay(jjm+1,l) = pplay(klon,l)
[175]	228	zzlay(jjm+1,l) = pzlay(klon,l)
[3]	229	temp_eq = ztemp((jjm+1)/2,:)
	230	press_eq = zplay((jjm+1)/2,:)/100. ! en mbar
	231	enddo
	232
[175]	233	ENDIF ! chimi or microfi=1
	234
	235	c -----------------------------
	236	c Gestion des variables de la microphysique :
	237	c
	238	c -------------------
	239	if (microfi.ge.1) then
	240
	241	c Traceurs microphysiques: passage en extensif: n/kg --> n/m^2 (2D ou 3D passage obligatoire)
	242	DO iq=1,nmicro
	243	c print*,tname(iq)
	244	DO l=1,llm
	245	DO i = 1, klon
	246	qaer(i,l,iq) = tr_seri(i,l,iq)*delp(i,l)/RG
	247	ENDDO
	248	ENDDO
	249	ENDDO
	250	c copie du tableau de traceur :
	251	qaer0(:,:,:)=qaer(:,:,:)
	252	c
	253	c -------------------
	254	c Extraction des gaz pour les nuages
	255	c
	256	c recuperation des indices des gaz qui nous interesse
	257	if (firstcall) then
	258	if (clouds.eq.1) then
	259	icldch4=-1
	260	icldc2h6=-1
	261	icldc2h2=-1
	262	do i=1,nqmax
	263	if (tname(i).eq."CH4") then
	264	icldch4=i
	265	c ich4=i
	266	elseif (tname(i).eq."C2H6") then
	267	icldc2h6=i
	268	elseif (tname(i).eq."C2H2") then
	269	icldc2h2=i
	270	endif
	271	enddo
	272	if (icldch4 .eq.-1 .or.
	273	& icldc2h6.eq.-1 .or.
	274	& icldc2h2.eq.-1 ) then
	275	print*, "Sacrebleu !!!"
	276	print*, "Vous voulez faire des nuages sans gaz."
	277	print*, "Mais vous etes inconscient. Je vais m'arreter la"
	278	print*, "pour vous laisser reflechir au probleme"
	279	STOP
	280	endif
	281	endif ! clouds=1
	282	endif ! firstcall
	283
	284	c Saturation et fraction molaire CLOUD
	285	c Calcul des saturations pour les esp chimique de la muphy des nuages.
	286	c On le fait ici pour les sortir dans physiq.F sans avoir a surcharger la routine.
	287	c Elles passent ensuite dans un common pour passer dans les I/O.
	288	c
	289	c-------------------------------------------
	290	IF (clouds.eq.1) THEN
	291	DO l=1,llm
	292	DO i = 1, klon
	293	call ch4sat(ptemp(i,l),pplay(i,l),tmp) !tmp en kg/kg !
	294	satch4(i,l) = tr_seri(i,l,icldch4)/(tmp*28./16.)
	295
	296	call c2h6sat(ptemp(i,l),pplay(i,l),tmp)
	297	satc2h6(i,l) =tr_seri(i,l,icldc2h6)/(tmp*28./30.)
	298
	299	call c2h2sat(ptemp(i,l),pplay(i,l),tmp)
	300	satc2h2(i,l) =tr_seri(i,l,icldc2h2)/(tmp*28./26.)
	301
	302	ENDDO
	303	ENDDO
	304
	305	c Copie des gaz (en 3D) <== UNIQUEMENT SI ON FAIT DES NUAGES
	306	gaz1(:,:) = tr_seri(:,:,icldch4)
	307	gaz2(:,:) = tr_seri(:,:,icldc2h6)
	308	gaz3(:,:) = tr_seri(:,:,icldc2h2)
	309
	310	ENDIF ! clouds=1
	311
	312	endif ! microfi.ge.1
	313
	314	c -------------------
	315	c Si microfi = 1 on est en 2D :
	316	c conversion des inputs de muphys
	317	IF (microfi.eq.1) THEN
	318
[3]	319	zmu0(1) = pmu0(1)
	320	zfract(1) = pfract(1)
	321	do j=2,jjm
	322	ig0=1+(j-2)*iim
	323	do i=1,iim
	324	zmu0(j) = zmu0(j) + pmu0(ig0+i)/iim
	325	zfract(j) = zfract(j) + pfract(ig0+i)/iim
	326	enddo
	327	enddo
	328	zmu0(jjm+1) = pmu0(klon)
	329	zfract(jjm+1) = pfract(klon)
[175]	330	c
	331	c traceurs 3D --> 2D
	332	c
	333	do iq=1,nqmax
[3]	334	do l=1,llm
	335	zqaer(1,l,iq) = qaer(1,l,iq)
	336	do j=2,jjm
	337	ig0=1+(j-2)*iim
	338	do i=1,iim
	339	zqaer(j,l,iq) = zqaer(j,l,iq) + qaer(ig0+i,l,iq)/iim
	340	enddo
	341	enddo
	342	zqaer(jjm+1,l,iq) = qaer(klon,l,iq)
	343	enddo
	344	enddo
[175]	345	c copie du tableau de traceur
	346	zqaer0(:,:,:) = zqaer(:,:,:)
	347	c
	348	c gaz 3D --> 2D <=== UNIQUEMENT SI ON FAIT DES NUAGES.
	349	c
	350	if (clouds.eq.1) then
	351	do l=1,llm
	352	zgaz1(1,l) = gaz1(1,l)
	353	zgaz2(1,l) = gaz2(1,l)
	354	zgaz3(1,l) = gaz3(1,l)
	355	do j=2,jjm
	356	ig0=1+(j-2)*iim
	357	do i=1,iim
	358	zgaz1(j,l) = zgaz1(j,l) + gaz1(ig0+i,l)/iim
	359	zgaz2(j,l) = zgaz2(j,l) + gaz2(ig0+i,l)/iim
	360	zgaz3(j,l) = zgaz3(j,l) + gaz3(ig0+i,l)/iim
	361	enddo
	362	enddo
	363	zgaz1(jjm+1,l) = gaz1(klon,l)
	364	zgaz2(jjm+1,l) = gaz2(klon,l)
	365	zgaz3(jjm+1,l) = gaz3(klon,l)
	366	enddo
	367
	368	zgaz10=zgaz1
	369	zgaz20=zgaz2
	370	zgaz30=zgaz3
	371	endif ! clouds=1
[3]	372
[175]	373	endif ! microfi=1
[3]	374
	375	c AUTRES TRACEURS
	376
	377	if (nqmax.gt.nmicro) then
	378	do iq=nmicro+1,nqmax
	379	do l=1,llm
	380	ychim(1,l,iq-nmicro) = tr_seri(1,l,iq)
	381	do j=2,jjm
	382	ig0=1+(j-2)*iim
	383	do i=1,iim
	384	ychim(j,l,iq-nmicro) = ychim(j,l,iq-nmicro)
	385	. + tr_seri(ig0+i,l,iq)/iim
	386	enddo
	387	enddo
	388	ychim(jjm+1,l,iq-nmicro) = tr_seri(klon,l,iq)
	389	enddo
[175]	390	nomqy(iq-nmicro) = tname(iq)
[3]	391	c print*,iq-nmicro,nomqy(iq-nmicro)
	392	enddo
	393	nomqy(nqmax-nmicro+1) = "HV"
	394	endif
	395
	396	c-----------------------------------------------------------------------
	397	c initialisation des qysat au premier appel:
	398	c ---------------------------------
	399
	400	c!! ATTENTION, qysat pris uniquement a l'equateur
	401	c!! justifie puisque dans cette region, les var de t et p sont faibles...
	402
	403	if (firstcall .and. chimi .and.(nqmax.gt.nmicro)) then
	404	call inicondens(nqmax-nmicro,press_eq,temp_eq,nomqy,qysat)
	405	endif
	406
	407	c-----------------------------------------------------------------------
[175]	408	c Appel de la microphysique en 2D/3D !!!!!!
[3]	409	c --------------------------
	410
	411	IF(firstcall) THEN
	412	print*,'MICROPHYSIQUE ',MICROFI
	413	ENDIF
	414
[175]	415	c call begintime(tt0)
	416	IF (MICROFI.eq.0) THEN
	417	c PAS DE MICROPHYSIQUE :
	418	c On appelle juste rdf pour creer la grille de rayons.
	419	IF (firstcall) THEN
	420	print*,'MICROPHYSIQUE OFF-LINE',MICROFI
	421	call rdf()
[3]	422	ENDIF
[175]	423	c NOTES :
	424	c L'appel de rdf ne sert a rien ici mis a part pour le TR. Si cet
	425	c appel a deja lieu dans le TR inutile de le refaire ici.
	426	c Je ne sais pas exactement comment marche les modules en F90
	427	c Mais je recopie les valeurs du common/part/ de rdf pour
	428	c les mettre dans un common interne a la microphysique (voir varmuphy.h)
	429	c DONC J'AI BESOIN D'AVOIR ACCES A L'ANCIEN COMMON !!!
	430	c
	431	ELSEIF (MICROFI.eq.1) THEN
	432	c MICROPHYSIQUE 2D :
	433	c Les input/output comportent le prefixe z pour 2D :)
	434	zdqmufi = 0. ! ne sert que pour chimi pour condensation
	435	call muphys(jjm+1,
	436	& zplev,zplay,zzlev,zzlay,
	437	& ztemp,zqaer,zgaz1,zgaz2,zgaz3,
	438	& nmicro,ptimestep,
	439	& zmu0,zfract,
	440	c -------- sorties diagnostiques
	441	& zflxesp_i,
	442	& ztau_drop,ztau_aer,
	443	& zsolesp,zprec)
[3]	444	ELSE
[175]	445	c MICROPHYSIQUE 3D :
	446	c Les input sont des champs 3D directement !
	447	call muphys(klon,
	448	& pplev,pplay,pzlev,pzlay,
	449	& ptemp,qaer,gaz1,gaz2,gaz3,
	450	& nmicro,ptimestep,
	451	& pmu0,pfract,
	452	c ------ sorties diagnostiques
	453	& flxesp_i,
	454	& tau_drop,tau_aer,
	455	& solesp,prec)
	456	c
	457	c NOTES :
	458	c Ici toutes nos sorties sont des champs 3D...(meme les diagnostiques)
	459	c On a rien a faire mis a part copier les dq dans les d_tr
	460	c
	461	ENDIF
	462	c call endtime(tt0,tt1)
	463	c ttmuphys=ttmuphys+tt1
[3]	464
[175]	465	c-----------------------------------------------------------------------
	466	c Mise a jour des sorties de muphys
	467	c -------------
	468	c En 2D on copie les sorties de muphys de la grille LATxALT
	469	c sur la grille complete.
	470	IF (microfi.eq.1) THEN
	471	c precipitations
	472	DO l=1,5
	473	prec(1,l) = zprec(1,l)
	474	ig0 = 2
	475	DO j=2,jjm
	476	DO i = 1, iim
	477	prec(ig0,l) = zprec(j,l)
	478	ig0 = ig0 + 1
	479	ENDDO
	480	ENDDO
	481	prec(ig0,l) = zprec(jjm+1,l)
	482	ENDDO
	483	c taux sedimentation
	484	DO l=1,llm
	485	c taux sed goutte
	486	IF (clouds.eq.1) THEN
	487	tau_drop(1,l) = ztau_drop(1,l)
	488	ig0 = 2
	489	DO j=2,jjm
	490	DO i = 1, iim
	491	tau_drop(ig0,l) = ztau_drop(j,l)
	492	ig0 = ig0 + 1
	493	ENDDO
	494	ENDDO
	495	tau_drop(ig0,l) = ztau_drop(jjm+1,l)
	496	ENDIF
	497	c taux sed aer
	498	DO iq=1,nrad
	499	tau_aer(1,l,iq) = ztau_aer(1,l,iq)
	500	ig0 = 2
	501	DO j=2,jjm
	502	DO i = 1, iim
	503	tau_aer(ig0,l,iq) = ztau_aer(j,l,iq)
	504	ig0 = ig0 + 1
	505	ENDDO
	506	ENDDO
	507	tau_aer(ig0,l,iq) = ztau_aer(jjm+1,l,iq)
	508	ENDDO
	509	ENDDO
	510	c flux glace / production glace
	511	IF (clouds.eq.1) THEN
	512	DO iq=1,3
[474]	513	DO l=1,llm
	514	flxesp_i(1,l,iq) = zflxesp_i(1,l,iq)
[175]	515	solesp(1,l,iq) = zsolesp(1,l,iq)
	516	ig0 = 2
	517	DO j=2,jjm
	518	DO i = 1, iim
[474]	519	flxesp_i(ig0,l,iq)=zflxesp_i(j,l,iq)
	520	solesp(ig0,l,iq) = zsolesp(j,l,iq)
[175]	521	ig0 = ig0 + 1
	522	ENDDO
	523	ENDDO
[474]	524	flxesp_i(ig0,l,iq)=zflxesp_i(jjm+1,l,iq)
[175]	525	solesp(ig0,l,iq) = zsolesp(jjm+1,l,iq)
	526	ENDDO
	527	ENDDO
	528	ENDIF
[3]	529	ENDIF
	530
	531	c-----------------------------------------------------------------------
[175]	532	c Gestion des sources
	533	c -------------
	534	c
	535	IF (clouds.eq.1) THEN
	536	IF (microfi.eq.1) THEN
	537	c On repasse les gaz en 3D si on a fait de la microphysique en 2D
	538	DO l=1,llm
	539	gaz1(1,l) = zgaz1(1,l)
	540	gaz2(1,l) = zgaz2(1,l)
	541	gaz3(1,l) = zgaz3(1,l)
	542	ig0 = 2
	543	DO j=2,jjm
	544	DO i = 1, iim
	545	gaz1(ig0,l) = zgaz1(j,l)* gaz1(ig0,l) /zgaz10(j,l)
	546	gaz2(ig0,l) = zgaz2(j,l)* gaz2(ig0,l) /zgaz20(j,l)
	547	gaz3(ig0,l) = zgaz3(j,l)* gaz3(ig0,l) /zgaz30(j,l)
	548	ig0 = ig0 + 1
	549	ENDDO
	550	ENDDO
	551	gaz1(ig0,l) = zgaz1(jjm+1,l)
	552	gaz2(ig0,l) = zgaz2(jjm+1,l)
	553	gaz3(ig0,l) = zgaz3(jjm+1,l)
	554	ENDDO
	555	ENDIF
	556	c Mise a jour du reservoir de CH4 (ie : seul le CH4 remplit le reservoir)
	557	DO i=1,klon
	558	reservoir(i) = reservoir(i)+prec(i,1)
	559	ENDDO
	560	c Calcul des sources :
	561	c ch4=0.
	562	c ch4(1) = gaz1(1,1)
	563	c do j=2,jjm
	564	c ig0=1+(j-2)*iim
	565	c do i=1,iim
	566	c ch4(j)= ch4(j) + gaz1(ig0+i,1)/iim
	567	c enddo
	568	c enddo
	569	c ch4(jjm+1) = gaz1(ig0,1)
	570
	571	CALL sources(klon,klev,ptimestep,z0,
	572	& pu,pv,pplev,pzlay,pzlev,
	573	& gaz1,gaz2,gaz3,
[474]	574	& ftsol,evapch4,reservoir)
[175]	575
	576	c ch4b=0.
	577	c ch4b(1) = gaz1(1,1)
	578	c do j=2,jjm
	579	c ig0=1+(j-2)*iim
	580	c do i=1,iim
	581	c ch4b(j)= ch4b(j) + gaz1(ig0+i,1)/iim
	582	c enddo
	583	c enddo
	584	c ch4b(jjm+1) = gaz1(ig0,1)
	585	c do j=1,jjm+1
	586	c write(499,*) j,ch4(j),ch4b(j)
	587	c enddo
	588	c write(499,*) ""
	589	ENDIF
	590	c-----------------------------------------------------------------------
[3]	591	c Condensation
	592	c -------------
	593
[175]	594	IF ((chimi).and.(nqmax.gt.nmicro)) then
[3]	595
[175]	596	c tendance (en /s) passee sur zdqmufi(nmicro+1 a nqmax)
[3]	597	c print*,'Condensation'
	598
	599	do iq=1,nqmax-nmicro
	600	do l=1,llm
	601	do j=1,jjm+1
	602	if (ychim(j,l,iq).gt.qysat(l,iq)) then
[175]	603	zdqmufi(j,l,nmicro+iq)= (-ychim(j,l,iq)+qysat(l,iq)) !delta y
	604	. / ptimestep ! / dt
[3]	605	endif
	606	enddo
	607	enddo
	608	enddo
	609
	610	ENDIF
	611
	612	c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	613	c eventuellement, modif initiale de la compo
	614	c
[175]	615	c tendance (en /s) passee sur zdqmufi(nmicro+1 a nqmax)
[3]	616	c
	617	c if (firstcall .and. chimi .and.(nqmax.gt.nmicro)) then
	618	c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	619	c!!!remise de CH4 a 1.5%!!!!!!!!!!!!!!!!!!!!!!
	620	c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	621	c do iq=1,nqmax-nmicro
	622	c if (nomqy(iq).eq."CH4") then
	623	c do l=1,llm
	624	c do j=1,jjm+1
	625	c if (ychim(j,l,iq).le.0.015) then
[175]	626	c zdqmufi(j,l,nmicro+iq)= (-ychim(j,l,iq)+0.015) !delta y
[3]	627	c . / ptimestep ! / dt
	628	c endif
	629	c enddo
	630	c enddo
	631	c endif
	632	c enddo
	633	c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	634	c
	635	c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	636	c!!!remise de C2H2 a 1.e-5 max !!!!!!!!!!!!!
	637	c!!!remise de C2H6 a 3.e-5 max !!!!!!!!!!!!!
	638	c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	639	c do iq=1,nqmax-nmicro
	640	c if (nomqy(iq).eq."C2H2") then
	641	c do l=1,llm
	642	c do j=1,jjm+1
	643	c if (ychim(j,l,iq).gt.1.e-5) then
[175]	644	c zdqmufi(j,l,nmicro+iq)= (-ychim(j,l,iq)+1.e-5) !delta y
[3]	645	c . / ptimestep ! / dt
	646	c endif
	647	c enddo
	648	c enddo
	649	c endif
	650	c if (nomqy(iq).eq."C2H6") then
	651	c do l=1,llm
	652	c do j=1,jjm+1
	653	c if (ychim(j,l,iq).gt.3.e-5) then
[175]	654	c zdqmufi(j,l,nmicro+iq)= (-ychim(j,l,iq)+3.e-5) !delta y
[3]	655	c . / ptimestep ! / dt
	656	c endif
	657	c enddo
	658	c enddo
	659	c endif
	660	c enddo
	661	c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	662	c endif
	663	c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	664
[175]	665	c ----- commentaire de fin (mise a jour des profil de fraction molaire)
	666
[3]	667	c-----------------------------------------------------------------------
	668	c Appel de la chimie
	669	c --------------------------
	670
	671	if((appkim.eq.1).and.(chimi)) then
	672	print*,'On passe dans la CHIMIE'
	673
	674	c do iq=1,nqmax-nmicro
	675	c if (nomqy(iq).eq."C2H2") then
	676	c print*,"C2H2top=",ychim(:,klev,iq)
	677	c endif
	678	c enddo
	679
	680	c Appel Chimie
	681	c ------------
	682	CALL calchim(nqmax-nmicro,ychim,nomqy,pdecli,lonsol,dtkim,
	683	. ztemp,zplay,zplev,
	684	. pdyfi)
	685	c ychim ne doit pas etre modifie, pdyfi en /s
	686
	687	endif
	688
	689	c-----------------------------------------------------------------------
	690	c retour des tendances vers 3D
	691	c ---------------------------------
	692
	693	c TRACEURS MICROPHYSIQUES
[175]	694	c
	695	c ---> pas de microphysique
	696	IF (microfi.eq.0) THEN
	697	DO iq=1,nmicro
	698	d_tr_mph(:,:,iq)=0.
	699	ENDDO
	700	ENDIF
	701	c ---> microphysique 2D
	702	IF (microfi.eq.1) THEN
	703	c on repasse le champ de traceurs en 3D (pas les tendances)
	704	DO iq=1,nmicro
	705	DO l=1,llm
	706	qaer(1,l,iq) = zqaer(1,l,iq)
	707	ig0 = 2
	708	DO j=2,jjm
[3]	709	DO i = 1, iim
[175]	710	c un petit patch :
	711	c Si la moyenne zonale au depart est "nulle" :
	712	c On a quand meme le droit de produire des traceurs dans la cellule.
	713	c On considere donc que la valeur de sortie 3D correspond a la valeur de sortie 2D.
[888]	714	c Cela permet aussi entre autre d eviter les NaN pour les traceurs des nuages !
[175]	715	c (au dessus de la tropo pas de nuages donc qaer(nrad+1:ntype*nrad) = 0 !!!)
	716	IF (zqaer0(j,l,iq).lt.1e-100) THEN
	717	qaer(ig0,l,iq) = zqaer(j,l,iq)
	718	ELSE
	719	qaer(ig0,l,iq) = zqaer(j,l,iq) *
	720	& qaer0(ig0,l,iq)/zqaer0(j,l,iq)
	721	ENDIF
	722	ig0 = ig0 + 1
[3]	723	ENDDO
[175]	724	ENDDO
	725	qaer(ig0,l,iq) = zqaer(jjm+1,l,iq)
	726	ENDDO
[3]	727	ENDDO
[175]	728	c La tendances correspond a (qaer-qaer0)/ptimestep
	729	DO iq=1,nmicro
	730	DO i=1,klon
	731	DO l=1,llm
	732	d_tr_mph(i,l,iq) = (qaer(i,l,iq)-qaer0(i,l,iq))/
	733	& ptimestep
	734	ENDDO
	735	ENDDO
	736	ENDDO
	737	c ---> microphysique 3D
	738	ELSEIF(microfi.gt.1) THEN
	739	DO iq=1,nmicro
	740	DO l=1,llm
	741	DO i = 1, klon
	742	d_tr_mph(i,l,iq)=(qaer(i,l,iq)-qaer0(i,l,iq))/ptimestep
	743	ENDDO
	744	ENDDO
	745	ENDDO
[3]	746
[175]	747	ENDIF ! microfi
[3]	748
[808]	749	DO iq=1,nmicro
	750	DO l=1,llm
	751	DO i = 1, klon
	752	c Traceurs microphysiques: passage en intensif: n/m^2 --> n/kg
	753	d_tr_mph(i,l,iq) = d_tr_mph(i,l,iq)*RG/delp(i,l)
	754	ENDDO
	755	ENDDO
	756	ENDDO
	757
[3]	758	c AUTRES TRACEURS
	759
	760	if ((chimi).and.(nqmax.gt.nmicro)) then
	761	c on passe de pdyfi (tendance chimique en /s calculee quand chimie appelee)
	762	c a d_tr_kim (tendance chimique 3D en /s, passee a physiq)
[175]	763	c et de zdqmufi a d_tr_mph (tendance condensation 3D en /s passee a physiq)
[3]	764
	765	DO iq=nmicro+1,nqmax
	766	DO l=1,llm
	767	d_tr_kim(1,l,iq) = pdyfi(1,l,iq-nmicro)
[175]	768	d_tr_mph(1,l,iq) = zdqmufi(1,l,iq)
[3]	769	ig0 = 2
	770	DO j=2,jjm
	771	DO i = 1, iim
	772	d_tr_kim(ig0,l,iq) = pdyfi(j,l,iq-nmicro)
	773	& *tr_seri(ig0,l,iq)/ychim(j,l,iq-nmicro)
[175]	774	d_tr_mph(ig0,l,iq) = zdqmufi(j,l,iq)
[3]	775	& *tr_seri(ig0,l,iq)/ychim(j,l,iq-nmicro)
	776	ig0 = ig0 + 1
	777	ENDDO
	778	ENDDO
	779	d_tr_kim(ig0,l,iq) = pdyfi(jjm+1,l,iq-nmicro)
[175]	780	d_tr_mph(ig0,l,iq) = zdqmufi(jjm+1,l,iq)
[3]	781	ENDDO
	782	ENDDO
	783
	784	endif ! chimi
	785
[175]	786	c--------------------------------------------------
	787	c CONDENSATION VIA MICROFI
	788	c----------------------
[474]	789	c La microphysique avec nuages doit se faire obligatoirement en 3D. (FAUX ACTUELLEMENT)
[888]	790	c Rien n empeche de faire la chimie en 2D. Cependant pour prendre en compte la
[175]	791	c condensation due a la microfi (en 3D) on recalcule la tendance finale pour
	792	c les especes concernees (CH4, C2H6 pour le moment).
	793	IF (microfi.ge.1.and.clouds.eq.1) THEN
	794	DO i=1,klon
	795	DO l=1,klev
	796	c condensation CH4
	797	d_tr_mph(i,l,icldch4)=(gaz1(i,l)-tr_seri(i,l,icldch4))
	798	& /ptimestep
	799	c condensation C2H6
	800	d_tr_mph(i,l,icldc2h6)=(gaz2(i,l)-tr_seri(i,l,icldc2h6))
	801	& /ptimestep
	802	c condensation C2H2
	803	d_tr_mph(i,l,icldc2h2)=(gaz3(i,l)-tr_seri(i,l,icldc2h2))
	804	& /ptimestep
	805	ENDDO
	806	ENDDO
	807	ENDIF
	808	c ch4c=0.
	809	c do l=1,llm
	810	c ch4c(1,l) = tr_seri(1,l,icldch4)
	811	c do j=2,jjm
	812	c ig0=1+(j-2)*iim
	813	c do i=1,iim
	814	c ch4c(j,l)= ch4c(j,l)+tr_seri(ig0+i,l,icldch4)/iim
	815	c enddo
	816	c enddo
	817	c ch4c(jjm+1,l) = tr_seri(klon,l,icldch4)
	818	c enddo
	819	c do l=1,llm
	820	c write(500,*) pplay(25,l),ch4c(25,l)
	821	c enddo
	822	c write(500,*) ""
	823
	824
	825	c--------------------------------------------------
	826	c MISE A JOUR CH4 : (pour refixer la fraction
	827	c molaire)
	828	c--------------------------------------------------
	829	c IF (firstcall) THEN
	830	c do i=1,klon
	831	c do j=1,llm
	832	c call ch4sat(ptemp(i,j),pplay(i,j),tmp) !tmp en kg/kg !
	833	c tmp=0.950.85tmp*28./16.
	834	c if (pplay(i,j).lt.20000.) then
	835	c dqch4 = 1.4e-2
	836	c else
	837	c dqch4 = tmp
	838	c endif
	839	c d_tr_mph(i,j,icldch4)=(-tr_seri(i,j,icldch4)+dqch4)/
	840	c & ptimestep
	841	c enddo
	842	c enddo
	843	c
	844	c ENDIF
	845
	846	c--------------------------------------------------
[474]	847	c CONVERSION PRECIPITATION :
	848	c en microns/secondes
	849	c--------------------------------------------------
	850	precip = prec * 1.e6 / ptimestep
	851
	852
	853	c--------------------------------------------------
[888]	854	c CALCUL DU FLUX DE CHALEUR LATENTE D EVAPORATION
[175]	855	c DU METHANE
	856	c--------------------------------------------------
	857	IF (clouds.eq.1) THEN
	858	DO i=1,klon
	859	fte= (1.-ftsol(i)/305.5)
	860	ftm= (1.-ftsol(i)/190.5)
	861	if(ftm.le.1.e-3) ftm=1.e-3
	862	if(fte.le.1.e-3) fte=1.e-3
	863	Lvch4 =8.314190.4
	864	& (7.08ftm0.354+10.951.1e-2ftm*0.456)
	865	& /mch4
[474]	866	! evapch4 en m3/m2 {ok}
[175]	867	! 425 en kg/m3
	868	! Lv en J/kg {ok}
	869	! ptimestep en s {ok}
[474]	870	fclat(i)=(evapch4(i)Lvch4rhoi_ch4) ! en J/m2/s
[175]	871	ENDDO
	872	ENDIF
	873
	874	c--------------------------------------------------
	875	c GESTION DES RAYONS DE GOUTTES POUR TR
	876	c--------------------------------------------------
	877	IF (clouds.eq.1) THEN
	878
	879	c Calcul du rayon des gouttes par bin ...
	880	c----------------------------------------
	881	DO i=1,klon
	882	DO j=1,klev
	883	DO iq=1,nrad
	884	* Rayon minimum selon la quantité de noyaux
	885	IF (qaer(i,j,iq+nrad) .le. 1.e-5) THEN
	886	rcloud(i,j,iq) = 1.e-10
	887	ELSE
	888	rcloud(i,j,iq)=
	889	& ((qaer(i,j,iq+2*nrad)/qaer(i,j,iq+nrad)+
	890	& qaer(i,j,iq+3*nrad)/qaer(i,j,iq+nrad) +
[474]	891	& v_e(iq))0.75/RPI)*(1./3.)
[175]	892	ENDIF
	893	ENDDO
	894	ENDDO
	895	ENDDO
	896
	897	c .... et de leur rayon moyen total (tt bins confondu)
	898	c------------------------------------------------------
	899	DO i=1,klon
	900	socccld=0.
	901	DO j=klev,1,-1 !de haut en bas pour le calcul des opacites
	902	vcl=0.
	903	nuc=0.
	904	xgsn=0.
	905	xmsn=0.
	906	xesn=0.
	907	xasn=0.
	908	DO iq=1,nrad
	909	vcl=vcl+qaer(i,j,iq+2*nrad)+
	910	& qaer(i,j,iq+3*nrad)+
	911	& qaer(i,j,iq+4*nrad)+
	912	& v_e(iq)*qaer(i,j,iq+nrad) ! volume des gouttes
	913	nuc=nuc+qaer(i,j,iq+nrad) ! nombre de noyaux
	914	xgsn=xgsn+qaer(i,j,iq+nrad)*v_e(iq) ! volume de noyaux
	915	xmsn=xmsn+qaer(i,j,iq+2*nrad) ! volume de methane
	916	xesn=xesn+qaer(i,j,iq+3*nrad) ! volume d' ethane
	917	xasn=xasn+qaer(i,j,iq+4*nrad) ! volume d' acethylene
	918	ENDDO
	919	IF (nuc .le. 1.e-5) THEN
	920	rmcloud(i,j)=1.0e-10
	921	xfrac(i,j,:)=0.
	922	ELSE
	923	IF(xgsn/vcl.lt.0. .or. xgsn/vcl.gt.1.001)
	924	& print*, 'PB AVEC XFRAC:', i,j,xgsn,vcl
	925	rmcloud(i,j)= ! rayon moyen des gouttes
[474]	926	& (vcl/nuc0.75/RPI)*(1./3.)
[175]	927	xfrac(i,j,1)=xgsn/vcl ! fraction volumique noyau/goutte
	928	xfrac(i,j,2)=xmsn/vcl ! fraction volumique CH4/goutte
	929	xfrac(i,j,3)=xesn/vcl ! fraction volumique C2H6/goutte
	930	xfrac(i,j,4)=xasn/vcl ! fraction volumique C2H2/goutte
	931	c calcul du rayon moyen (moyenne temporelle)
	932	rmcbar(i,j)=rmcbar(i,j)+rmcloud(i,j)
	933	xfbar(i,j,:)=xfbar(i,j,:)+xfrac(i,j,:)
	934	ncount(i,j) = ncount(i,j)+1
	935	ENDIF
[474]	936	socccld=socccld+RPI(rmcloud(i,j)2.)nuc
[175]	937	occcld(i,j)=socccld
	938	ENDDO
	939	ENDDO
	940	c
	941	c OCCCLD
	942	c Calcul le nombre d'occurence d'un nuage
[888]	943	c d opacité comprise en kmin et kmax
[175]	944	c k kmin kmax
	945	c 1 0.0000000 0.10000000
	946	c 2 0.10000000 0.17782794
	947	c 3 0.17782794 0.31622776
	948	c 4 0.31622776 0.56234139
	949	c 5 0.56234139 1.0000000
	950	c 6 1.0000000 1.7782795
	951	c 7 1.7782795 3.1622777
	952	c 8 3.1622777 5.6234136
	953	c 9 5.6234136 10.000000
	954	c 10 10.000000 17.782795
	955	c 11 17.782795 31.622778
	956	c 12 31.622778 100000.00
	957	c
[474]	958	c mise a zero de occld_m
	959	occcld_m=0.
[175]	960	DO i=1,klon
	961	DO j=1,klev
	962	DO k=1,12
	963	ex=10.**(0.25)
	964	kmin=0.
	965	kmax=1.e5
	966	if(k.ne.1) kmin=0.1ex*(k-2)
	967	if(k.ne.12) kmax=0.1ex*(k-1)
	968	if(occcld(i,j).ge.kmin .and. occcld(i,j).lt.kmax)
	969	& occcld_m(i,j,k)=1.
	970	ENDDO
	971	ENDDO
	972	ENDDO
	973	ENDIF ! fin condition clouds => pas besoin de calculer des rayons
	974
[3]	975	RETURN
	976	END
[808]	977

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

Download in other formats:

Original Format