Context Navigation

rain.F90 @ 1476

Last change on this file since 1476 was 1397, checked in by milmd, 10 years ago
In LMDZ.GENERIC replacement of all phystd .h files by module files.
File size: 14.1 KB

Rev	Line
[1308]	1	subroutine rain(ngrid,nlayer,nq,ptimestep,pplev,pplay,t,pdt,pq,pdq,d_t,dqrain,dqsrain,dqssnow,rneb)
[135]	2
[253]	3
[1315]	4	! to use 'getin'
	5	! use ioipsl_getincom
	6	use ioipsl_getincom_p
[728]	7	use watercommon_h, only: T_h2O_ice_liq,T_h2O_ice_clouds, RLVTT, RCPD, RCPV, RV, RVTMP2,Psat_water,Tsat_water,rhowater
	8	use radii_mod, only: h2o_cloudrad
[1283]	9	USE tracer_h, only: igcm_h2o_vap, igcm_h2o_ice
[1384]	10	use comcstfi_mod, only: g, r
[135]	11	implicit none
	12
	13	!==================================================================
	14	!
	15	! Purpose
	16	! -------
	17	! Calculates H2O precipitation using simplified microphysics.
	18	!
	19	! Authors
	20	! -------
	21	! Adapted from the LMDTERRE code by R. Wordsworth (2009)
[728]	22	! Added rain vaporization in case of T>Tsat
[135]	23	! Original author Z. X. Li (1993)
	24	!
	25	!==================================================================
	26
[858]	27	! Arguments
[1308]	28	integer,intent(in) :: ngrid ! number of atmospheric columns
	29	integer,intent(in) :: nlayer ! number of atmospheric layers
[858]	30	integer,intent(in) :: nq ! number of tracers
	31	real,intent(in) :: ptimestep ! time interval
[1308]	32	real,intent(in) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa)
	33	real,intent(in) :: pplay(ngrid,nlayer) ! mid-layer pressure (Pa)
	34	real,intent(in) :: t(ngrid,nlayer) ! input temperature (K)
	35	real,intent(in) :: pdt(ngrid,nlayer) ! input tendency on temperature (K/s)
	36	real,intent(in) :: pq(ngrid,nlayer,nq) ! tracers (kg/kg)
	37	real,intent(in) :: pdq(ngrid,nlayer,nq) ! input tendency on tracers
	38	real,intent(out) :: d_t(ngrid,nlayer) ! temperature tendency (K/s)
	39	real,intent(out) :: dqrain(ngrid,nlayer,nq) ! tendency of H2O precipitation (kg/kg.s-1)
