Context Navigation

rain.F90 @ 935

Last change on this file since 935 was 863, checked in by jleconte, 13 years ago

23/01/2013 == JL

Correction in largescale. a rneb factor was forgotten
Added some spectra in ave_stelpec
Corrected reevaporation in rain. Now conserve water better

File size: 13.8 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
[858]	68	logical,parameter :: evap_prec=.true. ! 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
[135]	145	firstcall = .false.
	146	ENDIF
	147
	148	! GCM -----> subroutine variables
	149	DO k = 1, nlayermx
[787]	150	DO i = 1, ngrid
[135]	151
[253]	152	zt(i,k) = t(i,k)+pdt(i,k)*ptimestep ! a big fat bug was here
	153	q(i,k) = pq(i,k,i_vap)+pdq(i,k,i_vap)*ptimestep
	154	ql(i,k) = pq(i,k,i_ice)+pdq(i,k,i_ice)*ptimestep
[135]	155
[253]	156	!q(i,k) = pq(i,k,i_vap)!+pdq(i,k,i_vap)
	157	!ql(i,k) = pq(i,k,i_ice)!+pdq(i,k,i_ice)
	158
[135]	159	if(q(i,k).lt.0.)then ! if this is not done, we don't conserve water
	160	q(i,k)=0.
	161	endif
	162	if(ql(i,k).lt.0.)then
	163	ql(i,k)=0.
	164	endif
	165
	166	ENDDO
	167	ENDDO
	168
	169	! Initialise the outputs
	170	DO k = 1, nlayermx
[787]	171	DO i = 1, ngrid
[135]	172	d_t(i,k) = 0.0
	173	d_q(i,k) = 0.0
	174	d_ql(i,k) = 0.0
	175	ENDDO
	176	ENDDO
[787]	177	DO i = 1, ngrid
[135]	178	zrfl(i) = 0.0
	179	zrfln(i) = 0.0
	180	ENDDO
	181
	182	! calculate saturation mixing ratio
	183	DO k = 1, nlayermx
[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
[135]	194	DO k = 1, nlayermx
[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
	201	! Vertical loop (from top to bottom)
	202	! We carry the rain with us and calculate that added by warm/cold precipitation
	203	! processes and that subtracted by evaporation at each level.
	204	DO 9999 k = nlayermx, 1, -1
	205
	206	IF (evap_prec) THEN ! note no rneb dependence!
[787]	207	DO i = 1, ngrid
[135]	208	IF (zrfl(i) .GT.0.) THEN
[253]	209
[728]	210	if(zt(i,k).gt.Tsat(i,k))then
[863]	211	!! treat the case where all liquid water should boil
	212	zqev=MIN((zt(i,k)-Tsat(i,k))RCPDl2c(i,k)/RLVTT/ptimestep,zrfl(i))
[728]	213	zrfl(i)=MAX(zrfl(i)-zqev,0.)
[863]	214	d_q(i,k)=zqev/l2c(i,k)*ptimestep
[728]	215	d_t(i,k) = - d_q(i,k) * RLVTT/RCPD
	216	else
[731]	217	zqev = MAX (0.0, (zqs(i,k)-q(i,k)))*l2c(i,k)/ptimestep !there was a bug here
[728]	218	zqevt= 2.0e-5*(1.0-q(i,k)/zqs(i,k)) & !default was 2.e-5
	219	sqrt(zrfl(i))l2c(i,k)/pplay(i,k)zt(i,k)R ! BC modif here
	220	zqevt = MAX (zqevt, 0.0)
	221	zqev = MIN (zqev, zqevt)
	222	zqev = MAX (zqev, 0.0)
	223	zrfln(i)= zrfl(i) - zqev
	224	zrfln(i)= max(zrfln(i),0.0)
[253]	225
[728]	226	d_q(i,k) = - (zrfln(i)-zrfl(i))/l2c(i,k)*ptimestep
	227	!d_t(i,k) = d_q(i,k) * RLVTT/RCPD!/(1.0+RVTMP2*q(i,k)) ! double BC modif here
	228	d_t(i,k) = - d_q(i,k) * RLVTT/RCPD ! was bugged!
	229	zrfl(i) = zrfln(i)
	230	end if
	231
[135]	232
	233	ENDIF
	234	ENDDO
	235	ENDIF
	236
[787]	237	DO i = 1, ngrid
[135]	238	zoliq(i) = 0.0
	239	ENDDO
	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)
[253]	260	zrfl(i) = zrfl(i) - d_ql(i,k)*l2c(i,k)/ptimestep
	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
[787]	280	call h2o_cloudrad(ngrid,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
	357	zrfl(i) = zrfl(i)+ MAX(ql(i,k)-zoliq(i),0.0)*l2c(i,k)/ptimestep
[135]	358	ENDIF
	359	ENDDO
	360
	361
[728]	362	endif ! if precip_scheme=1
	363
[135]	364	9999 continue
	365
	366	! Rain or snow on the ground
[787]	367	DO i = 1, ngrid
[253]	368	if(zrfl(i).lt.0.0)then
	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
[135]	373	dqssnow(i) = zrfl(i)
[253]	374	dqsrain(i) = 0.0
[135]	375	ELSE
[253]	376	dqssnow(i) = 0.0
[135]	377	dqsrain(i) = zrfl(i) ! liquid water = ice for now
	378	ENDIF
	379	ENDDO
	380
	381	! now subroutine -----> GCM variables
	382	DO k = 1, nlayermx
[787]	383	DO i = 1, ngrid
[135]	384
	385	if(evap_prec)then
	386	dqrain(i,k,i_vap) = d_q(i,k)/ptimestep
	387	d_t(i,k) = d_t(i,k)/ptimestep
	388	else
	389	dqrain(i,k,i_vap) = 0.0
	390	d_t(i,k) = 0.0
	391	endif
[253]	392	dqrain(i,k,i_ice) = d_ql(i,k)/ptimestep
[135]	393
	394	ENDDO
	395	ENDDO
	396
	397	RETURN
	398	end subroutine rain

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

Download in other formats:

Original Format