Context Navigation

rain.F90 @ 648

Last change on this file since 648 was 622, checked in by jleconte, 13 years ago
Added consistency checks for calculations including water and global1d+diurnal. Corrected small bugs in precipitation scheme
File size: 10.1 KB

Rev	Line
[135]	1	subroutine rain(ptimestep,pplev,pplay,t,pdt,pq,pdq,d_t,dqrain,dqsrain,dqssnow,rneb)
	2
[253]	3
[135]	4	use watercommon_h, only: To, RLVTT, RCPD, RCPV, RV, RVTMP2
	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 ztglace
	85	REAL zdz(ngridmx),zrho(ngridmx),ztot(ngridmx), zrhol(ngridmx)
	86	REAL zchau(ngridmx),zfroi(ngridmx),zfrac(ngridmx),zneb(ngridmx)
	87
	88	real ttemp, ptemp
	89	real tnext(ngridmx,nlayermx)
	90
[253]	91	real l2c(ngridmx,nlayermx)
	92	real dWtot
[135]	93
[253]	94
[135]	95	! Indices of water vapour and water ice tracers
	96	INTEGER, SAVE :: i_vap=0 ! water vapour
	97	INTEGER, SAVE :: i_ice=0 ! water ice
	98
	99	LOGICAL firstcall
	100	SAVE firstcall
	101
	102	! Online functions
	103	REAL fallv, zzz ! falling speed of ice crystals
	104	fallv (zzz) = 3.29 * ((zzz)**0.16)
	105
	106	DATA firstcall /.true./
	107
	108	IF (firstcall) THEN
	109
	110	i_vap=igcm_h2o_vap
	111	i_ice=igcm_h2o_ice
	112
	113	write(,) "rain: i_ice=",i_ice
	114	write(,) " i_vap=",i_vap
	115
	116	PRINT*, 'in rain.F, ninter=', ninter
	117	PRINT*, 'in rain.F, evap_prec=', evap_prec
	118
	119	firstcall = .false.
	120	ENDIF
	121
	122	! GCM -----> subroutine variables
	123	DO k = 1, nlayermx
	124	DO i = 1, ngridmx
	125
[253]	126	zt(i,k) = t(i,k)+pdt(i,k)*ptimestep ! a big fat bug was here
	127	q(i,k) = pq(i,k,i_vap)+pdq(i,k,i_vap)*ptimestep
	128	ql(i,k) = pq(i,k,i_ice)+pdq(i,k,i_ice)*ptimestep
[135]	129
[253]	130	!q(i,k) = pq(i,k,i_vap)!+pdq(i,k,i_vap)
	131	!ql(i,k) = pq(i,k,i_ice)!+pdq(i,k,i_ice)
	132
[135]	133	if(q(i,k).lt.0.)then ! if this is not done, we don't conserve water
	134	q(i,k)=0.
	135	endif
	136	if(ql(i,k).lt.0.)then
	137	ql(i,k)=0.
	138	endif
	139
	140	ENDDO
	141	ENDDO
	142
	143	! Determine the cold clouds by their temperature
	144	ztglace = To - 15.0
	145
	146	! Initialise the outputs
	147	DO k = 1, nlayermx
	148	DO i = 1, ngridmx
	149	d_t(i,k) = 0.0
	150	d_q(i,k) = 0.0
	151	d_ql(i,k) = 0.0
	152	ENDDO
	153	ENDDO
	154	DO i = 1, ngridmx
	155	zrfl(i) = 0.0
	156	zrfln(i) = 0.0
	157	ENDDO
	158
	159	! calculate saturation mixing ratio
	160	DO k = 1, nlayermx
	161	DO i = 1, ngridmx
[253]	162	ttemp = zt(i,k)
[135]	163	ptemp = pplay(i,k)
[253]	164	call watersat(ttemp,ptemp,zqs(i,k))
[135]	165	ENDDO
	166	ENDDO
	167
[253]	168	! get column / layer conversion factor
[135]	169	DO k = 1, nlayermx
	170	DO i = 1, ngridmx
[253]	171	!l2c(i,k)=(pplev(i,k)-pplev(i,k+1))/(g*ptimestep)
	172	l2c(i,k)=(pplev(i,k)-pplev(i,k+1))/g
[135]	173	ENDDO
	174	ENDDO
	175
	176
	177	! Vertical loop (from top to bottom)
	178	! We carry the rain with us and calculate that added by warm/cold precipitation
	179	! processes and that subtracted by evaporation at each level.
	180	DO 9999 k = nlayermx, 1, -1
	181
	182	IF (evap_prec) THEN ! note no rneb dependence!
	183	DO i = 1, ngridmx
	184	IF (zrfl(i) .GT.0.) THEN
[253]	185
	186	zqev = MAX (0.0, (zqs(i,k)-q(i,k)))/ptimestep! BC modif here
	187	zqevt = 2.0e-5*(1.0-q(i,k)/zqs(i,k)) &
	188	sqrt(zrfl(i))l2c(i,k)/pplay(i,k)zt(i,k)R ! BC modif here
	189	zqevt = MAX (zqevt, 0.0)
	190	zqev = MIN (zqev, zqevt)
	191	zqev = MAX (zqev, 0.0)
	192	zrfln(i) = zrfl(i) - zqev
	193	zrfln(i) = max(zrfln(i),0.0)
	194
[135]	195	!zqev = MAX (0.0, (zqs(i,k)-q(i,k))*eeff1 )
	196	!zqevt = (zrfl(i)/l2c(i,k))*eeff2
	197	!zqev = MIN (zqev, zqevt)
	198	!zrfln(i) = zrfl(i) - zqev*l2c(i,k)
