Context Navigation

rain.F90 @ 1243

Last change on this file since 1243 was 1016, checked in by jleconte, 11 years ago

07/08/2013 == JL

Water cycle in double precision (largescale+moistadj)
Improved wate rayleigh.
First step for rayleigh with variable species. Now, just need to change optcv.
changed some interpolation indices in callcorrk to limit dependency of OLR on the number of layers

File size: 14.0 KB

Rev	Line
[787]	1	subroutine rain(ngrid,nq,ptimestep,pplev,pplay,t,pdt,pq,pdq,d_t,dqrain,dqsrain,dqssnow,rneb)
[135]	2
[253]	3
[728]	4	! to use 'getin'
	5	use ioipsl_getincom
	6	use watercommon_h, only: T_h2O_ice_liq,T_h2O_ice_clouds, RLVTT, RCPD, RCPV, RV, RVTMP2,Psat_water,Tsat_water,rhowater
	7	use radii_mod, only: h2o_cloudrad
[787]	8	USE tracer_h
[135]	9	implicit none
	10
	11	!==================================================================
	12	!
	13	! Purpose
	14	! -------
	15	! Calculates H2O precipitation using simplified microphysics.
	16	!
	17	! Authors
	18	! -------
	19	! Adapted from the LMDTERRE code by R. Wordsworth (2009)
[728]	20	! Added rain vaporization in case of T>Tsat
[135]	21	! Original author Z. X. Li (1993)
	22	!
	23	!==================================================================
	24
	25	#include "dimensions.h"
	26	#include "dimphys.h"
	27	#include "comcstfi.h"
	28	#include "callkeys.h"
	29
[858]	30	! Arguments
	31	integer,intent(in) :: ngrid ! number of atmospherci columns
	32	integer,intent(in) :: nq ! number of tracers
	33	real,intent(in) :: ptimestep ! time interval
	34	real,intent(in) :: pplev(ngrid,nlayermx+1) ! inter-layer pressure (Pa)
	35	real,intent(in) :: pplay(ngrid,nlayermx) ! mid-layer pressure (Pa)
	36	real,intent(in) :: t(ngrid,nlayermx) ! input temperature (K)
	37	real,intent(in) :: pdt(ngrid,nlayermx) ! input tendency on temperature (K/s)
	38	real,intent(in) :: pq(ngrid,nlayermx,nq) ! tracers (kg/kg)
	39	real,intent(in) :: pdq(ngrid,nlayermx,nq) ! input tendency on tracers
	40	real,intent(out) :: d_t(ngrid,nlayermx) ! temperature tendency (K/s)
	41	real,intent(out) :: dqrain(ngrid,nlayermx,nq) ! tendency of H2O precipitation (kg/kg.s-1)
	42	real,intent(out) :: dqsrain(ngrid) ! rain flux at the surface (kg.m-2.s-1)
	43	real,intent(out) :: dqssnow(ngrid) ! snow flux at the surface (kg.m-2.s-1)
	44	real,intent(in) :: rneb(ngrid,nlayermx) ! cloud fraction
[787]	45
	46	REAL zt(ngrid,nlayermx) ! working temperature (K)
	47	REAL ql(ngrid,nlayermx) ! liquid water (Kg/Kg)
	48	REAL q(ngrid,nlayermx) ! specific humidity (Kg/Kg)
	49	REAL d_q(ngrid,nlayermx) ! water vapor increment
	50	REAL d_ql(ngrid,nlayermx) ! liquid water / ice increment
[135]	51
	52	! Subroutine options
[858]	53	REAL,PARAMETER :: seuil_neb=0.001 ! Nebulosity threshold
[135]	54
[728]	55	INTEGER,save :: precip_scheme ! id number for precipitaion scheme
	56	! for simple scheme (precip_scheme=1)
	57	REAL,SAVE :: rainthreshold ! Precipitation threshold in simple scheme
	58	! for sundquist scheme (precip_scheme=2-3)
	59	REAL,SAVE :: cloud_sat ! Precipitation threshold in non simple scheme
	60	REAL,SAVE :: precip_timescale ! Precipitation timescale
	61	! for Boucher scheme (precip_scheme=4)
	62	REAL,SAVE :: Cboucher ! Precipitation constant in Boucher 95 scheme
	63	REAL,PARAMETER :: Kboucher=1.19E8
	64	REAL,SAVE :: c1
