Context Navigation

rain.F90 @ 735

Last change on this file since 735 was 731, checked in by jleconte, 13 years ago
Corrected precipitation evaporation scheme + snow fall (JL+BC) Commented some unecessary writediagfi calls that where used for test.
File size: 13.8 KB

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

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

Download in other formats:

Original Format