Context Navigation

rain.F90 @ 2613

Last change on this file since 2613 was 1995, checked in by jleconte, 6 years ago
fix a small bug introduced in previous commit
File size: 14.3 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
[1995]	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
[135]	185	ptemp = pplay(i,k)
[1993]	186	call Psat_water(zt(i,k) ,ptemp,psat_tmp,zqs(i,k))
[728]	187	call Tsat_water(ptemp,Tsat(i,k))
[135]	188	ENDDO
	189	ENDDO
	190
[253]	191	! get column / layer conversion factor
[1308]	192	DO k = 1, nlayer
[787]	193	DO i = 1, ngrid
[253]	194	l2c(i,k)=(pplev(i,k)-pplev(i,k+1))/g
[135]	195	ENDDO
	196	ENDDO
	197
	198	! Vertical loop (from top to bottom)
	199	! We carry the rain with us and calculate that added by warm/cold precipitation
	200	! processes and that subtracted by evaporation at each level.
[1308]	201	DO k = nlayer, 1, -1
[135]	202
	203	IF (evap_prec) THEN ! note no rneb dependence!
[787]	204	DO i = 1, ngrid
[135]	205	IF (zrfl(i) .GT.0.) THEN
[253]	206
[728]	207	if(zt(i,k).gt.Tsat(i,k))then
[863]	208	!! treat the case where all liquid water should boil
	209	zqev=MIN((zt(i,k)-Tsat(i,k))RCPDl2c(i,k)/RLVTT/ptimestep,zrfl(i))
[728]	210	zrfl(i)=MAX(zrfl(i)-zqev,0.)
[863]	211	d_q(i,k)=zqev/l2c(i,k)*ptimestep
[728]	212	d_t(i,k) = - d_q(i,k) * RLVTT/RCPD
	213	else
[731]	214	zqev = MAX (0.0, (zqs(i,k)-q(i,k)))*l2c(i,k)/ptimestep !there was a bug here
[728]	215	zqevt= 2.0e-5*(1.0-q(i,k)/zqs(i,k)) & !default was 2.e-5
	216	sqrt(zrfl(i))l2c(i,k)/pplay(i,k)zt(i,k)R ! BC modif here
	217	zqevt = MAX (zqevt, 0.0)
	218	zqev = MIN (zqev, zqevt)
	219	zqev = MAX (zqev, 0.0)
	220	zrfln(i)= zrfl(i) - zqev
	221	zrfln(i)= max(zrfln(i),0.0)
[253]	222
[728]	223	d_q(i,k) = - (zrfln(i)-zrfl(i))/l2c(i,k)*ptimestep
	224	!d_t(i,k) = d_q(i,k) * RLVTT/RCPD!/(1.0+RVTMP2*q(i,k)) ! double BC modif here
	225	d_t(i,k) = - d_q(i,k) * RLVTT/RCPD ! was bugged!
	226	zrfl(i) = zrfln(i)
	227	end if
	228
[135]	229
[1283]	230	ENDIF ! of IF (zrfl(i) .GT.0.)
[135]	231	ENDDO
[1283]	232	ENDIF ! of IF (evap_prec)
[135]	233
[1283]	234	zoliq(1:ngrid) = 0.0
[135]	235
	236
[728]	237	if(precip_scheme.eq.1)then
[135]	238
[787]	239	DO i = 1, ngrid
[253]	240	ttemp = zt(i,k)
[650]	241	IF (ttemp .ge. T_h2O_ice_liq) THEN
[253]	242	lconvert=rainthreshold
	243	ELSEIF (ttemp .gt. t_crit) THEN
	244	lconvert=rainthreshold*(1.- t_crit/ttemp)
	245	lconvert=MAX(0.0,lconvert)
	246	ELSE
	247	lconvert=0.
	248	ENDIF
[135]	249
[253]	250
	251	IF (ql(i,k).gt.1.e-9) then
	252	zneb(i) = MAX(rneb(i,k), seuil_neb)
	253	IF ((ql(i,k)/zneb(i)).gt.lconvert)THEN ! precipitate!
[622]	254	d_ql(i,k) = -MAX((ql(i,k)-lconvert*zneb(i)),0.0)
[253]	255	zrfl(i) = zrfl(i) - d_ql(i,k)*l2c(i,k)/ptimestep
	256	ENDIF
	257	ENDIF
[135]	258	ENDDO
	259
[728]	260	elseif (precip_scheme.ge.2) then
	261
[787]	262	DO i = 1, ngrid
[135]	263	IF (rneb(i,k).GT.0.0) THEN
	264	zoliq(i) = ql(i,k)
[253]	265	zrho(i) = pplay(i,k) / ( zt(i,k) * R )
[135]	266	zdz(i) = (pplev(i,k)-pplev(i,k+1)) / (zrho(i)*g)
[650]	267	zfrac(i) = (zt(i,k)-T_h2O_ice_clouds) / (T_h2O_ice_liq-T_h2O_ice_clouds)
[135]	268	zfrac(i) = MAX(zfrac(i), 0.0)
	269	zfrac(i) = MIN(zfrac(i), 1.0)
	270	zneb(i) = MAX(rneb(i,k), seuil_neb)
	271	ENDIF
[731]	272	ENDDO
[135]	273
[731]	274	!recalculate liquid water particle radii
[1308]	275	call h2o_cloudrad(ngrid,nlayer,ql,reffh2oliq,reffh2oice)
[731]	276
[728]	277	SELECT CASE(precip_scheme)
	278	!precip scheme from Sundquist 78
	279	CASE(2)
	280