[135]	65
[858]	66	INTEGER,PARAMETER :: ninter=5
[135]	67
[1016]	68	logical,save :: evap_prec ! Does the rain evaporate?
[135]	69
	70	! for simple scheme
[858]	71	real,parameter :: t_crit=218.0
[135]	72	real lconvert
	73
	74	! Local variables
	75	INTEGER i, k, n
[787]	76	REAL zqs(ngrid,nlayermx),Tsat(ngrid,nlayermx), zdelta, zcor
	77	REAL zrfl(ngrid), zrfln(ngrid), zqev, zqevt
[135]	78
[787]	79	REAL zoliq(ngrid)
	80	REAL zdz(ngrid),zrho(ngrid),ztot(ngrid), zrhol(ngrid)
	81	REAL zchau(ngrid),zfroi(ngrid),zfrac(ngrid),zneb(ngrid)
[135]	82
[787]	83	real reffh2oliq(ngrid,nlayermx),reffh2oice(ngrid,nlayermx)
[728]	84
	85	real ttemp, ptemp, psat_tmp
[787]	86	real tnext(ngrid,nlayermx)
[135]	87
[787]	88	real l2c(ngrid,nlayermx)
[253]	89	real dWtot
[135]	90
[253]	91
[135]	92	! Indices of water vapour and water ice tracers
	93	INTEGER, SAVE :: i_vap=0 ! water vapour
	94	INTEGER, SAVE :: i_ice=0 ! water ice
	95
[858]	96	LOGICAL,SAVE :: firstcall=.true.
[135]	97
	98	! Online functions
[731]	99	REAL fallv, fall2v, zzz ! falling speed of ice crystals
[135]	100	fallv (zzz) = 3.29 * ((zzz)**0.16)
[731]	101	fall2v (zzz) =10.6 * ((zzz)**0.31) !for use with radii
[135]	102
	103
	104	IF (firstcall) THEN
	105
	106	i_vap=igcm_h2o_vap
	107	i_ice=igcm_h2o_ice
	108
	109	write(,) "rain: i_ice=",i_ice
	110	write(,) " i_vap=",i_vap
	111
	112	PRINT*, 'in rain.F, ninter=', ninter
	113	PRINT*, 'in rain.F, evap_prec=', evap_prec
	114
[728]	115	write(,) "Precipitation scheme to use?"
	116	precip_scheme=1 ! default value
	117	call getin("precip_scheme",precip_scheme)
	118	write(,) " precip_scheme = ",precip_scheme
	119
	120	if (precip_scheme.eq.1) then
	121	write(,) "rainthreshold in simple scheme?"
	122	rainthreshold=0. ! default value
	123	call getin("rainthreshold",rainthreshold)
	124	write(,) " rainthreshold = ",rainthreshold
	125
	126	else if (precip_scheme.eq.2.or.precip_scheme.eq.3) then
	127	write(,) "cloud water saturation level in non simple scheme?"
	128	cloud_sat=2.6e-4 ! default value
	129	call getin("cloud_sat",cloud_sat)
	130	write(,) " cloud_sat = ",cloud_sat
	131	write(,) "precipitation timescale in non simple scheme?"
	132	precip_timescale=3600. ! default value
	133	call getin("precip_timescale",precip_timescale)
	134	write(,) " precip_timescale = ",precip_timescale
	135
	136	else if (precip_scheme.eq.4) then
	137	write(,) "multiplicative constant in Boucher 95 precip scheme"
	138	Cboucher=1. ! default value
	139	call getin("Cboucher",Cboucher)
	140	write(,) " Cboucher = ",Cboucher
	141	c1=1.001.097/rhowaterCboucher*Kboucher
	142
	143	endif
	144
[1016]	145	write(,) "re-evaporate precipitations?"
	146	evap_prec=.true. ! default value
	147	call getin("evap_prec",evap_prec)
	148	write(,) " evap_prec = ",evap_prec
	149
[135]	150	firstcall = .false.
	151	ENDIF
	152
	153	! GCM -----> subroutine variables
	154	DO k = 1, nlayermx
