Context Navigation

rain.F90 @ 724

Last change on this file since 724 was 650, checked in by jleconte, 13 years ago

Corrected the temperature used to differentiate sublimation and evaporation in watersat_grad
Minor name changes in watercommon
Better physical parametrization of the effective radius of liquid and icy water cloud particles in callcorrk

(for radfixed=true)

Added consistency check in inifis
Moved 1d water initialization from physiqu to rcm1d

File size: 10.0 KB

Rev	Line
[135]	1	subroutine rain(ptimestep,pplev,pplay,t,pdt,pq,pdq,d_t,dqrain,dqsrain,dqssnow,rneb)
	2
[253]	3
[650]	4	use watercommon_h, only: T_h2O_ice_liq,T_h2O_ice_clouds, RLVTT, RCPD, RCPV, RV, RVTMP2
[135]	5
	6	implicit none
	7
	8	!==================================================================
	9	!
	10	! Purpose
	11	! -------
	12	! Calculates H2O precipitation using simplified microphysics.
	13	!
	14	! Authors
	15	! -------
	16	! Adapted from the LMDTERRE code by R. Wordsworth (2009)
	17	! Original author Z. X. Li (1993)
	18	!
	19	!==================================================================
	20
	21	#include "dimensions.h"
	22	#include "dimphys.h"
	23	#include "tracer.h"
	24	#include "comcstfi.h"
	25	#include "callkeys.h"
	26
	27	! Pre-arguments (for universal model)
	28	real pq(ngridmx,nlayermx,nqmx) ! tracer (kg/kg)
	29	real qsurf(ngridmx,nqmx) ! tracer at the surface (kg.m-2)
	30	REAL pdt(ngridmx,nlayermx),pdq(ngridmx,nlayermx,nqmx)
	31
	32	real dqrain(ngridmx,nlayermx,nqmx) ! tendency of H2O precipitation (kg/kg.s-1)
	33	real dqsrain(ngridmx) ! rain flux at the surface (kg.m-2.s-1)
	34	real dqssnow(ngridmx) ! snow flux at the surface (kg.m-2.s-1)
	35	REAL d_t(ngridmx,nlayermx) ! temperature increment
	36
	37	! Arguments
	38	REAL ptimestep ! time interval
	39	REAL pplev(ngridmx,nlayermx+1) ! inter-layer pressure
	40	REAL pplay(ngridmx,nlayermx) ! mid-layer pressure
[253]	41	REAL t(ngridmx,nlayermx) ! input temperature (K)
	42	REAL zt(ngridmx,nlayermx) ! working temperature (K)
[135]	43	REAL ql(ngridmx,nlayermx) ! liquid water (Kg/Kg)
	44	REAL q(ngridmx,nlayermx) ! specific humidity (Kg/Kg)
	45	REAL rneb(ngridmx,nlayermx) ! cloud fraction
	46	REAL d_q(ngridmx,nlayermx) ! water vapor increment
	47	REAL d_ql(ngridmx,nlayermx) ! liquid water / ice increment
	48
	49	! Subroutine options
	50	REAL seuil_neb ! Nebulosity threshold
	51	PARAMETER (seuil_neb=0.001)
	52
[622]	53	REAL,PARAMETER :: cl=2.0e-4 ! Precipitation threshold
	54	REAL,PARAMETER :: ct=1./1800. ! Precipitation rate
[135]	55
	56	INTEGER ninter
	57	PARAMETER (ninter=5)
	58
	59	logical simple ! Use very simple Emanuel scheme?
[622]	60	parameter(simple=.true.)
[135]	61
	62	logical evap_prec ! Does the rain evaporate?
	63	parameter(evap_prec=.true.)
	64
	65	! for simple scheme
	66	real t_crit
	67	PARAMETER (t_crit=218.0)
	68	real lconvert
	69
	70	! for precipitation evaporation (old scheme)
	71	real eeff1
	72	real eeff2
	73	! parameter (eeff1=0.95)
	74	! parameter (eeff2=0.98)
	75	parameter (eeff1=0.5)
	76	parameter (eeff2=0.8)
	77
	78	! Local variables
	79	INTEGER i, k, n
	80	REAL zqs(ngridmx,nlayermx), zdelta, zcor
	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
	87	real ttemp, ptemp
	88	real tnext(ngridmx,nlayermx)
	89
[253]	90	real l2c(ngridmx,nlayermx)
	91	real dWtot