[135]	281	DO n = 1, ninter
[787]	282	DO i = 1, ngrid
[135]	283	IF (rneb(i,k).GT.0.0) THEN
[728]	284	! this is the ONLY place where zneb, precip_timescale and cloud_sat are used
[253]	285
[728]	286	zchau(i) = (ptimestep/(FLOAT(ninter)precip_timescale)) zoliq(i) &
	287	* (1.0-EXP(-(zoliq(i)/zneb(i)/cloud_sat)*2)) zfrac(i)
[135]	288	zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
	289	zfroi(i) = ptimestep/FLOAT(ninter)/zdz(i)*zoliq(i) &
[731]	290	fall2v(reffh2oice(i,k)) (1.0-zfrac(i)) ! zfroi behaves oddly...
[135]	291	ztot(i) = zchau(i) + zfroi(i)
	292
	293	IF (zneb(i).EQ.seuil_neb) ztot(i) = 0.0
	294	ztot(i) = MIN(MAX(ztot(i),0.0),zoliq(i))
	295	zoliq(i) = MAX(zoliq(i)-ztot(i), 0.0)
[253]	296
[135]	297	ENDIF
	298	ENDDO
	299	ENDDO
	300
[728]	301	!precip scheme modified from Sundquist 78 (in q**3)
	302	CASE(3)
	303
	304	DO n = 1, ninter
[787]	305	DO i = 1, ngrid
[728]	306	IF (rneb(i,k).GT.0.0) THEN
	307	! this is the ONLY place where zneb, precip_timescale and cloud_sat are used
	308
	309	zchau(i) = (ptimestep/(FLOAT(ninter)precip_timescalecloud_sat*2)) (zoliq(i)/zneb(i))**3
	310	zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
	311	zfroi(i) = ptimestep/FLOAT(ninter)/zdz(i)*zoliq(i) &
[731]	312	fall2v(reffh2oice(i,k)) (1.0-zfrac(i)) ! zfroi behaves oddly...
[728]	313	ztot(i) = zchau(i) + zfroi(i)
	314
	315	IF (zneb(i).EQ.seuil_neb) ztot(i) = 0.0
	316	ztot(i) = MIN(MAX(ztot(i),0.0),zoliq(i))
	317	zoliq(i) = MAX(zoliq(i)-ztot(i), 0.0)
	318
	319	ENDIF
	320	ENDDO
	321	ENDDO
	322
	323	!precip scheme modified from Boucher 95
	324	CASE(4)
	325
	326	DO n = 1, ninter
[787]	327	DO i = 1, ngrid
[728]	328	IF (rneb(i,k).GT.0.0) THEN
	329	! this is the ONLY place where zneb and c1 are used
	330
	331	zchau(i) = ptimestep/FLOAT(ninter) c1 zrho(i) &
	332	(zoliq(i)/zneb(i))2reffh2oliq(i,k)zneb(i) zfrac(i)
	333	zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
	334	zfroi(i) = ptimestep/FLOAT(ninter)/zdz(i)*zoliq(i) &
[731]	335	fall2v(reffh2oice(i,k)) (1.0-zfrac(i)) ! zfroi behaves oddly...
[728]	336	ztot(i) = zchau(i) + zfroi(i)
	337
	338	IF (zneb(i).EQ.seuil_neb) ztot(i) = 0.0
	339	ztot(i) = MIN(MAX(ztot(i),0.0),zoliq(i))
	340	zoliq(i) = MAX(zoliq(i)-ztot(i), 0.0)
	341
	342	ENDIF
	343	ENDDO
	344	ENDDO
	345
	346	END SELECT ! precip_scheme
	347
[135]	348	! Change in cloud density and surface H2O values
[787]	349	DO i = 1, ngrid
[135]	350	IF (rneb(i,k).GT.0.0) THEN
[253]	351	d_ql(i,k) = (zoliq(i) - ql(i,k))!/ptimestep
	352	zrfl(i) = zrfl(i)+ MAX(ql(i,k)-zoliq(i),0.0)*l2c(i,k)/ptimestep
[135]	353	ENDIF
	354	ENDDO
	355
	356
[728]	357	endif ! if precip_scheme=1
	358
[1308]	359	ENDDO ! of DO k = nlayer, 1, -1
[135]	360
	361	! Rain or snow on the ground
[787]	362	DO i = 1, ngrid
[253]	363	if(zrfl(i).lt.0.0)then
	364	print*,'Droplets of negative rain are falling...'
	365	call abort
	366	endif
[650]	367	IF (t(i,1) .LT. T_h2O_ice_liq) THEN
[135]	368	dqssnow(i) = zrfl(i)
[253]	369	dqsrain(i) = 0.0
[135]	370	ELSE
[253]	371	dqssnow(i) = 0.0
[135]	372	dqsrain(i) = zrfl(i) ! liquid water = ice for now
	373	ENDIF
	374	ENDDO
	375
	376	! now subroutine -----> GCM variables
[1283]	377	if (evap_prec) then
[1308]	378	dqrain(1:ngrid,1:nlayer,i_vap)=d_q(1:ngrid,1:nlayer)/ptimestep
	379	d_t(1:ngrid,1:nlayer)=d_t(1:ngrid,1:nlayer)/ptimestep
[1859]	380	do i=1,ngrid
	381	reevap_precip(i)=0.
	382	do k=1,nlayer
	383	reevap_precip(i)=reevap_precip(i)+dqrain(i,k,i_vap)*l2c(i,k)
	384	enddo
	385	enddo
[1283]	386	else
[1308]	387	dqrain(1:ngrid,1:nlayer,i_vap)=0.0
	388	d_t(1:ngrid,1:nlayer)=0.0
[1283]	389	endif
[1308]	390	dqrain(1:ngrid,1:nlayer,i_ice) = d_ql(1:ngrid,1:nlayer)/ptimestep
[135]	391
	392	end subroutine rain

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

Download in other formats:

Original Format