[787]	155	DO i = 1, ngrid
[135]	156
[253]	157	zt(i,k) = t(i,k)+pdt(i,k)*ptimestep ! a big fat bug was here
	158	q(i,k) = pq(i,k,i_vap)+pdq(i,k,i_vap)*ptimestep
	159	ql(i,k) = pq(i,k,i_ice)+pdq(i,k,i_ice)*ptimestep
[135]	160
[253]	161	!q(i,k) = pq(i,k,i_vap)!+pdq(i,k,i_vap)
	162	!ql(i,k) = pq(i,k,i_ice)!+pdq(i,k,i_ice)
	163
[135]	164	if(q(i,k).lt.0.)then ! if this is not done, we don't conserve water
	165	q(i,k)=0.
	166	endif
	167	if(ql(i,k).lt.0.)then
	168	ql(i,k)=0.
	169	endif
	170
	171	ENDDO
	172	ENDDO
	173
	174	! Initialise the outputs
	175	DO k = 1, nlayermx
[787]	176	DO i = 1, ngrid
[135]	177	d_t(i,k) = 0.0
	178	d_q(i,k) = 0.0
	179	d_ql(i,k) = 0.0
	180	ENDDO
	181	ENDDO
[787]	182	DO i = 1, ngrid
[135]	183	zrfl(i) = 0.0
	184	zrfln(i) = 0.0
	185	ENDDO
	186
	187	! calculate saturation mixing ratio
	188	DO k = 1, nlayermx
[787]	189	DO i = 1, ngrid
[253]	190	ttemp = zt(i,k)
[135]	191	ptemp = pplay(i,k)
[728]	192	! call watersat(ttemp,ptemp,zqs(i,k))
	193	call Psat_water(ttemp,ptemp,psat_tmp,zqs(i,k))
	194	call Tsat_water(ptemp,Tsat(i,k))
[135]	195	ENDDO
	196	ENDDO
	197
[253]	198	! get column / layer conversion factor
[135]	199	DO k = 1, nlayermx
[787]	200	DO i = 1, ngrid
[253]	201	l2c(i,k)=(pplev(i,k)-pplev(i,k+1))/g
[135]	202	ENDDO
	203	ENDDO
	204
	205
	206	! Vertical loop (from top to bottom)
	207	! We carry the rain with us and calculate that added by warm/cold precipitation
	208	! processes and that subtracted by evaporation at each level.
	209	DO 9999 k = nlayermx, 1, -1
	210
	211	IF (evap_prec) THEN ! note no rneb dependence!
[787]	212	DO i = 1, ngrid
[135]	213	IF (zrfl(i) .GT.0.) THEN
[253]	214
[728]	215	if(zt(i,k).gt.Tsat(i,k))then
[863]	216	!! treat the case where all liquid water should boil
	217	zqev=MIN((zt(i,k)-Tsat(i,k))RCPDl2c(i,k)/RLVTT/ptimestep,zrfl(i))
[728]	218	zrfl(i)=MAX(zrfl(i)-zqev,0.)
[863]	219	d_q(i,k)=zqev/l2c(i,k)*ptimestep
[728]	220	d_t(i,k) = - d_q(i,k) * RLVTT/RCPD
	221	else
[731]	222	zqev = MAX (0.0, (zqs(i,k)-q(i,k)))*l2c(i,k)/ptimestep !there was a bug here
[728]	223	zqevt= 2.0e-5*(1.0-q(i,k)/zqs(i,k)) & !default was 2.e-5
	224	sqrt(zrfl(i))l2c(i,k)/pplay(i,k)zt(i,k)R ! BC modif here
	225	zqevt = MAX (zqevt, 0.0)
	226	zqev = MIN (zqev, zqevt)
	227	zqev = MAX (zqev, 0.0)
	228	zrfln(i)= zrfl(i) - zqev
	229	zrfln(i)= max(zrfln(i),0.0)
