Context Navigation

rain.F90 @ 650

Last change on this file since 650 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

Line
1	subroutine rain(ptimestep,pplev,pplay,t,pdt,pq,pdq,d_t,dqrain,dqsrain,dqssnow,rneb)
2
3
4	use watercommon_h, only: T_h2O_ice_liq,T_h2O_ice_clouds, 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) ! input temperature (K)
42	REAL zt(ngridmx,nlayermx) ! working temperature (K)
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
53	REAL,PARAMETER :: cl=2.0e-4 ! Precipitation threshold
54	REAL,PARAMETER :: ct=1./1800. ! Precipitation rate
55
56	INTEGER ninter
57	PARAMETER (ninter=5)
58
59	logical simple ! Use very simple Emanuel scheme?
60	parameter(simple=.true.)
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
90	real l2c(ngridmx,nlayermx)
91	real dWtot
92
93
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
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
128
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
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
158	ttemp = zt(i,k)
159	ptemp = pplay(i,k)
160	call watersat(ttemp,ptemp,zqs(i,k))
161	ENDDO
162	ENDDO
163
164	! get column / layer conversion factor
165	DO k = 1, nlayermx
166	DO i = 1, ngridmx
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
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
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
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	!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))
197	! this is what is actually written in the manual
198
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
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 = zt(i,k)
217	IF (ttemp .ge. T_h2O_ice_liq) THEN
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
225
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!
230	d_ql(i,k) = -MAX((ql(i,k)-lconvert*zneb(i)),0.0)
231	zrfl(i) = zrfl(i) - d_ql(i,k)*l2c(i,k)/ptimestep
232	ENDIF
233	ENDIF
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)
241	zrho(i) = pplay(i,k) / ( zt(i,k) * R )
242	zdz(i) = (pplev(i,k)-pplev(i,k+1)) / (zrho(i)*g)
243	zfrac(i) = (zt(i,k)-T_h2O_ice_clouds) / (T_h2O_ice_liq-T_h2O_ice_clouds)
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
253	zchau(i) = (ctptimestep/FLOAT(ninter)) zoliq(i) &
254	* (1.0-EXP(-(zoliq(i)/zneb(i)/cl)*2)) zfrac(i)
255
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)
266
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
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
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
285	if(zrfl(i).lt.0.0)then
286	print*,'Droplets of negative rain are falling...'
287	call abort
288	endif
289	IF (t(i,1) .LT. T_h2O_ice_liq) THEN
290	dqssnow(i) = zrfl(i)
291	dqsrain(i) = 0.0
292	ELSE
293	dqssnow(i) = 0.0
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
309	dqrain(i,k,i_ice) = d_ql(i,k)/ptimestep
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