Context Navigation

rain.F90 @ 937

Last change on this file since 937 was 863, checked in by jleconte, 13 years ago

23/01/2013 == JL

Correction in largescale. a rneb factor was forgotten
Added some spectra in ave_stelpec
Corrected reevaporation in rain. Now conserve water better

File size: 13.8 KB

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

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

Download in other formats:

Original Format