[135]	92
[253]	93
[135]	94	! Indices of water vapour and water ice tracers
	95	INTEGER, SAVE :: i_vap=0 ! water vapour
	96	INTEGER, SAVE :: i_ice=0 ! water ice
	97
	98	LOGICAL firstcall
	99	SAVE firstcall
	100
	101	! Online functions
	102	REAL fallv, zzz ! falling speed of ice crystals
	103	fallv (zzz) = 3.29 * ((zzz)**0.16)
	104
	105	DATA firstcall /.true./
	106
	107	IF (firstcall) THEN
	108
	109	i_vap=igcm_h2o_vap
	110	i_ice=igcm_h2o_ice
	111
	112	write(,) "rain: i_ice=",i_ice
	113	write(,) " i_vap=",i_vap
	114
	115	PRINT*, 'in rain.F, ninter=', ninter
	116	PRINT*, 'in rain.F, evap_prec=', evap_prec
	117
	118	firstcall = .false.
	119	ENDIF
	120
	121	! GCM -----> subroutine variables
	122	DO k = 1, nlayermx
	123	DO i = 1, ngridmx
	124
[253]	125	zt(i,k) = t(i,k)+pdt(i,k)*ptimestep ! a big fat bug was here
	126	q(i,k) = pq(i,k,i_vap)+pdq(i,k,i_vap)*ptimestep
	127	ql(i,k) = pq(i,k,i_ice)+pdq(i,k,i_ice)*ptimestep
[135]	128
[253]	129	!q(i,k) = pq(i,k,i_vap)!+pdq(i,k,i_vap)
	130	!ql(i,k) = pq(i,k,i_ice)!+pdq(i,k,i_ice)
	131
[135]	132	if(q(i,k).lt.0.)then ! if this is not done, we don't conserve water
	133	q(i,k)=0.
	134	endif
	135	if(ql(i,k).lt.0.)then
	136	ql(i,k)=0.
	137	endif
	138
	139	ENDDO
	140	ENDDO
	141
	142	! Initialise the outputs
	143	DO k = 1, nlayermx
	144	DO i = 1, ngridmx
	145	d_t(i,k) = 0.0
	146	d_q(i,k) = 0.0
	147	d_ql(i,k) = 0.0
	148	ENDDO
	149	ENDDO
	150	DO i = 1, ngridmx
	151	zrfl(i) = 0.0
	152	zrfln(i) = 0.0
	153	ENDDO
	154
	155	! calculate saturation mixing ratio
	156	DO k = 1, nlayermx
	157	DO i = 1, ngridmx
[253]	158	ttemp = zt(i,k)
[135]	159	ptemp = pplay(i,k)
[253]	160	call watersat(ttemp,ptemp,zqs(i,k))
[135]	161	ENDDO
	162	ENDDO
	163
[253]	164	! get column / layer conversion factor
[135]	165	DO k = 1, nlayermx
	166	DO i = 1, ngridmx
[253]	167	!l2c(i,k)=(pplev(i,k)-pplev(i,k+1))/(g*ptimestep)
	168	l2c(i,k)=(pplev(i,k)-pplev(i,k+1))/g
[135]	169	ENDDO
	170	ENDDO
	171
	172
	173	! Vertical loop (from top to bottom)
	174	! We carry the rain with us and calculate that added by warm/cold precipitation
	175	! processes and that subtracted by evaporation at each level.
	176	DO 9999 k = nlayermx, 1, -1
	177
	178	IF (evap_prec) THEN ! note no rneb dependence!
	179	DO i = 1, ngridmx
	180	IF (zrfl(i) .GT.0.) THEN
[253]	181
	182	zqev = MAX (0.0, (zqs(i,k)-q(i,k)))/ptimestep! BC modif here
	183	zqevt = 2.0e-5*(1.0-q(i,k)/zqs(i,k)) &
	184	sqrt(zrfl(i))l2c(i,k)/pplay(i,k)zt(i,k)R ! BC modif here
	185	zqevt = MAX (zqevt, 0.0)
	186	zqev = MIN (zqev, zqevt)
	187	zqev = MAX (zqev, 0.0)
	188	zrfln(i) = zrfl(i) - zqev
	189	zrfln(i) = max(zrfln(i),0.0)
	190
[135]	191	!zqev = MAX (0.0, (zqs(i,k)-q(i,k))*eeff1 )
	192	!zqevt = (zrfl(i)/l2c(i,k))*eeff2
	193	!zqev = MIN (zqev, zqevt)
	194	!zrfln(i) = zrfl(i) - zqev*l2c(i,k)
[253]	195	!zrfln(i) = zrfl(i) - 1.5e-5(1.0-q(i,k)/zqs(i,k))sqrt(zrfl(i))
	196	!zrfln(i) = min(zrfln(i),zrfl(i))
[135]	197	! this is what is actually written in the manual
	198
[253]	199	d_q(i,k) = - (zrfln(i)-zrfl(i))/l2c(i,k)*ptimestep
	200	!d_t(i,k) = d_q(i,k) * RLVTT/RCPD!/(1.0+RVTMP2*q(i,k)) ! double BC modif here
	201	d_t(i,k) = - d_q(i,k) * RLVTT/RCPD ! was bugged!
	202
[135]	203	zrfl(i) = zrfln(i)
	204	ENDIF
	205	ENDDO
	206	ENDIF
	207
	208	DO i = 1, ngridmx
	209	zoliq(i) = 0.0
	210	ENDDO
	211
	212
	213	if(simple)then
	214
	215	DO i = 1, ngridmx
