Context Navigation

rain.F90 @ 1243

Last change on this file since 1243 was 1016, checked in by jleconte, 11 years ago

07/08/2013 == JL

Water cycle in double precision (largescale+moistadj)
Improved wate rayleigh.
First step for rayleigh with variable species. Now, just need to change optcv.
changed some interpolation indices in callcorrk to limit dependency of OLR on the number of layers

File size: 14.0 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,save :: evap_prec ! 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	write(,) "re-evaporate precipitations?"
146	evap_prec=.true. ! default value
147	call getin("evap_prec",evap_prec)
148	write(,) " evap_prec = ",evap_prec
149
150	firstcall = .false.
151	ENDIF
152
153	! GCM -----> subroutine variables
154	DO k = 1, nlayermx
155	DO i = 1, ngrid
156
157	zt(i,k) = t(i,k)+pdt(i,k)*ptimestep ! a big fat bug was here
158	q(i,k) = pq(i,k,i_vap)+pdq(i,k,i_vap)*ptimestep
159	ql(i,k) = pq(i,k,i_ice)+pdq(i,k,i_ice)*ptimestep
160
161	!q(i,k) = pq(i,k,i_vap)!+pdq(i,k,i_vap)
162	!ql(i,k) = pq(i,k,i_ice)!+pdq(i,k,i_ice)
163
164	if(q(i,k).lt.0.)then ! if this is not done, we don't conserve water
165	q(i,k)=0.
166	endif
167	if(ql(i,k).lt.0.)then
168	ql(i,k)=0.
169	endif
170
171	ENDDO
172	ENDDO
173
174	! Initialise the outputs
175	DO k = 1, nlayermx
176	DO i = 1, ngrid
177	d_t(i,k) = 0.0
178	d_q(i,k) = 0.0
179	d_ql(i,k) = 0.0
180	ENDDO
181	ENDDO
182	DO i = 1, ngrid
183	zrfl(i) = 0.0
184	zrfln(i) = 0.0
185	ENDDO
186
187	! calculate saturation mixing ratio
188	DO k = 1, nlayermx
189	DO i = 1, ngrid
190	ttemp = zt(i,k)
191	ptemp = pplay(i,k)
192	! call watersat(ttemp,ptemp,zqs(i,k))
193	call Psat_water(ttemp,ptemp,psat_tmp,zqs(i,k))
194	call Tsat_water(ptemp,Tsat(i,k))
195	ENDDO
196	ENDDO
197
198	! get column / layer conversion factor
199	DO k = 1, nlayermx
200	DO i = 1, ngrid
201	l2c(i,k)=(pplev(i,k)-pplev(i,k+1))/g
202	ENDDO
203	ENDDO
204
205
206	! Vertical loop (from top to bottom)
207	! We carry the rain with us and calculate that added by warm/cold precipitation
208	! processes and that subtracted by evaporation at each level.
209	DO 9999 k = nlayermx, 1, -1
210
211	IF (evap_prec) THEN ! note no rneb dependence!
212	DO i = 1, ngrid
213	IF (zrfl(i) .GT.0.) THEN
214
215	if(zt(i,k).gt.Tsat(i,k))then
216	!! treat the case where all liquid water should boil
217	zqev=MIN((zt(i,k)-Tsat(i,k))RCPDl2c(i,k)/RLVTT/ptimestep,zrfl(i))
218	zrfl(i)=MAX(zrfl(i)-zqev,0.)
219	d_q(i,k)=zqev/l2c(i,k)*ptimestep
220	d_t(i,k) = - d_q(i,k) * RLVTT/RCPD
221	else
222	zqev = MAX (0.0, (zqs(i,k)-q(i,k)))*l2c(i,k)/ptimestep !there was a bug here
223	zqevt= 2.0e-5*(1.0-q(i,k)/zqs(i,k)) & !default was 2.e-5
224	sqrt(zrfl(i))l2c(i,k)/pplay(i,k)zt(i,k)R ! BC modif here
225	zqevt = MAX (zqevt, 0.0)
226	zqev = MIN (zqev, zqevt)
227	zqev = MAX (zqev, 0.0)
228	zrfln(i)= zrfl(i) - zqev
229	zrfln(i)= max(zrfln(i),0.0)
230
231	d_q(i,k) = - (zrfln(i)-zrfl(i))/l2c(i,k)*ptimestep
232	!d_t(i,k) = d_q(i,k) * RLVTT/RCPD!/(1.0+RVTMP2*q(i,k)) ! double BC modif here
233	d_t(i,k) = - d_q(i,k) * RLVTT/RCPD ! was bugged!
234	zrfl(i) = zrfln(i)
235	end if
236
237
238	ENDIF
239	ENDDO
240	ENDIF
241
242	DO i = 1, ngrid
243	zoliq(i) = 0.0
244	ENDDO
245
246
247	if(precip_scheme.eq.1)then
248
249	DO i = 1, ngrid
250	ttemp = zt(i,k)
251	IF (ttemp .ge. T_h2O_ice_liq) THEN
252	lconvert=rainthreshold
253	ELSEIF (ttemp .gt. t_crit) THEN
254	lconvert=rainthreshold*(1.- t_crit/ttemp)
255	lconvert=MAX(0.0,lconvert)
256	ELSE
257	lconvert=0.
258	ENDIF
259
260
261	IF (ql(i,k).gt.1.e-9) then
262	zneb(i) = MAX(rneb(i,k), seuil_neb)
263	IF ((ql(i,k)/zneb(i)).gt.lconvert)THEN ! precipitate!
264	d_ql(i,k) = -MAX((ql(i,k)-lconvert*zneb(i)),0.0)
265	zrfl(i) = zrfl(i) - d_ql(i,k)*l2c(i,k)/ptimestep
266	ENDIF
267	ENDIF
268	ENDDO
269
270	elseif (precip_scheme.ge.2) then
271
272	DO i = 1, ngrid
273	IF (rneb(i,k).GT.0.0) THEN
274	zoliq(i) = ql(i,k)
275	zrho(i) = pplay(i,k) / ( zt(i,k) * R )
276	zdz(i) = (pplev(i,k)-pplev(i,k+1)) / (zrho(i)*g)
277	zfrac(i) = (zt(i,k)-T_h2O_ice_clouds) / (T_h2O_ice_liq-T_h2O_ice_clouds)
278	zfrac(i) = MAX(zfrac(i), 0.0)
279	zfrac(i) = MIN(zfrac(i), 1.0)
280	zneb(i) = MAX(rneb(i,k), seuil_neb)
281	ENDIF
282	ENDDO
283
284	!recalculate liquid water particle radii
285	call h2o_cloudrad(ngrid,ql,reffh2oliq,reffh2oice)
286
287	SELECT CASE(precip_scheme)
288	!precip scheme from Sundquist 78
289	CASE(2)
290
291	DO n = 1, ninter
292	DO i = 1, ngrid
293	IF (rneb(i,k).GT.0.0) THEN
294	! this is the ONLY place where zneb, precip_timescale and cloud_sat are used
295
296	zchau(i) = (ptimestep/(FLOAT(ninter)precip_timescale)) zoliq(i) &
297	* (1.0-EXP(-(zoliq(i)/zneb(i)/cloud_sat)*2)) zfrac(i)
298	zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
299	zfroi(i) = ptimestep/FLOAT(ninter)/zdz(i)*zoliq(i) &
300	fall2v(reffh2oice(i,k)) (1.0-zfrac(i)) ! zfroi behaves oddly...
301	ztot(i) = zchau(i) + zfroi(i)
302
303	IF (zneb(i).EQ.seuil_neb) ztot(i) = 0.0
304	ztot(i) = MIN(MAX(ztot(i),0.0),zoliq(i))
305	zoliq(i) = MAX(zoliq(i)-ztot(i), 0.0)
306
307	ENDIF
308	ENDDO
309	ENDDO
310
311	!precip scheme modified from Sundquist 78 (in q**3)
312	CASE(3)
313
314	DO n = 1, ninter
315	DO i = 1, ngrid
316	IF (rneb(i,k).GT.0.0) THEN
317	! this is the ONLY place where zneb, precip_timescale and cloud_sat are used
318
319	zchau(i) = (ptimestep/(FLOAT(ninter)precip_timescalecloud_sat*2)) (zoliq(i)/zneb(i))**3
320	zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
321	zfroi(i) = ptimestep/FLOAT(ninter)/zdz(i)*zoliq(i) &
322	fall2v(reffh2oice(i,k)) (1.0-zfrac(i)) ! zfroi behaves oddly...
323	ztot(i) = zchau(i) + zfroi(i)
324
325	IF (zneb(i).EQ.seuil_neb) ztot(i) = 0.0
326	ztot(i) = MIN(MAX(ztot(i),0.0),zoliq(i))
327	zoliq(i) = MAX(zoliq(i)-ztot(i), 0.0)
328
329	ENDIF
330	ENDDO
331	ENDDO
332
333	!precip scheme modified from Boucher 95
334	CASE(4)
335
336	DO n = 1, ninter
337	DO i = 1, ngrid
338	IF (rneb(i,k).GT.0.0) THEN
339	! this is the ONLY place where zneb and c1 are used
340
341	zchau(i) = ptimestep/FLOAT(ninter) c1 zrho(i) &
342	(zoliq(i)/zneb(i))2reffh2oliq(i,k)zneb(i) zfrac(i)
343	zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
344	zfroi(i) = ptimestep/FLOAT(ninter)/zdz(i)*zoliq(i) &
345	fall2v(reffh2oice(i,k)) (1.0-zfrac(i)) ! zfroi behaves oddly...
346	ztot(i) = zchau(i) + zfroi(i)
347
348	IF (zneb(i).EQ.seuil_neb) ztot(i) = 0.0
349	ztot(i) = MIN(MAX(ztot(i),0.0),zoliq(i))
350	zoliq(i) = MAX(zoliq(i)-ztot(i), 0.0)
351
352	ENDIF
353	ENDDO
354	ENDDO
355
356	END SELECT ! precip_scheme
357
358	! Change in cloud density and surface H2O values
359	DO i = 1, ngrid
360	IF (rneb(i,k).GT.0.0) THEN
361	d_ql(i,k) = (zoliq(i) - ql(i,k))!/ptimestep
362	zrfl(i) = zrfl(i)+ MAX(ql(i,k)-zoliq(i),0.0)*l2c(i,k)/ptimestep
363	ENDIF
364	ENDDO
365
366
367	endif ! if precip_scheme=1
368
369	9999 continue
370
371	! Rain or snow on the ground
372	DO i = 1, ngrid
373	if(zrfl(i).lt.0.0)then
374	print*,'Droplets of negative rain are falling...'
375	call abort
376	endif
377	IF (t(i,1) .LT. T_h2O_ice_liq) THEN
378	dqssnow(i) = zrfl(i)
379	dqsrain(i) = 0.0
380	ELSE
381	dqssnow(i) = 0.0
382	dqsrain(i) = zrfl(i) ! liquid water = ice for now
383	ENDIF
384	ENDDO
385
386	! now subroutine -----> GCM variables
387	DO k = 1, nlayermx
388	DO i = 1, ngrid
389
390	if(evap_prec)then
391	dqrain(i,k,i_vap) = d_q(i,k)/ptimestep
392	d_t(i,k) = d_t(i,k)/ptimestep
393	else
394	dqrain(i,k,i_vap) = 0.0
395	d_t(i,k) = 0.0
396	endif
397	dqrain(i,k,i_ice) = d_ql(i,k)/ptimestep
398
399	ENDDO
400	ENDDO
401
402	RETURN
403	end subroutine rain

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

Download in other formats:

Original Format