[253]	199	!zrfln(i) = zrfl(i) - 1.5e-5(1.0-q(i,k)/zqs(i,k))sqrt(zrfl(i))
	200	!zrfln(i) = min(zrfln(i),zrfl(i))
[135]	201	! this is what is actually written in the manual
	202
[253]	203	d_q(i,k) = - (zrfln(i)-zrfl(i))/l2c(i,k)*ptimestep
	204	!d_t(i,k) = d_q(i,k) * RLVTT/RCPD!/(1.0+RVTMP2*q(i,k)) ! double BC modif here
	205	d_t(i,k) = - d_q(i,k) * RLVTT/RCPD ! was bugged!
	206
[135]	207	zrfl(i) = zrfln(i)
	208	ENDIF
	209	ENDDO
	210	ENDIF
	211
	212	DO i = 1, ngridmx
	213	zoliq(i) = 0.0
	214	ENDDO
	215
	216
	217	if(simple)then
	218
	219	DO i = 1, ngridmx
[253]	220	ttemp = zt(i,k)
	221	IF (ttemp .ge. To) THEN
	222	lconvert=rainthreshold
	223	ELSEIF (ttemp .gt. t_crit) THEN
	224	lconvert=rainthreshold*(1.- t_crit/ttemp)
	225	lconvert=MAX(0.0,lconvert)
	226	ELSE
	227	lconvert=0.
	228	ENDIF
[135]	229
[253]	230
	231	IF (ql(i,k).gt.1.e-9) then
	232	zneb(i) = MAX(rneb(i,k), seuil_neb)
	233	IF ((ql(i,k)/zneb(i)).gt.lconvert)THEN ! precipitate!
[622]	234	d_ql(i,k) = -MAX((ql(i,k)-lconvert*zneb(i)),0.0)
[253]	235	zrfl(i) = zrfl(i) - d_ql(i,k)*l2c(i,k)/ptimestep
	236	ENDIF
	237	ENDIF
[135]	238	ENDDO
	239
	240	else
	241
	242	DO i = 1, ngridmx
	243	IF (rneb(i,k).GT.0.0) THEN
	244	zoliq(i) = ql(i,k)
[253]	245	zrho(i) = pplay(i,k) / ( zt(i,k) * R )
[135]	246	zdz(i) = (pplev(i,k)-pplev(i,k+1)) / (zrho(i)*g)
[253]	247	zfrac(i) = (zt(i,k)-ztglace) / (To-ztglace)
[135]	248	zfrac(i) = MAX(zfrac(i), 0.0)
	249	zfrac(i) = MIN(zfrac(i), 1.0)
	250	zneb(i) = MAX(rneb(i,k), seuil_neb)
	251	ENDIF
	252	ENDDO
	253
	254	DO n = 1, ninter
	255	DO i = 1, ngridmx
	256	IF (rneb(i,k).GT.0.0) THEN
[622]	257	zchau(i) = (ctptimestep/FLOAT(ninter)) zoliq(i) &
[135]	258	* (1.0-EXP(-(zoliq(i)/zneb(i)/cl)*2)) zfrac(i)
[253]	259
[135]	260	! this is the ONLY place where zneb, ct and cl are used
	261	zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
	262	zfroi(i) = ptimestep/FLOAT(ninter)/zdz(i)*zoliq(i) &
	263	fallv(zrhol(i)) (1.0-zfrac(i)) ! zfroi behaves oddly...
	264	! * 0.1 * (1.0-zfrac(i))
	265	ztot(i) = zchau(i) + zfroi(i)
	266
	267	IF (zneb(i).EQ.seuil_neb) ztot(i) = 0.0
	268	ztot(i) = MIN(MAX(ztot(i),0.0),zoliq(i))
	269	zoliq(i) = MAX(zoliq(i)-ztot(i), 0.0)
[253]	270
[135]	271	ENDIF
	272	ENDDO
	273	ENDDO
	274
	275	! Change in cloud density and surface H2O values
	276	DO i = 1, ngridmx
	277	IF (rneb(i,k).GT.0.0) THEN
[253]	278	d_ql(i,k) = (zoliq(i) - ql(i,k))!/ptimestep
	279	zrfl(i) = zrfl(i)+ MAX(ql(i,k)-zoliq(i),0.0)*l2c(i,k)/ptimestep
[135]	280	ENDIF
	281	ENDDO
	282
	283	endif ! if simple
	284
	285	9999 continue
	286
	287	! Rain or snow on the ground
	288	DO i = 1, ngridmx
[253]	289	if(zrfl(i).lt.0.0)then
	290	print*,'Droplets of negative rain are falling...'
	291	call abort
	292	endif
[135]	293	IF (t(i,1) .LT. To) THEN
	294	dqssnow(i) = zrfl(i)
[253]	295	dqsrain(i) = 0.0
[135]	296	ELSE
[253]	297	dqssnow(i) = 0.0
[135]	298	dqsrain(i) = zrfl(i) ! liquid water = ice for now
	299	ENDIF
	300	ENDDO
	301
	302	! now subroutine -----> GCM variables
	303	DO k = 1, nlayermx
	304	DO i = 1, ngridmx
	305
	306	if(evap_prec)then
	307	dqrain(i,k,i_vap) = d_q(i,k)/ptimestep
	308	d_t(i,k) = d_t(i,k)/ptimestep
	309	else
	310	dqrain(i,k,i_vap) = 0.0
	311	d_t(i,k) = 0.0
	312	endif
[253]	313	dqrain(i,k,i_ice) = d_ql(i,k)/ptimestep
[135]	314
	315	ENDDO
	316	ENDDO
	317
	318	RETURN
	319	end subroutine rain

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

Download in other formats:

Original Format