Context Navigation

rain.F90 @ 1834

Last change on this file since 1834 was 1521, checked in by emillour, 9 years ago

All GCMs: Updates to make planetary codes (+Earth) setups converge.

Made a "phy_common" directory in libf, to contain routines common (wrt structural nature of underlying code/grid) to all LMDZ-related physics packages.
moved all "mod_phys_*" and "mod_grid_phy_lmdz" files from dynlonlat_phylonlat to "phy_common"
moved "ioipsl_getincom_p.F90 from "misc" to "phy_common" and modified it to match Earth GCM version and renamed it ioipsl_getin_p_mod.F90
added an "abort_physics" (as in Earth GCM) in "phy_common"
added a "print_control_mod.F90 (as in Earth GCM) in phy_common
made similar changes in LMDZ.GENERIC and LMDZ.MARS

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	use ioipsl_getin_p_mod, only: getin_p
5	use watercommon_h, only: T_h2O_ice_liq,T_h2O_ice_clouds, RLVTT, RCPD, RCPV, RV, RVTMP2,Psat_water,Tsat_water,rhowater
6	use radii_mod, only: h2o_cloudrad
7	USE tracer_h, only: igcm_h2o_vap, igcm_h2o_ice
8	use comcstfi_mod, only: g, r
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	! Arguments
26	integer,intent(in) :: ngrid ! number of atmospheric columns
27	integer,intent(in) :: nlayer ! number of atmospheric layers
28	integer,intent(in) :: nq ! number of tracers
29	real,intent(in) :: ptimestep ! time interval
30	real,intent(in) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa)
31	real,intent(in) :: pplay(ngrid,nlayer) ! mid-layer pressure (Pa)
32	real,intent(in) :: t(ngrid,nlayer) ! input temperature (K)
33	real,intent(in) :: pdt(ngrid,nlayer) ! input tendency on temperature (K/s)
34	real,intent(in) :: pq(ngrid,nlayer,nq) ! tracers (kg/kg)
35	real,intent(in) :: pdq(ngrid,nlayer,nq) ! input tendency on tracers
36	real,intent(out) :: d_t(ngrid,nlayer) ! temperature tendency (K/s)
37	real,intent(out) :: dqrain(ngrid,nlayer,nq) ! tendency of H2O precipitation (kg/kg.s-1)
38	real,intent(out) :: dqsrain(ngrid) ! rain flux at the surface (kg.m-2.s-1)
39	real,intent(out) :: dqssnow(ngrid) ! snow flux at the surface (kg.m-2.s-1)
40	real,intent(in) :: rneb(ngrid,nlayer) ! cloud fraction
41
42	REAL zt(ngrid,nlayer) ! working temperature (K)
43	REAL ql(ngrid,nlayer) ! liquid water (Kg/Kg)
44	REAL q(ngrid,nlayer) ! specific humidity (Kg/Kg)
45	REAL d_q(ngrid,nlayer) ! water vapor increment
46	REAL d_ql(ngrid,nlayer) ! liquid water / ice increment
47
48	! Subroutine options
49	REAL,PARAMETER :: seuil_neb=0.001 ! Nebulosity threshold
50
51	INTEGER,save :: precip_scheme ! id number for precipitaion scheme
52	! for simple scheme (precip_scheme=1)
53	REAL,SAVE :: rainthreshold ! Precipitation threshold in simple scheme
54	! for sundquist scheme (precip_scheme=2-3)
55	REAL,SAVE :: cloud_sat ! Precipitation threshold in non simple scheme
56	REAL,SAVE :: precip_timescale ! Precipitation timescale
57	! for Boucher scheme (precip_scheme=4)
58	REAL,SAVE :: Cboucher ! Precipitation constant in Boucher 95 scheme
59	REAL,PARAMETER :: Kboucher=1.19E8
60	REAL,SAVE :: c1
61	!$OMP THREADPRIVATE(precip_scheme,rainthreshold,cloud_sat,precip_timescale,Cboucher,c1)
62
63	INTEGER,PARAMETER :: ninter=5
64
65	logical,save :: evap_prec ! Does the rain evaporate?
66	!$OMP THREADPRIVATE(evap_prec)
67
68	! for simple scheme
69	real,parameter :: t_crit=218.0
70	real lconvert
71
72	! Local variables
73	INTEGER i, k, n
74	REAL zqs(ngrid,nlayer),Tsat(ngrid,nlayer), zdelta, zcor
75	REAL zrfl(ngrid), zrfln(ngrid), zqev, zqevt
76
77	REAL zoliq(ngrid)
78	REAL zdz(ngrid),zrho(ngrid),ztot(ngrid), zrhol(ngrid)
79	REAL zchau(ngrid),zfroi(ngrid),zfrac(ngrid),zneb(ngrid)
80
81	real reffh2oliq(ngrid,nlayer),reffh2oice(ngrid,nlayer)
82
83	real ttemp, ptemp, psat_tmp
84	real tnext(ngrid,nlayer)
85
86	real l2c(ngrid,nlayer)
87	real dWtot
88
89
90	! Indices of water vapour and water ice tracers
91	INTEGER, SAVE :: i_vap=0 ! water vapour
92	INTEGER, SAVE :: i_ice=0 ! water ice
93	!$OMP THREADPRIVATE(i_vap,i_ice)
94
95	LOGICAL,SAVE :: firstcall=.true.
96	!$OMP THREADPRIVATE(firstcall)
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_p("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_p("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_p("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_p("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_p("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_p("evap_prec",evap_prec)
148	write(,) " evap_prec = ",evap_prec
149
150	firstcall = .false.
151	ENDIF ! of IF (firstcall)
152
153	! GCM -----> subroutine variables
154	DO k = 1, nlayer
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	d_t(1:ngrid,1:nlayer) = 0.0
176	d_q(1:ngrid,1:nlayer) = 0.0
177	d_ql(1:ngrid,1:nlayer) = 0.0
178	zrfl(1:ngrid) = 0.0
179	zrfln(1:ngrid) = 0.0
180
181	! calculate saturation mixing ratio
182	DO k = 1, nlayer
183	DO i = 1, ngrid
184	ttemp = zt(i,k)
185	ptemp = pplay(i,k)
186	! call watersat(ttemp,ptemp,zqs(i,k))
187	call Psat_water(ttemp,ptemp,psat_tmp,zqs(i,k))
188	call Tsat_water(ptemp,Tsat(i,k))
189	ENDDO
190	ENDDO
191
192	! get column / layer conversion factor
193	DO k = 1, nlayer
194	DO i = 1, ngrid
195	l2c(i,k)=(pplev(i,k)-pplev(i,k+1))/g
196	ENDDO
197	ENDDO
198
199	! Vertical loop (from top to bottom)
200	! We carry the rain with us and calculate that added by warm/cold precipitation
201	! processes and that subtracted by evaporation at each level.
202	DO k = nlayer, 1, -1
203
204	IF (evap_prec) THEN ! note no rneb dependence!
205	DO i = 1, ngrid
206	IF (zrfl(i) .GT.0.) THEN
207
208	if(zt(i,k).gt.Tsat(i,k))then
209	!! treat the case where all liquid water should boil
210	zqev=MIN((zt(i,k)-Tsat(i,k))RCPDl2c(i,k)/RLVTT/ptimestep,zrfl(i))
211	zrfl(i)=MAX(zrfl(i)-zqev,0.)
212	d_q(i,k)=zqev/l2c(i,k)*ptimestep
213	d_t(i,k) = - d_q(i,k) * RLVTT/RCPD
214	else
215	zqev = MAX (0.0, (zqs(i,k)-q(i,k)))*l2c(i,k)/ptimestep !there was a bug here
216	zqevt= 2.0e-5*(1.0-q(i,k)/zqs(i,k)) & !default was 2.e-5
217	sqrt(zrfl(i))l2c(i,k)/pplay(i,k)zt(i,k)R ! BC modif here
218	zqevt = MAX (zqevt, 0.0)
219	zqev = MIN (zqev, zqevt)
220	zqev = MAX (zqev, 0.0)
221	zrfln(i)= zrfl(i) - zqev
222	zrfln(i)= max(zrfln(i),0.0)
223
224	d_q(i,k) = - (zrfln(i)-zrfl(i))/l2c(i,k)*ptimestep
225	!d_t(i,k) = d_q(i,k) * RLVTT/RCPD!/(1.0+RVTMP2*q(i,k)) ! double BC modif here
226	d_t(i,k) = - d_q(i,k) * RLVTT/RCPD ! was bugged!
227	zrfl(i) = zrfln(i)
228	end if
229
230
231	ENDIF ! of IF (zrfl(i) .GT.0.)
232	ENDDO
233	ENDIF ! of IF (evap_prec)
234
235	zoliq(1:ngrid) = 0.0
236
237
238	if(precip_scheme.eq.1)then
239
240	DO i = 1, ngrid
241	ttemp = zt(i,k)
242	IF (ttemp .ge. T_h2O_ice_liq) THEN
243	lconvert=rainthreshold
244	ELSEIF (ttemp .gt. t_crit) THEN
245	lconvert=rainthreshold*(1.- t_crit/ttemp)
246	lconvert=MAX(0.0,lconvert)
247	ELSE
248	lconvert=0.
249	ENDIF
250
251
252	IF (ql(i,k).gt.1.e-9) then
253	zneb(i) = MAX(rneb(i,k), seuil_neb)
254	IF ((ql(i,k)/zneb(i)).gt.lconvert)THEN ! precipitate!
255	d_ql(i,k) = -MAX((ql(i,k)-lconvert*zneb(i)),0.0)
256	zrfl(i) = zrfl(i) - d_ql(i,k)*l2c(i,k)/ptimestep
257	ENDIF
258	ENDIF
259	ENDDO
260
261	elseif (precip_scheme.ge.2) then
262
263	DO i = 1, ngrid
264	IF (rneb(i,k).GT.0.0) THEN
265	zoliq(i) = ql(i,k)
266	zrho(i) = pplay(i,k) / ( zt(i,k) * R )
267	zdz(i) = (pplev(i,k)-pplev(i,k+1)) / (zrho(i)*g)
268	zfrac(i) = (zt(i,k)-T_h2O_ice_clouds) / (T_h2O_ice_liq-T_h2O_ice_clouds)
269	zfrac(i) = MAX(zfrac(i), 0.0)
270	zfrac(i) = MIN(zfrac(i), 1.0)
271	zneb(i) = MAX(rneb(i,k), seuil_neb)
272	ENDIF
273	ENDDO
274
275	!recalculate liquid water particle radii
276	call h2o_cloudrad(ngrid,nlayer,ql,reffh2oliq,reffh2oice)
277
278	SELECT CASE(precip_scheme)
279	!precip scheme from Sundquist 78
280	CASE(2)
281
282	DO n = 1, ninter
283	DO i = 1, ngrid
284	IF (rneb(i,k).GT.0.0) THEN
285	! this is the ONLY place where zneb, precip_timescale and cloud_sat are used
286
287	zchau(i) = (ptimestep/(FLOAT(ninter)precip_timescale)) zoliq(i) &
288	* (1.0-EXP(-(zoliq(i)/zneb(i)/cloud_sat)*2)) zfrac(i)
289	zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
290	zfroi(i) = ptimestep/FLOAT(ninter)/zdz(i)*zoliq(i) &
291	fall2v(reffh2oice(i,k)) (1.0-zfrac(i)) ! zfroi behaves oddly...
292	ztot(i) = zchau(i) + zfroi(i)
293
294	IF (zneb(i).EQ.seuil_neb) ztot(i) = 0.0
295	ztot(i) = MIN(MAX(ztot(i),0.0),zoliq(i))
296	zoliq(i) = MAX(zoliq(i)-ztot(i), 0.0)
297
298	ENDIF
299	ENDDO
300	ENDDO
301
302	!precip scheme modified from Sundquist 78 (in q**3)
303	CASE(3)
304
305	DO n = 1, ninter
306	DO i = 1, ngrid
307	IF (rneb(i,k).GT.0.0) THEN
308	! this is the ONLY place where zneb, precip_timescale and cloud_sat are used
309
310	zchau(i) = (ptimestep/(FLOAT(ninter)precip_timescalecloud_sat*2)) (zoliq(i)/zneb(i))**3
311	zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
312	zfroi(i) = ptimestep/FLOAT(ninter)/zdz(i)*zoliq(i) &
313	fall2v(reffh2oice(i,k)) (1.0-zfrac(i)) ! zfroi behaves oddly...
314	ztot(i) = zchau(i) + zfroi(i)
315
316	IF (zneb(i).EQ.seuil_neb) ztot(i) = 0.0
317	ztot(i) = MIN(MAX(ztot(i),0.0),zoliq(i))
318	zoliq(i) = MAX(zoliq(i)-ztot(i), 0.0)
319
320	ENDIF
321	ENDDO
322	ENDDO
323
324	!precip scheme modified from Boucher 95
325	CASE(4)
326
327	DO n = 1, ninter
328	DO i = 1, ngrid
329	IF (rneb(i,k).GT.0.0) THEN
330	! this is the ONLY place where zneb and c1 are used
331
332	zchau(i) = ptimestep/FLOAT(ninter) c1 zrho(i) &
333	(zoliq(i)/zneb(i))2reffh2oliq(i,k)zneb(i) zfrac(i)
334	zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
335	zfroi(i) = ptimestep/FLOAT(ninter)/zdz(i)*zoliq(i) &
336	fall2v(reffh2oice(i,k)) (1.0-zfrac(i)) ! zfroi behaves oddly...
337	ztot(i) = zchau(i) + zfroi(i)
338
339	IF (zneb(i).EQ.seuil_neb) ztot(i) = 0.0
340	ztot(i) = MIN(MAX(ztot(i),0.0),zoliq(i))
341	zoliq(i) = MAX(zoliq(i)-ztot(i), 0.0)
342
343	ENDIF
344	ENDDO
345	ENDDO
346
347	END SELECT ! precip_scheme
348
349	! Change in cloud density and surface H2O values
350	DO i = 1, ngrid
351	IF (rneb(i,k).GT.0.0) THEN
352	d_ql(i,k) = (zoliq(i) - ql(i,k))!/ptimestep
353	zrfl(i) = zrfl(i)+ MAX(ql(i,k)-zoliq(i),0.0)*l2c(i,k)/ptimestep
354	ENDIF
355	ENDDO
356
357
358	endif ! if precip_scheme=1
359
360	ENDDO ! of DO k = nlayer, 1, -1
361
362	! Rain or snow on the ground
363	DO i = 1, ngrid
364	if(zrfl(i).lt.0.0)then
365	print*,'Droplets of negative rain are falling...'
366	call abort
367	endif
368	IF (t(i,1) .LT. T_h2O_ice_liq) THEN
369	dqssnow(i) = zrfl(i)
370	dqsrain(i) = 0.0
371	ELSE
372	dqssnow(i) = 0.0
373	dqsrain(i) = zrfl(i) ! liquid water = ice for now
374	ENDIF
375	ENDDO
376
377	! now subroutine -----> GCM variables
378	if (evap_prec) then
379	dqrain(1:ngrid,1:nlayer,i_vap)=d_q(1:ngrid,1:nlayer)/ptimestep
380	d_t(1:ngrid,1:nlayer)=d_t(1:ngrid,1:nlayer)/ptimestep
381	else
382	dqrain(1:ngrid,1:nlayer,i_vap)=0.0
383	d_t(1:ngrid,1:nlayer)=0.0
384	endif
385	dqrain(1:ngrid,1:nlayer,i_ice) = d_ql(1:ngrid,1:nlayer)/ptimestep
386
387	end subroutine rain

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

Download in other formats:

Original Format