[253]	230
[728]	231	d_q(i,k) = - (zrfln(i)-zrfl(i))/l2c(i,k)*ptimestep
	232	!d_t(i,k) = d_q(i,k) * RLVTT/RCPD!/(1.0+RVTMP2*q(i,k)) ! double BC modif here
	233	d_t(i,k) = - d_q(i,k) * RLVTT/RCPD ! was bugged!
	234	zrfl(i) = zrfln(i)
	235	end if
	236
[135]	237
	238	ENDIF
	239	ENDDO
	240	ENDIF
	241
[787]	242	DO i = 1, ngrid
[135]	243	zoliq(i) = 0.0
	244	ENDDO
	245
	246
[728]	247	if(precip_scheme.eq.1)then
[135]	248
[787]	249	DO i = 1, ngrid
[253]	250	ttemp = zt(i,k)
[650]	251	IF (ttemp .ge. T_h2O_ice_liq) THEN
[253]	252	lconvert=rainthreshold
	253	ELSEIF (ttemp .gt. t_crit) THEN
	254	lconvert=rainthreshold*(1.- t_crit/ttemp)
	255	lconvert=MAX(0.0,lconvert)
	256	ELSE
	257	lconvert=0.
	258	ENDIF
[135]	259
[253]	260
	261	IF (ql(i,k).gt.1.e-9) then
	262	zneb(i) = MAX(rneb(i,k), seuil_neb)
	263	IF ((ql(i,k)/zneb(i)).gt.lconvert)THEN ! precipitate!
[622]	264	d_ql(i,k) = -MAX((ql(i,k)-lconvert*zneb(i)),0.0)
[253]	265	zrfl(i) = zrfl(i) - d_ql(i,k)*l2c(i,k)/ptimestep
	266	ENDIF
	267	ENDIF
[135]	268	ENDDO
	269
[728]	270	elseif (precip_scheme.ge.2) then
	271
[787]	272	DO i = 1, ngrid
[135]	273	IF (rneb(i,k).GT.0.0) THEN
	274	zoliq(i) = ql(i,k)
[253]	275	zrho(i) = pplay(i,k) / ( zt(i,k) * R )
[135]	276	zdz(i) = (pplev(i,k)-pplev(i,k+1)) / (zrho(i)*g)
[650]	277	zfrac(i) = (zt(i,k)-T_h2O_ice_clouds) / (T_h2O_ice_liq-T_h2O_ice_clouds)
[135]	278	zfrac(i) = MAX(zfrac(i), 0.0)
	279	zfrac(i) = MIN(zfrac(i), 1.0)
	280	zneb(i) = MAX(rneb(i,k), seuil_neb)
	281	ENDIF
[731]	282	ENDDO
[135]	283
[731]	284	!recalculate liquid water particle radii
[787]	285	call h2o_cloudrad(ngrid,ql,reffh2oliq,reffh2oice)
[731]	286
[728]	287	SELECT CASE(precip_scheme)
	288	!precip scheme from Sundquist 78
	289	CASE(2)
	290
[135]	291	DO n = 1, ninter
[787]	292	DO i = 1, ngrid
[135]	293	IF (rneb(i,k).GT.0.0) THEN
[728]	294	! this is the ONLY place where zneb, precip_timescale and cloud_sat are used
[253]	295
[728]	296	zchau(i) = (ptimestep/(FLOAT(ninter)precip_timescale)) zoliq(i) &
	297	* (1.0-EXP(-(zoliq(i)/zneb(i)/cloud_sat)*2)) zfrac(i)
[135]	298	zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
	299	zfroi(i) = ptimestep/FLOAT(ninter)/zdz(i)*zoliq(i) &
[731]	300	fall2v(reffh2oice(i,k)) (1.0-zfrac(i)) ! zfroi behaves oddly...
[135]	301	ztot(i) = zchau(i) + zfroi(i)
	302
	303	IF (zneb(i).EQ.seuil_neb) ztot(i) = 0.0
	304	ztot(i) = MIN(MAX(ztot(i),0.0),zoliq(i))
	305	zoliq(i) = MAX(zoliq(i)-ztot(i), 0.0)
[253]	306
[135]	307	ENDIF
	308	ENDDO
	309	ENDDO
	310
[728]	311	!precip scheme modified from Sundquist 78 (in q**3)
	312	CASE(3)
	313
	314	DO n = 1, ninter
