Context Navigation

rain.F90 @ 775

Last change on this file since 775 was 731, checked in by jleconte, 13 years ago
Corrected precipitation evaporation scheme + snow fall (JL+BC) Commented some unecessary writediagfi calls that where used for test.
File size: 13.8 KB

Line
1	subroutine rain(ptimestep,pplev,pplay,t,pdt,pq,pdq,d_t,dqrain,dqsrain,dqssnow,rneb)
2
3
4	! to use 'getin'
5	use ioipsl_getincom
6	use watercommon_h, only: T_h2O_ice_liq,T_h2O_ice_clouds, RLVTT, RCPD, RCPV, RV, RVTMP2,Psat_water,Tsat_water,rhowater
7	use radii_mod, only: h2o_cloudrad
8	implicit none
9
10	!==================================================================
11	!
12	! Purpose
13	! -------
14	! Calculates H2O precipitation using simplified microphysics.
15	!
16	! Authors
17	! -------
18	! Adapted from the LMDTERRE code by R. Wordsworth (2009)
19	! Added rain vaporization in case of T>Tsat
20	! Original author Z. X. Li (1993)
21	!
22	!==================================================================
23
24	#include "dimensions.h"
25	#include "dimphys.h"
26	#include "tracer.h"
27	#include "comcstfi.h"
28	#include "callkeys.h"
29
30	! Pre-arguments (for universal model)
31	real pq(ngridmx,nlayermx,nqmx) ! tracer (kg/kg)
32	real qsurf(ngridmx,nqmx) ! tracer at the surface (kg.m-2)
33	REAL pdt(ngridmx,nlayermx),pdq(ngridmx,nlayermx,nqmx)
34
35	real dqrain(ngridmx,nlayermx,nqmx) ! tendency of H2O precipitation (kg/kg.s-1)
36	real dqsrain(ngridmx) ! rain flux at the surface (kg.m-2.s-1)
37	real dqssnow(ngridmx) ! snow flux at the surface (kg.m-2.s-1)
38	REAL d_t(ngridmx,nlayermx) ! temperature increment
39
40	! Arguments
41	REAL ptimestep ! time interval
42	REAL pplev(ngridmx,nlayermx+1) ! inter-layer pressure
43	REAL pplay(ngridmx,nlayermx) ! mid-layer pressure
44	REAL t(ngridmx,nlayermx) ! input temperature (K)
45	REAL zt(ngridmx,nlayermx) ! working temperature (K)
46	REAL ql(ngridmx,nlayermx) ! liquid water (Kg/Kg)
47	REAL q(ngridmx,nlayermx) ! specific humidity (Kg/Kg)
48	REAL rneb(ngridmx,nlayermx) ! cloud fraction
49	REAL d_q(ngridmx,nlayermx) ! water vapor increment
50	REAL d_ql(ngridmx,nlayermx) ! liquid water / ice increment
51
52	! Subroutine options
53	REAL seuil_neb ! Nebulosity threshold
54	PARAMETER (seuil_neb=0.001)
55
56	INTEGER,save :: precip_scheme ! id number for precipitaion scheme
57	! for simple scheme (precip_scheme=1)
58	REAL,SAVE :: rainthreshold ! Precipitation threshold in simple scheme
59	! for sundquist scheme (precip_scheme=2-3)
60	REAL,SAVE :: cloud_sat ! Precipitation threshold in non simple scheme
61	REAL,SAVE :: precip_timescale ! Precipitation timescale
62	! for Boucher scheme (precip_scheme=4)
63	REAL,SAVE :: Cboucher ! Precipitation constant in Boucher 95 scheme
64	REAL,PARAMETER :: Kboucher=1.19E8
65	REAL,SAVE :: c1
66
67	INTEGER ninter
68	PARAMETER (ninter=5)
69
70	logical evap_prec ! Does the rain evaporate?
71	parameter(evap_prec=.true.)
72
73	! for simple scheme
74	real t_crit
75	PARAMETER (t_crit=218.0)
76	real lconvert
77
78	! Local variables
79	INTEGER i, k, n
80	REAL zqs(ngridmx,nlayermx),Tsat(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 reffh2oliq(ngridmx,nlayermx),reffh2oice(ngridmx,nlayermx)
88
89	real ttemp, ptemp, psat_tmp
90	real tnext(ngridmx,nlayermx)
91
92	real l2c(ngridmx,nlayermx)
93	real dWtot
94
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, fall2v, zzz ! falling speed of ice crystals
105	fallv (zzz) = 3.29 * ((zzz)**0.16)
106	fall2v (zzz) =10.6 * ((zzz)**0.31) !for use with radii
107
108	DATA firstcall /.true./
109
110	IF (firstcall) THEN
111
112	i_vap=igcm_h2o_vap
113	i_ice=igcm_h2o_ice
114
115	write(,) "rain: i_ice=",i_ice
116	write(,) " i_vap=",i_vap
117
118	PRINT*, 'in rain.F, ninter=', ninter
119	PRINT*, 'in rain.F, evap_prec=', evap_prec
120
121	write(,) "Precipitation scheme to use?"
122	precip_scheme=1 ! default value
123	call getin("precip_scheme",precip_scheme)
124	write(,) " precip_scheme = ",precip_scheme
125
126	if (precip_scheme.eq.1) then
127	write(,) "rainthreshold in simple scheme?"
128	rainthreshold=0. ! default value
129	call getin("rainthreshold",rainthreshold)
130	write(,) " rainthreshold = ",rainthreshold
131
132	else if (precip_scheme.eq.2.or.precip_scheme.eq.3) then
133	write(,) "cloud water saturation level in non simple scheme?"
134	cloud_sat=2.6e-4 ! default value
135	call getin("cloud_sat",cloud_sat)
136	write(,) " cloud_sat = ",cloud_sat
137	write(,) "precipitation timescale in non simple scheme?"
138	precip_timescale=3600. ! default value
139	call getin("precip_timescale",precip_timescale)
140	write(,) " precip_timescale = ",precip_timescale
141
142	else if (precip_scheme.eq.4) then
143	write(,) "multiplicative constant in Boucher 95 precip scheme"
144	Cboucher=1. ! default value
145	call getin("Cboucher",Cboucher)
146	write(,) " Cboucher = ",Cboucher
147	c1=1.001.097/rhowaterCboucher*Kboucher
148
149	endif
150
151	firstcall = .false.
152	ENDIF
153
154	! GCM -----> subroutine variables
155	DO k = 1, nlayermx
156	DO i = 1, ngridmx
157
158	zt(i,k) = t(i,k)+pdt(i,k)*ptimestep ! a big fat bug was here
159	q(i,k) = pq(i,k,i_vap)+pdq(i,k,i_vap)*ptimestep
160	ql(i,k) = pq(i,k,i_ice)+pdq(i,k,i_ice)*ptimestep
161
162	!q(i,k) = pq(i,k,i_vap)!+pdq(i,k,i_vap)
163	!ql(i,k) = pq(i,k,i_ice)!+pdq(i,k,i_ice)
164
165	if(q(i,k).lt.0.)then ! if this is not done, we don't conserve water
166	q(i,k)=0.
167	endif
168	if(ql(i,k).lt.0.)then
169	ql(i,k)=0.
170	endif
171
172	ENDDO
173	ENDDO
174
175	! Initialise the outputs
176	DO k = 1, nlayermx
177	DO i = 1, ngridmx
178	d_t(i,k) = 0.0
179	d_q(i,k) = 0.0
180	d_ql(i,k) = 0.0
181	ENDDO
182	ENDDO
183	DO i = 1, ngridmx
184	zrfl(i) = 0.0
185	zrfln(i) = 0.0
186	ENDDO
187
188	! calculate saturation mixing ratio
189	DO k = 1, nlayermx
190	DO i = 1, ngridmx
191	ttemp = zt(i,k)
192	ptemp = pplay(i,k)
193	! call watersat(ttemp,ptemp,zqs(i,k))
194	call Psat_water(ttemp,ptemp,psat_tmp,zqs(i,k))
195	call Tsat_water(ptemp,Tsat(i,k))
196	ENDDO
197	ENDDO
198
199	! get column / layer conversion factor
200	DO k = 1, nlayermx
201	DO i = 1, ngridmx
202	l2c(i,k)=(pplev(i,k)-pplev(i,k+1))/g
203	ENDDO
204	ENDDO
205
206
207	! Vertical loop (from top to bottom)
208	! We carry the rain with us and calculate that added by warm/cold precipitation
209	! processes and that subtracted by evaporation at each level.
210	DO 9999 k = nlayermx, 1, -1
211
212	IF (evap_prec) THEN ! note no rneb dependence!
213	DO i = 1, ngridmx
214	IF (zrfl(i) .GT.0.) THEN
215
216	if(zt(i,k).gt.Tsat(i,k))then
217	! treat the case where all liquid water should boil
218	zqev=MIN((zt(i,k)-Tsat(i,k))RCPDl2c(i,k)/RLVTT,zrfl(i))
219	zrfl(i)=MAX(zrfl(i)-zqev,0.)
220	d_q(i,k)=zqev/l2c(i,k)
221	d_t(i,k) = - d_q(i,k) * RLVTT/RCPD
222	else
223	zqev = MAX (0.0, (zqs(i,k)-q(i,k)))*l2c(i,k)/ptimestep !there was a bug here
224	zqevt= 2.0e-5*(1.0-q(i,k)/zqs(i,k)) & !default was 2.e-5
225	sqrt(zrfl(i))l2c(i,k)/pplay(i,k)zt(i,k)R ! BC modif here
226	zqevt = MAX (zqevt, 0.0)
227	zqev = MIN (zqev, zqevt)
228	zqev = MAX (zqev, 0.0)
229	zrfln(i)= zrfl(i) - zqev
230	zrfln(i)= max(zrfln(i),0.0)
231
232	d_q(i,k) = - (zrfln(i)-zrfl(i))/l2c(i,k)*ptimestep
233	!d_t(i,k) = d_q(i,k) * RLVTT/RCPD!/(1.0+RVTMP2*q(i,k)) ! double BC modif here
234	d_t(i,k) = - d_q(i,k) * RLVTT/RCPD ! was bugged!
235	zrfl(i) = zrfln(i)
236	end if
237
238
239	ENDIF
240	ENDDO
241	ENDIF
242
243	DO i = 1, ngridmx
244	zoliq(i) = 0.0
245	ENDDO
246
247
248	if(precip_scheme.eq.1)then
249
250	DO i = 1, ngridmx
251	ttemp = zt(i,k)
252	IF (ttemp .ge. T_h2O_ice_liq) THEN
253	lconvert=rainthreshold
254	ELSEIF (ttemp .gt. t_crit) THEN
255	lconvert=rainthreshold*(1.- t_crit/ttemp)
256	lconvert=MAX(0.0,lconvert)
257	ELSE
258	lconvert=0.
259	ENDIF
260
261
262	IF (ql(i,k).gt.1.e-9) then
263	zneb(i) = MAX(rneb(i,k), seuil_neb)
264	IF ((ql(i,k)/zneb(i)).gt.lconvert)THEN ! precipitate!
265	d_ql(i,k) = -MAX((ql(i,k)-lconvert*zneb(i)),0.0)
266	zrfl(i) = zrfl(i) - d_ql(i,k)*l2c(i,k)/ptimestep
267	ENDIF
268	ENDIF
269	ENDDO
270
271	elseif (precip_scheme.ge.2) then
272
273	DO i = 1, ngridmx
274	IF (rneb(i,k).GT.0.0) THEN
275	zoliq(i) = ql(i,k)
276	zrho(i) = pplay(i,k) / ( zt(i,k) * R )
277	zdz(i) = (pplev(i,k)-pplev(i,k+1)) / (zrho(i)*g)
278	zfrac(i) = (zt(i,k)-T_h2O_ice_clouds) / (T_h2O_ice_liq-T_h2O_ice_clouds)
279	zfrac(i) = MAX(zfrac(i), 0.0)
280	zfrac(i) = MIN(zfrac(i), 1.0)
281	zneb(i) = MAX(rneb(i,k), seuil_neb)
282	ENDIF
283	ENDDO
284
285	!recalculate liquid water particle radii
286	call h2o_cloudrad(ql,reffh2oliq,reffh2oice)
287
288	SELECT CASE(precip_scheme)
289	!precip scheme from Sundquist 78
290	CASE(2)
291
292	DO n = 1, ninter
293	DO i = 1, ngridmx
294	IF (rneb(i,k).GT.0.0) THEN
295	! this is the ONLY place where zneb, precip_timescale and cloud_sat are used
296
297	zchau(i) = (ptimestep/(FLOAT(ninter)precip_timescale)) zoliq(i) &
298	* (1.0-EXP(-(zoliq(i)/zneb(i)/cloud_sat)*2)) zfrac(i)
299	zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
300	zfroi(i) = ptimestep/FLOAT(ninter)/zdz(i)*zoliq(i) &
301	fall2v(reffh2oice(i,k)) (1.0-zfrac(i)) ! zfroi behaves oddly...
302	ztot(i) = zchau(i) + zfroi(i)
303
304	IF (zneb(i).EQ.seuil_neb) ztot(i) = 0.0
305	ztot(i) = MIN(MAX(ztot(i),0.0),zoliq(i))
306	zoliq(i) = MAX(zoliq(i)-ztot(i), 0.0)
307
308	ENDIF
309	ENDDO
310	ENDDO
311
312	!precip scheme modified from Sundquist 78 (in q**3)
313	CASE(3)
314
315	DO n = 1, ninter
316	DO i = 1, ngridmx
317	IF (rneb(i,k).GT.0.0) THEN
318	! this is the ONLY place where zneb, precip_timescale and cloud_sat are used
319
320	zchau(i) = (ptimestep/(FLOAT(ninter)precip_timescalecloud_sat*2)) (zoliq(i)/zneb(i))**3
321	zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
322	zfroi(i) = ptimestep/FLOAT(ninter)/zdz(i)*zoliq(i) &
323	fall2v(reffh2oice(i,k)) (1.0-zfrac(i)) ! zfroi behaves oddly...
324	ztot(i) = zchau(i) + zfroi(i)
325
326	IF (zneb(i).EQ.seuil_neb) ztot(i) = 0.0
327	ztot(i) = MIN(MAX(ztot(i),0.0),zoliq(i))
328	zoliq(i) = MAX(zoliq(i)-ztot(i), 0.0)
329
330	ENDIF
331	ENDDO
332	ENDDO
333
334	!precip scheme modified from Boucher 95
335	CASE(4)
336
337	DO n = 1, ninter
338	DO i = 1, ngridmx
339	IF (rneb(i,k).GT.0.0) THEN
340	! this is the ONLY place where zneb and c1 are used
341
342	zchau(i) = ptimestep/FLOAT(ninter) c1 zrho(i) &
343	(zoliq(i)/zneb(i))2reffh2oliq(i,k)zneb(i) zfrac(i)
344	zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
345	zfroi(i) = ptimestep/FLOAT(ninter)/zdz(i)*zoliq(i) &
346	fall2v(reffh2oice(i,k)) (1.0-zfrac(i)) ! zfroi behaves oddly...
347	ztot(i) = zchau(i) + zfroi(i)
348
349	IF (zneb(i).EQ.seuil_neb) ztot(i) = 0.0
350	ztot(i) = MIN(MAX(ztot(i),0.0),zoliq(i))
351	zoliq(i) = MAX(zoliq(i)-ztot(i), 0.0)
352
353	ENDIF
354	ENDDO
355	ENDDO
356
357	END SELECT ! precip_scheme
358
359	! Change in cloud density and surface H2O values
360	DO i = 1, ngridmx
361	IF (rneb(i,k).GT.0.0) THEN
362	d_ql(i,k) = (zoliq(i) - ql(i,k))!/ptimestep
363	zrfl(i) = zrfl(i)+ MAX(ql(i,k)-zoliq(i),0.0)*l2c(i,k)/ptimestep
364	ENDIF
365	ENDDO
366
367
368	endif ! if precip_scheme=1
369
370	9999 continue
371
372	! Rain or snow on the ground
373	DO i = 1, ngridmx
374	if(zrfl(i).lt.0.0)then
375	print*,'Droplets of negative rain are falling...'
376	call abort
377	endif
378	IF (t(i,1) .LT. T_h2O_ice_liq) THEN
379	dqssnow(i) = zrfl(i)
380	dqsrain(i) = 0.0
381	ELSE
382	dqssnow(i) = 0.0
383	dqsrain(i) = zrfl(i) ! liquid water = ice for now
384	ENDIF
385	ENDDO
386
387	! now subroutine -----> GCM variables
388	DO k = 1, nlayermx
389	DO i = 1, ngridmx
390
391	if(evap_prec)then
392	dqrain(i,k,i_vap) = d_q(i,k)/ptimestep
393	d_t(i,k) = d_t(i,k)/ptimestep
394	else
395	dqrain(i,k,i_vap) = 0.0
396	d_t(i,k) = 0.0
397	endif
398	dqrain(i,k,i_ice) = d_ql(i,k)/ptimestep
399
400	ENDDO
401	ENDDO
402
403	RETURN
404	end subroutine rain

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

Download in other formats:

Original Format