Context Navigation

rain.F90 @ 135

Last change on this file since 135 was 135, checked in by aslmd, 14 years ago
CHANGEMENT ARBORESCENCE ETAPE 2 -- NON COMPLET
File size: 9.5 KB

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

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

Download in other formats:

Original Format