[787]	315	DO i = 1, ngrid
[728]	316	IF (rneb(i,k).GT.0.0) THEN
	317	! this is the ONLY place where zneb, precip_timescale and cloud_sat are used
	318
	319	zchau(i) = (ptimestep/(FLOAT(ninter)precip_timescalecloud_sat*2)) (zoliq(i)/zneb(i))**3
	320	zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
	321	zfroi(i) = ptimestep/FLOAT(ninter)/zdz(i)*zoliq(i) &
[731]	322	fall2v(reffh2oice(i,k)) (1.0-zfrac(i)) ! zfroi behaves oddly...
[728]	323	ztot(i) = zchau(i) + zfroi(i)
	324
	325	IF (zneb(i).EQ.seuil_neb) ztot(i) = 0.0
	326	ztot(i) = MIN(MAX(ztot(i),0.0),zoliq(i))
	327	zoliq(i) = MAX(zoliq(i)-ztot(i), 0.0)
	328
	329	ENDIF
	330	ENDDO
	331	ENDDO
	332
	333	!precip scheme modified from Boucher 95
	334	CASE(4)
	335
	336	DO n = 1, ninter
[787]	337	DO i = 1, ngrid
[728]	338	IF (rneb(i,k).GT.0.0) THEN
	339	! this is the ONLY place where zneb and c1 are used
	340
	341	zchau(i) = ptimestep/FLOAT(ninter) c1 zrho(i) &
	342	(zoliq(i)/zneb(i))2reffh2oliq(i,k)zneb(i) zfrac(i)
	343	zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
	344	zfroi(i) = ptimestep/FLOAT(ninter)/zdz(i)*zoliq(i) &
[731]	345	fall2v(reffh2oice(i,k)) (1.0-zfrac(i)) ! zfroi behaves oddly...
[728]	346	ztot(i) = zchau(i) + zfroi(i)
	347
	348	IF (zneb(i).EQ.seuil_neb) ztot(i) = 0.0
	349	ztot(i) = MIN(MAX(ztot(i),0.0),zoliq(i))
	350	zoliq(i) = MAX(zoliq(i)-ztot(i), 0.0)
	351
	352	ENDIF
	353	ENDDO
	354	ENDDO
	355
	356	END SELECT ! precip_scheme
	357
[135]	358	! Change in cloud density and surface H2O values
[787]	359	DO i = 1, ngrid
[135]	360	IF (rneb(i,k).GT.0.0) THEN
[253]	361	d_ql(i,k) = (zoliq(i) - ql(i,k))!/ptimestep
	362	zrfl(i) = zrfl(i)+ MAX(ql(i,k)-zoliq(i),0.0)*l2c(i,k)/ptimestep
[135]	363	ENDIF
	364	ENDDO
	365
	366
[728]	367	endif ! if precip_scheme=1
	368
[135]	369	9999 continue
	370
	371	! Rain or snow on the ground
[787]	372	DO i = 1, ngrid
[253]	373	if(zrfl(i).lt.0.0)then
	374	print*,'Droplets of negative rain are falling...'
	375	call abort
	376	endif
[650]	377	IF (t(i,1) .LT. T_h2O_ice_liq) THEN
[135]	378	dqssnow(i) = zrfl(i)
[253]	379	dqsrain(i) = 0.0
[135]	380	ELSE
[253]	381	dqssnow(i) = 0.0
[135]	382	dqsrain(i) = zrfl(i) ! liquid water = ice for now
	383	ENDIF
	384	ENDDO
	385
	386	! now subroutine -----> GCM variables
	387	DO k = 1, nlayermx
[787]	388	DO i = 1, ngrid
[135]	389
	390	if(evap_prec)then
	391	dqrain(i,k,i_vap) = d_q(i,k)/ptimestep
	392	d_t(i,k) = d_t(i,k)/ptimestep
	393	else
	394	dqrain(i,k,i_vap) = 0.0
	395	d_t(i,k) = 0.0
	396	endif
[253]	397	dqrain(i,k,i_ice) = d_ql(i,k)/ptimestep
[135]	398
	399	ENDDO
	400	ENDDO
	401
	402	RETURN
	403	end subroutine rain

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

Download in other formats:

Original Format