[253]	216	ttemp = zt(i,k)
[650]	217	IF (ttemp .ge. T_h2O_ice_liq) THEN
[253]	218	lconvert=rainthreshold
	219	ELSEIF (ttemp .gt. t_crit) THEN
	220	lconvert=rainthreshold*(1.- t_crit/ttemp)
	221	lconvert=MAX(0.0,lconvert)
	222	ELSE
	223	lconvert=0.
	224	ENDIF
[135]	225
[253]	226
	227	IF (ql(i,k).gt.1.e-9) then
	228	zneb(i) = MAX(rneb(i,k), seuil_neb)
	229	IF ((ql(i,k)/zneb(i)).gt.lconvert)THEN ! precipitate!
[622]	230	d_ql(i,k) = -MAX((ql(i,k)-lconvert*zneb(i)),0.0)
[253]	231	zrfl(i) = zrfl(i) - d_ql(i,k)*l2c(i,k)/ptimestep
	232	ENDIF
	233	ENDIF
[135]	234	ENDDO
	235
	236	else
	237
	238	DO i = 1, ngridmx
	239	IF (rneb(i,k).GT.0.0) THEN
	240	zoliq(i) = ql(i,k)
[253]	241	zrho(i) = pplay(i,k) / ( zt(i,k) * R )
[135]	242	zdz(i) = (pplev(i,k)-pplev(i,k+1)) / (zrho(i)*g)
[650]	243	zfrac(i) = (zt(i,k)-T_h2O_ice_clouds) / (T_h2O_ice_liq-T_h2O_ice_clouds)
[135]	244	zfrac(i) = MAX(zfrac(i), 0.0)
	245	zfrac(i) = MIN(zfrac(i), 1.0)
	246	zneb(i) = MAX(rneb(i,k), seuil_neb)
	247	ENDIF
	248	ENDDO
	249
	250	DO n = 1, ninter
	251	DO i = 1, ngridmx
	252	IF (rneb(i,k).GT.0.0) THEN
[622]	253	zchau(i) = (ctptimestep/FLOAT(ninter)) zoliq(i) &
[135]	254	* (1.0-EXP(-(zoliq(i)/zneb(i)/cl)*2)) zfrac(i)
[253]	255
[135]	256	! this is the ONLY place where zneb, ct and cl are used
	257	zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
	258	zfroi(i) = ptimestep/FLOAT(ninter)/zdz(i)*zoliq(i) &
	259	fallv(zrhol(i)) (1.0-zfrac(i)) ! zfroi behaves oddly...
	260	! * 0.1 * (1.0-zfrac(i))
	261	ztot(i) = zchau(i) + zfroi(i)
	262
	263	IF (zneb(i).EQ.seuil_neb) ztot(i) = 0.0
	264	ztot(i) = MIN(MAX(ztot(i),0.0),zoliq(i))
	265	zoliq(i) = MAX(zoliq(i)-ztot(i), 0.0)
[253]	266
[135]	267	ENDIF
	268	ENDDO
	269	ENDDO
	270
	271	! Change in cloud density and surface H2O values
	272	DO i = 1, ngridmx
	273	IF (rneb(i,k).GT.0.0) THEN
[253]	274	d_ql(i,k) = (zoliq(i) - ql(i,k))!/ptimestep
	275	zrfl(i) = zrfl(i)+ MAX(ql(i,k)-zoliq(i),0.0)*l2c(i,k)/ptimestep
[135]	276	ENDIF
	277	ENDDO
	278
	279	endif ! if simple
	280
	281	9999 continue
	282
	283	! Rain or snow on the ground
	284	DO i = 1, ngridmx
[253]	285	if(zrfl(i).lt.0.0)then
	286	print*,'Droplets of negative rain are falling...'
	287	call abort
	288	endif
[650]	289	IF (t(i,1) .LT. T_h2O_ice_liq) THEN
[135]	290	dqssnow(i) = zrfl(i)
[253]	291	dqsrain(i) = 0.0
[135]	292	ELSE
[253]	293	dqssnow(i) = 0.0
[135]	294	dqsrain(i) = zrfl(i) ! liquid water = ice for now
	295	ENDIF
	296	ENDDO
	297
	298	! now subroutine -----> GCM variables
	299	DO k = 1, nlayermx
	300	DO i = 1, ngridmx
	301
	302	if(evap_prec)then
	303	dqrain(i,k,i_vap) = d_q(i,k)/ptimestep
	304	d_t(i,k) = d_t(i,k)/ptimestep
	305	else
	306	dqrain(i,k,i_vap) = 0.0
	307	d_t(i,k) = 0.0
	308	endif
[253]	309	dqrain(i,k,i_ice) = d_ql(i,k)/ptimestep
[135]	310
	311	ENDDO
	312	ENDDO
	313
	314	RETURN
	315	end subroutine rain

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

Download in other formats:

Original Format