Context Navigation

rain.F90 @ 524

Last change on this file since 524 was 253, checked in by emillour, 14 years ago

Generic GCM

Massive update to version 0.7

EM+RW

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

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

Download in other formats:

Original Format