Context Navigation

rain.F90 @ 3562

Last change on this file since 3562 was 2871, checked in by jleconte, 2 years ago
Improved rain evaporation for les + correction in callkeys
File size: 15.0 KB

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

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

Download in other formats:

Original Format