Context Navigation

ch4cloud.F90 @ 3607

Last change on this file since 3607 was 3585, checked in by debatzbr, 3 weeks ago
Connecting microphysics to radiative transfer + miscellaneous cleans
File size: 11.9 KB

Rev	Line
[3585]	1	SUBROUTINE ch4cloud(ngrid,nlay,ptimestep, &
[3247]	2	pplev,pplay,pdpsrf,pzlev,pzlay,pt,pdt, &
	3	pq,pdq,pdqcloud,pdqscloud,pdtcloud, &
	4	nq,rice_ch4)
	5
	6	use comgeomfi_h
	7	use comcstfi_mod, only: pi, g, cpp
	8	use tracer_h, only: igcm_ch4_gas, igcm_ch4_ice, rho_ch4_ice, lw_ch4
	9	use callkeys_mod, only: Nmix_ch4
	10
	11	IMPLICIT NONE
	12
	13	!=======================================================================
	14	! Treatment of saturation of METHANE
	15	!
	16	!
	17	! Modif de zq si saturation dans l'atmosphere
	18	! si zq(ig,l)> zqsat(ig,l) -> zq(ig,l)=zqsat(ig,l)
	19	! Le test est effectue de bas en haut. CH4 condensee
	20	! (si saturation) est remise dans la couche en dessous.
	21	! Le methane condensee dans la couche du bas est deposee a la surface
	22	!
	23	! Melanie Vangvichith (adapted from Mars water clouds)
	24	! Tanguy Bertrand
	25	! Completely rewritten by Francois Forget (to properly estimate the
	26	! latent heat release.
	27	!
	28	!=======================================================================
	29
	30	!-----------------------------------------------------------------------
	31	! declarations:
	32	! -------------
	33
	34	! Inputs:
	35	! ------
	36
	37	INTEGER ngrid,nlay
	38	REAL ptimestep ! pas de temps physique (s)
	39	REAL pplev(ngrid,nlay+1) ! pression aux inter-couches (Pa)
	40	REAL pplay(ngrid,nlay) ! pression au milieu des couches (Pa)
	41	REAL pdpsrf(ngrid) ! tendance surf pressure
	42	REAL pzlev(ngrid,nlay+1) ! altitude at layer boundaries
	43	REAL pzlay(ngrid,nlay) ! altitude at the middle of the layers
	44	REAL pt(ngrid,nlay) ! temperature at the middle of the layers (K)
	45	REAL pdt(ngrid,nlay) ! tendance temperature des autres param.
	46
	47	real pq(ngrid,nlay,nq) ! traceur (kg/kg)
	48	real pdq(ngrid,nlay,nq) ! tendance avant condensation (kg/kg.s-1)
	49	integer nq ! nombre de traceurs
	50
	51	! Outputs:
	52	! -------
	53
	54	real pdqcloud(ngrid,nlay,nq) ! tendance de la condensation CH4(kg/kg.s-1)
	55	real pdqscloud(ngrid,nq) ! flux en surface (kg.m-2.s-1)
	56	REAL pdtcloud(ngrid,nlay) ! tendance temperature due
	57	! a la chaleur latente
	58	REAL rice_ch4(ngrid,nlay) ! CH4 Ice mass mean radius (m)
	59	! (r_c in montmessin_2004)
	60
	61	! local:
	62	! ------
	63	! REAL Nmix ! Cloud condensation nuclei
	64	! parameter (Nmix=1.E5) ! /kg
	65	real rnuclei ! Nuclei geometric mean radius (m)
	66	parameter (rnuclei=2.E-7) ! m
	67
	68	REAL CBRT
	69	EXTERNAL CBRT
	70	INTEGER ig,l
	71
	72
	73	REAL zq(ngrid,nlay,nq) ! local value of tracers
	74	REAL zqsat(ngrid,nlay) ! saturation
	75	REAL zt(ngrid,nlay) ! local value of temperature
	76
	77	REAL vecnull(ngrid*nlay)
	78
	79	LOGICAL,SAVE :: firstcall=.true.
	80
	81	! indexes of methane gas, methane ice and dust tracers:
	82	INTEGER,SAVE :: i_ch4=0 ! methane gas
	83	INTEGER,SAVE :: i_ice=0 ! methane ice
	84
	85	REAL Tfin,qch4_fin
	86	REAL temp_fin ! function used to compute Tfin at the end of the timestep
	87
	88	! TEMPORAIRE :
	89
	90	real pdtcloudmax,pdtcloudmin
	91	integer igmin, igmax,lmin,lmax
	92
	93
	94
	95	! ** un petit test de coherence
	96	! --------------------------
	97
	98	IF (firstcall) THEN
	99	IF(ngrid.NE.ngrid) THEN
	100	PRINT*,'STOP dans ch4cloud'
	101	PRINT*,'probleme de dimensions :'
	102	PRINT*,'ngrid =',ngrid
	103	PRINT*,'ngrid =',ngrid
	104	STOP
	105	ENDIF
	106
	107	if (nq.gt.nq) then
	108	write(,) 'stop in ch4cloud (nq.gt.nq)!'
	109	write(,) 'nq=',nq,' nq=',nq
	110	stop
	111	endif
	112
	113	i_ch4=igcm_ch4_gas
	114	i_ice=igcm_ch4_ice
	115
[3252]	116	if (i_ch4.eq.0) then
	117	write(,) "stop in ch4cloud: i_ch4=0. "
	118	write(,) "Maybe put ch4 in traceur.def?"
	119	stop
	120	endif
	121
[3247]	122	write(,) "methanecloud: i_ch4=",i_ch4
	123	write(,) " i_ice=",i_ice
	124
	125	firstcall=.false.
	126	ENDIF ! of IF (firstcall)
	127
	128
	129	!-----------------------------------------------------------------------
	130	! 1. initialisation
	131	! -----------------
	132
	133	! On "update" la valeur de q(nq) (methane vapor) et temperature.
	134	! On effectue qqes calculs preliminaires sur les couches :
	135
	136	do l=1,nlay
	137	do ig=1,ngrid
	138	zt(ig,l)=pt(ig,l)+ pdt(ig,l)*ptimestep
	139	zq(ig,l,i_ch4)=pq(ig,l,i_ch4)+pdq(ig,l,i_ch4)*ptimestep
	140	zq(ig,l,i_ch4)=max(zq(ig,l,i_ch4),1.E-30)
	141	zq(ig,l,i_ice)=pq(ig,l,i_ice)+pdq(ig,l,i_ice)*ptimestep
	142	zq(ig,l,i_ice)=max(zq(ig,l,i_ice),0.)
	143	enddo
	144	enddo
	145
	146	pdqscloud(1:ngrid,1:nq)=0
	147	pdqcloud(1:ngrid,1:nlay,1:nq)=0
	148	pdtcloud(1:ngrid,1:nlay)=0
	149
	150	! ----------------------------------------------
	151	!
	152	!
	153	! q à saturation dans la couche l à temperature prédite zt
	154	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	155	call methanesat(ngrid*nlay,zt,pplay,zqsat,vecnull)
	156
	157
	158	! Loop to work where CH4 condensation/sublimation is occuring
	159	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	160	do l=1,nlay
	161	do ig=1,ngrid
	162	if ((zq(ig,l,i_ch4).gt.zqsat(ig,l)).or. &
	163	(zq(ig,l,i_ice).gt.1.E-10) ) then ! phase change
	164
	165	if(zq(ig,l,i_ice).lt.(zqsat(ig,l)-zq(ig,l,i_ch4))) then
	166	!case when everything is sublimed:
	167	qch4_fin = zq(ig,l,i_ch4) + zq(ig,l,i_ice)
	168	Tfin = zt(ig,l) - (qch4_fin-zq(ig,l,i_ch4))*lw_ch4/cpp
	169	else
	170	!case when CH4 ice is present at the end
	171	qch4_fin = zqsat(ig,l)
	172	Tfin = zt(ig,l) - (qch4_fin-zq(ig,l,i_ch4))*lw_ch4/cpp
	173	if(abs(Tfin-zt(ig,l)).gt.1.) then
	174	! Improved calculation if the temperature change is significant (1K?)
	175
	176	! Estimating the temperature Tfin and condensation at the end of the
	177	! timestep due to condensation-sublimation process, taking into acount
	178	! latent heat, given by: Tfin - zt = (zq - qsat(Tfin))*L/cp
	179	! Tfin, solution, of the equation, is estimated by iteration
	180	! We assume that improved Tfin is between zt and the 1st estimation of Tfin:
	181	Tfin = temp_fin(zt(ig,l),Tfin,0.01,pplay(ig,l), &
	182	zt(ig,l),zq(ig,l,i_ch4),qch4_fin)
	183	! Security check for small changes
	184	if(abs(zq(ig,l,i_ch4)-qch4_fin).gt. &
	185	abs(zq(ig,l,i_ch4)- zqsat(ig,l)) ) then
	186	write(,) "warning: inaccuracy in CH4 clouds"
	187	write(,) 'zq(ig,l,i_ch4),qch4_fin, zqsat(ig,l)', &
	188	zq(ig,l,i_ch4),qch4_fin, zqsat(ig,l)
	189
	190	qch4_fin = zqsat(ig,l)
	191
	192	end if
	193	end if
	194	end if
	195	! Final tendancies
	196	pdqcloud(ig,l,i_ch4)=(qch4_fin-zq(ig,l,i_ch4))/ptimestep
	197
	198	!TB16
	199	IF (zq(ig,l,i_ch4)+pdqcloud(ig,l,i_ch4)*ptimestep.lt.0.) &
	200	then
	201	pdqcloud(ig,l,i_ch4)=-zq(ig,l,i_ch4)/ptimestep
	202	END IF
	203
	204	pdqcloud(ig,l,i_ice) = - pdqcloud(ig,l,i_ch4)
	205	pdtcloud(ig,l)=(Tfin - zt(ig,l))/ptimestep
	206
	207	! Ice particle radius
	208	zq(ig,l,i_ice)= &
	209	zq(ig,l,i_ice)+pdqcloud(ig,l,i_ice)*ptimestep
	210	rice_ch4(ig,l)=max( CBRT ( (zq(ig,l,i_ice)/rho_ch4_ice &
	211	+Nmix_ch4(4./3.)pirnuclei3.)/(Nmix_ch44./3.*pi)), &
	212	rnuclei)
	213	end if
	214	enddo ! of do ig=1,ngrid
	215	enddo ! of do l=1,nlay
	216
	217	! A correction if a lot of subliming N2 fills the 1st layer FF04/2005
	218	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	219	! Then that should not affect the ice particle radius
	220	do ig=1,ngrid
	221	if(pdpsrf(ig)ptimestep.gt.0.9(pplev(ig,1)-pplev(ig,2)))then
	222	if(pdpsrf(ig)ptimestep.gt.0.9(pplev(ig,1)-pplev(ig,3))) &
	223	rice_ch4(ig,2)=rice_ch4(ig,3)
	224	rice_ch4(ig,1)=rice_ch4(ig,2)
	225	end if
	226	end do
	227
	228	!********** ANALYSE pour Icarus **************
	229	!
	230	! pdtcloudmax =0.
	231	! pdtcloudmin =0.
	232	! igmin=999
	233	! igmax=999
	234	! lmin =999
	235	! lmax =999
	236	! do l=1, nlay
	237	! do ig=1,ngrid
	238	! if(pdtcloud(ig,l).gt.pdtcloudmax) then
	239	! igmax=ig
	240	! lmax = l
	241	! pdtcloudmax=pdtcloud(ig,l)
	242	! else if(pdtcloud(ig,l).lt.pdtcloudmin) then
	243	! igmin=ig
	244	! lmin = l
	245	! pdtcloudmin=pdtcloud(ig,l)
	246	! end if
	247	! end do
	248	! end do
	249	!
	250	! write(,)
	251	! write(,) "MAX , ig, l, Tbefore, Tafter, Tfin_ini"
	252	! write(,) igmax,lmax, zt(igmax,lmax),
	253	! & zt(igmax,lmax)+ pdtcloud(igmax,lmax)*ptimestep,
	254	! &zt(igmax,lmax)-(zqsat(igmax,lmax)-zq(igmax,lmax,i_ch4))*lw_ch4/cpp
	255	! write(,) "MAX , Qch4gas_before, Qch4gas_after, Qsat ini"
	256	! write(,) pq(igmax,lmax,i_ch4)+pdq(igmax,lmax,i_ch4)*ptimestep,
	257	! & pq(igmax,lmax,i_ch4)+
	258	! & (pdq(igmax,lmax,i_ch4)+ pdqcloud(igmax,lmax,i_ch4))*ptimestep
	259	! & ,zqsat(igmax,lmax)
	260	! write(,) "MAX , Qch4ice_before, Qch4ice_after"
	261	! write(,) pq(igmax,lmax,i_ice)+pdq(igmax,lmax,i_ice)*ptimestep,
	262	! & pq(igmax,lmax,i_ice)+
	263	! & (pdq(igmax,lmax,i_ice)+ pdqcloud(igmax,lmax,i_ice))*ptimestep
	264	! write(,) 'Rice=', rice_ch4(igmax,lmax)
	265	!
	266	! write(,)
	267	! write(,) "MIN , ig, l, Tbefore, Tafter, Tfin_ini"
	268	! write(,) igmin,lmin, zt(igmin,lmin),
	269	! & zt(igmin,lmin)+ pdtcloud(igmin,lmin)*ptimestep,
	270	! &zt(igmin,lmin)-(zqsat(igmin,lmin)-zq(igmin,lmin,i_ch4))*lw_ch4/cpp
	271	! write(,) "MIN , Qch4gas_before, Qch4gas_after, Qsat ini"
	272	! write(,) pq(igmin,lmin,i_ch4)+pdq(igmin,lmin,i_ch4)*ptimestep,
	273	! & pq(igmin,lmin,i_ch4)+
	274	! & (pdq(igmin,lmin,i_ch4)+ pdqcloud(igmin,lmin,i_ch4))*ptimestep
	275	! & ,zqsat(igmin,lmin)
	276	! write(,) "MIN , Qch4ice_before, Qch4ice_after"
	277	! write(,) pq(igmin,lmin,i_ice)+pdq(igmin,lmin,i_ice)*ptimestep,
	278	! & pq(igmin,lmin,i_ice)+
	279	! & (pdq(igmin,lmin,i_ice)+ pdqcloud(igmin,lmin,i_ice))*ptimestep
	280	! write(,) 'Rice=', rice_ch4(igmin,lmin)
	281	! write(,) "pdtcloud(igmin,lmin) " , pdtcloud(igmin,lmin)
	282	! write(,) "pdqcloud(i_ch4)", pdqcloud(igmin,lmin,i_ch4)
	283	! write(,) "pdqcloud(i_ice)", pdqcloud(igmin,lmin,i_ice)
	284	!
	285
	286	!********** FIN ANALYSE pour Icarus **************
	287
	288	!**************************************************
	289	! Output --- removed
	290	!**************************************************
	291	! NB: for diagnostics use zq(), the updated value of tracers
	292	! Computing ext visible optical depth in each layer
	293
	294	RETURN
	295	END
	296
	297	!================================================================================
	298
	299	FUNCTION temp_fin(X1,X2,XACC,pres,temp,zq,qsat)
	300	use comcstfi_mod, only: pi, g, cpp
	301	use tracer_h, only: igcm_ch4_gas, igcm_ch4_ice, rho_ch4_ice, lw_ch4
	302	implicit none
	303
	304	real pres, temp, zq,qsat
	305	real X1,X2,XACC, F
	306	real FMID,D,temp_fin,DX, XMID
	307	integer JMAX,J
	308	real func
	309
	310	PARAMETER (JMAX=40)
	311
	312	call methanesat(1,X2,pres,qsat,0.)
	313	FMID = X2 -temp -(zq-qsat)*lw_ch4/cpp
	314
	315	call methanesat(1,X1,pres,qsat,0.)
	316	F = X1 -temp -(zq-qsat)*lw_ch4/cpp
	317
	318	IF(FFMID.GE.0.) write(,*) 'Fix Tfin firt guesses in ch4cloud'
	319	IF(F.LT.0.)THEN
	320	temp_fin=X1
	321	DX=X2-X1
	322	ELSE
	323	temp_fin=X2
	324	DX=X1-X2
	325	ENDIF
	326	DO 11 J=1,JMAX
	327	DX=DX*.5
	328	XMID=temp_fin+DX
	329
	330	call methanesat(1,XMID,pres,qsat,0.)
	331	FMID = XMID -temp -(zq-qsat)*lw_ch4/cpp
	332
	333	IF(FMID.LE.0.)temp_fin=XMID
	334	IF(ABS(DX).LT.XACC .OR. FMID.EQ.0.) RETURN
	335	11 CONTINUE
	336	! PAUSE 'too many bisections'
	337	END
	338
	339

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

Download in other formats:

Original Format