Context Navigation

rain.F90 @ 1477

Last change on this file since 1477 was 1397, checked in by milmd, 10 years ago
In LMDZ.GENERIC replacement of all phystd .h files by module files.
File size: 14.1 KB

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

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

Download in other formats:

Original Format