Context Navigation

rain.F90 @ 1974

Last change on this file since 1974 was 1859, checked in by mturbet, 7 years ago
add diagnostic of reevaporation of precipitation
File size: 14.4 KB

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

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

Download in other formats:

Original Format