[858]	40	real,intent(out) :: dqsrain(ngrid) ! rain flux at the surface (kg.m-2.s-1)
	41	real,intent(out) :: dqssnow(ngrid) ! snow flux at the surface (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
[787]	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
[1308]	83	real reffh2oliq(ngrid,nlayer),reffh2oice(ngrid,nlayer)
[728]	84
	85	real ttemp, ptemp, psat_tmp
[1308]	86	real tnext(ngrid,nlayer)
[135]	87
[1308]	88	real l2c(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
[1283]	180	zrfl(1:ngrid) = 0.0
	181	zrfln(1:ngrid) = 0.0
[135]	182
	183	! calculate saturation mixing ratio
[1308]	184	DO k = 1, nlayer
[787]	185	DO i = 1, ngrid
[253]	186	ttemp = zt(i,k)
[135]	187	ptemp = pplay(i,k)
[728]	188	! call watersat(ttemp,ptemp,zqs(i,k))
	189	call Psat_water(ttemp,ptemp,psat_tmp,zqs(i,k))
	190	call Tsat_water(ptemp,Tsat(i,k))
[135]	191	ENDDO
	192	ENDDO
	193
[253]	194	! get column / layer conversion factor
[1308]	195	DO k = 1, nlayer
[787]	196	DO i = 1, ngrid
[253]	197	l2c(i,k)=(pplev(i,k)-pplev(i,k+1))/g
[135]	198	ENDDO
	199	ENDDO
	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.
[1308]	204	DO k = nlayer, 1, -1
[135]	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
[1283]	233	ENDIF ! of IF (zrfl(i) .GT.0.)
[135]	234	ENDDO
[1283]	235	ENDIF ! of IF (evap_prec)
[135]	236
[1283]	237	zoliq(1:ngrid) = 0.0
[135]	238
	239
[728]	240	if(precip_scheme.eq.1)then
[135]	241
[787]	242	DO i = 1, ngrid
[253]	243	ttemp = zt(i,k)
[650]	244	IF (ttemp .ge. T_h2O_ice_liq) THEN
[253]	245	lconvert=rainthreshold
	246	ELSEIF (ttemp .gt. t_crit) THEN
	247	lconvert=rainthreshold*(1.- t_crit/ttemp)
	248	lconvert=MAX(0.0,lconvert)
	249	ELSE
	250	lconvert=0.
	251	ENDIF
[135]	252
[253]	253
	254	IF (ql(i,k).gt.1.e-9) then
	255	zneb(i) = MAX(rneb(i,k), seuil_neb)
	256	IF ((ql(i,k)/zneb(i)).gt.lconvert)THEN ! precipitate!
[622]	257	d_ql(i,k) = -MAX((ql(i,k)-lconvert*zneb(i)),0.0)
[253]	258	zrfl(i) = zrfl(i) - d_ql(i,k)*l2c(i,k)/ptimestep
	259	ENDIF
	260	ENDIF
[135]	261	ENDDO
	262
[728]	263	elseif (precip_scheme.ge.2) then
	264
[787]	265	DO i = 1, ngrid
[135]	266	IF (rneb(i,k).GT.0.0) THEN
	267	zoliq(i) = ql(i,k)
[253]	268	zrho(i) = pplay(i,k) / ( zt(i,k) * R )
[135]	269	zdz(i) = (pplev(i,k)-pplev(i,k+1)) / (zrho(i)*g)
[650]	270	zfrac(i) = (zt(i,k)-T_h2O_ice_clouds) / (T_h2O_ice_liq-T_h2O_ice_clouds)
[135]	271	zfrac(i) = MAX(zfrac(i), 0.0)
	272	zfrac(i) = MIN(zfrac(i), 1.0)
	273	zneb(i) = MAX(rneb(i,k), seuil_neb)
	274	ENDIF
[731]	275	ENDDO
[135]	276
[731]	277	!recalculate liquid water particle radii
[1308]	278	call h2o_cloudrad(ngrid,nlayer,ql,reffh2oliq,reffh2oice)
[731]	279
[728]	280	SELECT CASE(precip_scheme)
	281	!precip scheme from Sundquist 78
	282	CASE(2)
	283
[135]	284	DO n = 1, ninter
[787]	285	DO i = 1, ngrid
[135]	286	IF (rneb(i,k).GT.0.0) THEN
[728]	287	! this is the ONLY place where zneb, precip_timescale and cloud_sat are used
[253]	288
[728]	289	zchau(i) = (ptimestep/(FLOAT(ninter)precip_timescale)) zoliq(i) &
	290	* (1.0-EXP(-(zoliq(i)/zneb(i)/cloud_sat)*2)) zfrac(i)
[135]	291	zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
	292	zfroi(i) = ptimestep/FLOAT(ninter)/zdz(i)*zoliq(i) &
[731]	293	fall2v(reffh2oice(i,k)) (1.0-zfrac(i)) ! zfroi behaves oddly...
[135]	294	ztot(i) = zchau(i) + zfroi(i)
	295
	296	IF (zneb(i).EQ.seuil_neb) ztot(i) = 0.0
	297	ztot(i) = MIN(MAX(ztot(i),0.0),zoliq(i))
	298	zoliq(i) = MAX(zoliq(i)-ztot(i), 0.0)
[253]	299
[135]	300	ENDIF
	301	ENDDO
	302	ENDDO
	303
[728]	304	!precip scheme modified from Sundquist 78 (in q**3)
	305	CASE(3)
	306
	307	DO n = 1, ninter
[787]	308	DO i = 1, ngrid
[728]	309	IF (rneb(i,k).GT.0.0) THEN
	310	! this is the ONLY place where zneb, precip_timescale and cloud_sat are used
	311
	312	zchau(i) = (ptimestep/(FLOAT(ninter)precip_timescalecloud_sat*2)) (zoliq(i)/zneb(i))**3
	313	zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
	314	zfroi(i) = ptimestep/FLOAT(ninter)/zdz(i)*zoliq(i) &
[731]	315	fall2v(reffh2oice(i,k)) (1.0-zfrac(i)) ! zfroi behaves oddly...
[728]	316	ztot(i) = zchau(i) + zfroi(i)
	317
	318	IF (zneb(i).EQ.seuil_neb) ztot(i) = 0.0
	319	ztot(i) = MIN(MAX(ztot(i),0.0),zoliq(i))
	320	zoliq(i) = MAX(zoliq(i)-ztot(i), 0.0)
	321
	322	ENDIF
	323	ENDDO
	324	ENDDO
	325
	326	!precip scheme modified from Boucher 95
	327	CASE(4)
	328
	329	DO n = 1, ninter
[787]	330	DO i = 1, ngrid
[728]	331	IF (rneb(i,k).GT.0.0) THEN
	332	! this is the ONLY place where zneb and c1 are used
	333
	334	zchau(i) = ptimestep/FLOAT(ninter) c1 zrho(i) &
	335	(zoliq(i)/zneb(i))2reffh2oliq(i,k)zneb(i) zfrac(i)
	336	zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
	337	zfroi(i) = ptimestep/FLOAT(ninter)/zdz(i)*zoliq(i) &
[731]	338	fall2v(reffh2oice(i,k)) (1.0-zfrac(i)) ! zfroi behaves oddly...
[728]	339	ztot(i) = zchau(i) + zfroi(i)
	340
	341	IF (zneb(i).EQ.seuil_neb) ztot(i) = 0.0
	342	ztot(i) = MIN(MAX(ztot(i),0.0),zoliq(i))
	343	zoliq(i) = MAX(zoliq(i)-ztot(i), 0.0)
	344
	345	ENDIF
	346	ENDDO
	347	ENDDO
	348
	349	END SELECT ! precip_scheme
	350
[135]	351	! Change in cloud density and surface H2O values
[787]	352	DO i = 1, ngrid
[135]	353	IF (rneb(i,k).GT.0.0) THEN
[253]	354	d_ql(i,k) = (zoliq(i) - ql(i,k))!/ptimestep
	355	zrfl(i) = zrfl(i)+ MAX(ql(i,k)-zoliq(i),0.0)*l2c(i,k)/ptimestep
[135]	356	ENDIF
	357	ENDDO
	358
	359
[728]	360	endif ! if precip_scheme=1
	361
[1308]	362	ENDDO ! of DO k = nlayer, 1, -1
[135]	363
	364	! Rain or snow on the ground
[787]	365	DO i = 1, ngrid
[253]	366	if(zrfl(i).lt.0.0)then
	367	print*,'Droplets of negative rain are falling...'
	368	call abort
	369	endif
[650]	370	IF (t(i,1) .LT. T_h2O_ice_liq) THEN
[135]	371	dqssnow(i) = zrfl(i)
[253]	372	dqsrain(i) = 0.0
[135]	373	ELSE
[253]	374	dqssnow(i) = 0.0
[135]	375	dqsrain(i) = zrfl(i) ! liquid water = ice for now
	376	ENDIF
	377	ENDDO
	378
	379	! now subroutine -----> GCM variables
[1283]	380	if (evap_prec) then
[1308]	381	dqrain(1:ngrid,1:nlayer,i_vap)=d_q(1:ngrid,1:nlayer)/ptimestep
	382	d_t(1:ngrid,1:nlayer)=d_t(1:ngrid,1:nlayer)/ptimestep
[1283]	383	else
[1308]	384	dqrain(1:ngrid,1:nlayer,i_vap)=0.0
	385	d_t(1:ngrid,1:nlayer)=0.0
[1283]	386	endif
[1308]	387	dqrain(1:ngrid,1:nlayer,i_ice) = d_ql(1:ngrid,1:nlayer)/ptimestep
[135]	388
	389	end subroutine rain

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

Download in other formats:

Original Format