Context Navigation

rain.F90 @ 730

Last change on this file since 730 was 728, checked in by jleconte, 13 years ago

18/07/2012 == JL

New water cycle scheme:
- largescale now in F90. Robustness increased by i) including evap inside largescale ii) computing the

condensed water amount iteratively

same improvements in moistadj.
Water thermodynamical data and saturation curves centralized in module watercommn_h
- The saturation curves used are now Tetens formula as they are analyticaly inversible (Ts(P)-> Ps(T)).

New saturation curve yields very good agreement with the former one.

Saturation curves are now generalized for arbitrary water amount (not just q<<1)
The old watersat should be removed soon.

The effect of water vapor on total (surface) pressure can be taken into account by setting

mass_redistrib=.true. in callphys.def (routine mass_redistribution inspired from co2_condense in martian
model but with a different scheme as many routines evaporate/condense water vapor).

New cloud and precipitation scheme (JL + BC):
- The default recovery assumption for computing the total cloud fraction has been changed (total random gave too

large cloud fractions). See totalcloudfrac.F90 for details and to change this.

Totalcloudfraction now set the total cloud fraction to the fraction of the

optically thickest cloud and totalcloudfrac is thus called in aeropacity.

Only the total cloud fraction is used to compute optical depth in aeropacity (no more effective

optical depth with exponential formula).

4 precipitation schemes are now available (see rain.F90 for details). The choice can be made using precip_scheme

in callphys.def. Usage of the more physically based model of Boucher et al 95 (precip_scheme=4) is recommended.
default behavior is set to the former "simple scheme" (precip_scheme=1).

See rain.f90 to determine the parameter to be defined in callphys.def as a function of the precipitation scheme used.

Physiq.F90 now written in a matricial (more F90) way.
Radii (H2O and CO2 cloud particles, aerosols, duts, ...) calculations now centralized in module radii_mod.F90

and work with the new aerosol scheme implemented by Laura K. Some inconsistency may remain in callsedim.

Implementation compiled with ifort and pgf90.
gcm.e runs in Earth and Early Mars case with CO2 and H2O cycle + dust.

File size: 14.0 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, zzz ! falling speed of ice crystals
105	fallv (zzz) = 3.29 * ((zzz)**0.16)
106
107	DATA firstcall /.true./
108
109	IF (firstcall) THEN
110
111	i_vap=igcm_h2o_vap
112	i_ice=igcm_h2o_ice
113
114	write(,) "rain: i_ice=",i_ice
115	write(,) " i_vap=",i_vap
116
117	PRINT*, 'in rain.F, ninter=', ninter
118	PRINT*, 'in rain.F, evap_prec=', evap_prec
119
120	write(,) "Precipitation scheme to use?"
121	precip_scheme=1 ! default value
122	call getin("precip_scheme",precip_scheme)
123	write(,) " precip_scheme = ",precip_scheme
124
125	if (precip_scheme.eq.1) then
126	write(,) "rainthreshold in simple scheme?"
127	rainthreshold=0. ! default value
128	call getin("rainthreshold",rainthreshold)
129	write(,) " rainthreshold = ",rainthreshold
130
131	else if (precip_scheme.eq.2.or.precip_scheme.eq.3) then
132	write(,) "cloud water saturation level in non simple scheme?"
133	cloud_sat=2.6e-4 ! default value
134	call getin("cloud_sat",cloud_sat)
135	write(,) " cloud_sat = ",cloud_sat
136	write(,) "precipitation timescale in non simple scheme?"
137	precip_timescale=3600. ! default value
138	call getin("precip_timescale",precip_timescale)
139	write(,) " precip_timescale = ",precip_timescale
140
141	else if (precip_scheme.eq.4) then
142	write(,) "multiplicative constant in Boucher 95 precip scheme"
143	Cboucher=1. ! default value
144	call getin("Cboucher",Cboucher)
145	write(,) " Cboucher = ",Cboucher
146	c1=1.001.097/rhowaterCboucher*Kboucher
147
148	endif
149
150	firstcall = .false.
151	ENDIF
152
153	! GCM -----> subroutine variables
154	DO k = 1, nlayermx
155	DO i = 1, ngridmx
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, ngridmx
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, ngridmx
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, ngridmx
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, ngridmx
201	!l2c(i,k)=(pplev(i,k)-pplev(i,k+1))/(g*ptimestep)
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	! print*,'in rain',i,k,zrfl(i),zt(i,k),Tsat(i,k)
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	! print*,zqev,zrfl(i),d_q(i,k),d_t(i,k)
223	else
224	zqev = MAX (0.0, (zqs(i,k)-q(i,k)))/ptimestep! there is a bug here
225	zqevt= 2.0e-5*(1.0-q(i,k)/zqs(i,k)) & !default was 2.e-5
226	sqrt(zrfl(i))l2c(i,k)/pplay(i,k)zt(i,k)R ! BC modif here
227	zqevt = MAX (zqevt, 0.0)
228	zqev = MIN (zqev, zqevt)
229	zqev = MAX (zqev, 0.0)
230	zrfln(i)= zrfl(i) - zqev
231	zrfln(i)= max(zrfln(i),0.0)
232
233	d_q(i,k) = - (zrfln(i)-zrfl(i))/l2c(i,k)*ptimestep
234	!d_t(i,k) = d_q(i,k) * RLVTT/RCPD!/(1.0+RVTMP2*q(i,k)) ! double BC modif here
235	d_t(i,k) = - d_q(i,k) * RLVTT/RCPD ! was bugged!
236	zrfl(i) = zrfln(i)
237	end if
238
239
240	ENDIF
241	ENDDO
242	ENDIF
243
244	DO i = 1, ngridmx
245	zoliq(i) = 0.0
246	ENDDO
247
248
249	if(precip_scheme.eq.1)then
250
251	DO i = 1, ngridmx
252	ttemp = zt(i,k)
253	IF (ttemp .ge. T_h2O_ice_liq) THEN
254	lconvert=rainthreshold
255	ELSEIF (ttemp .gt. t_crit) THEN
256	lconvert=rainthreshold*(1.- t_crit/ttemp)
257	lconvert=MAX(0.0,lconvert)
258	ELSE
259	lconvert=0.
260	ENDIF
261
262
263	IF (ql(i,k).gt.1.e-9) then
264	zneb(i) = MAX(rneb(i,k), seuil_neb)
265	IF ((ql(i,k)/zneb(i)).gt.lconvert)THEN ! precipitate!
266	d_ql(i,k) = -MAX((ql(i,k)-lconvert*zneb(i)),0.0)
267	zrfl(i) = zrfl(i) - d_ql(i,k)*l2c(i,k)/ptimestep
268	ENDIF
269	ENDIF
270	ENDDO
271
272	elseif (precip_scheme.ge.2) then
273
274	DO i = 1, ngridmx
275	IF (rneb(i,k).GT.0.0) THEN
276	zoliq(i) = ql(i,k)
277	zrho(i) = pplay(i,k) / ( zt(i,k) * R )
278	zdz(i) = (pplev(i,k)-pplev(i,k+1)) / (zrho(i)*g)
279	zfrac(i) = (zt(i,k)-T_h2O_ice_clouds) / (T_h2O_ice_liq-T_h2O_ice_clouds)
280	zfrac(i) = MAX(zfrac(i), 0.0)
281	zfrac(i) = MIN(zfrac(i), 1.0)
282	zneb(i) = MAX(rneb(i,k), seuil_neb)
283	ENDIF
284	ENDDO
285
286	SELECT CASE(precip_scheme)
287	!precip scheme from Sundquist 78
288	CASE(2)
289
290	DO n = 1, ninter
291	DO i = 1, ngridmx
292	IF (rneb(i,k).GT.0.0) THEN
293	! this is the ONLY place where zneb, precip_timescale and cloud_sat are used
294
295	zchau(i) = (ptimestep/(FLOAT(ninter)precip_timescale)) zoliq(i) &
296	* (1.0-EXP(-(zoliq(i)/zneb(i)/cloud_sat)*2)) zfrac(i)
297	zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
298	zfroi(i) = ptimestep/FLOAT(ninter)/zdz(i)*zoliq(i) &
299	fallv(zrhol(i)) (1.0-zfrac(i)) ! zfroi behaves oddly...
300	! * 0.1 * (1.0-zfrac(i))
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, ngridmx
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	fallv(zrhol(i)) (1.0-zfrac(i)) ! zfroi behaves oddly...
323	! * 0.1 * (1.0-zfrac(i))
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	!recalculate liquid water particle radii
338	call h2o_cloudrad(ql,reffh2oliq,reffh2oice)
339
340	DO n = 1, ninter
341	DO i = 1, ngridmx
342	IF (rneb(i,k).GT.0.0) THEN
343	! this is the ONLY place where zneb and c1 are used
344
345	zchau(i) = ptimestep/FLOAT(ninter) c1 zrho(i) &
346	(zoliq(i)/zneb(i))2reffh2oliq(i,k)zneb(i) zfrac(i)
347	zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
348	zfroi(i) = ptimestep/FLOAT(ninter)/zdz(i)*zoliq(i) &
349	fallv(zrhol(i)) (1.0-zfrac(i)) ! zfroi behaves oddly...
350	! * 0.1 * (1.0-zfrac(i))
351	ztot(i) = zchau(i) + zfroi(i)
352
353	IF (zneb(i).EQ.seuil_neb) ztot(i) = 0.0
354	ztot(i) = MIN(MAX(ztot(i),0.0),zoliq(i))
355	zoliq(i) = MAX(zoliq(i)-ztot(i), 0.0)
356
357	ENDIF
358	ENDDO
359	ENDDO
360
361	END SELECT ! precip_scheme
362
363	! Change in cloud density and surface H2O values
364	DO i = 1, ngridmx
365	IF (rneb(i,k).GT.0.0) THEN
366	d_ql(i,k) = (zoliq(i) - ql(i,k))!/ptimestep
367	zrfl(i) = zrfl(i)+ MAX(ql(i,k)-zoliq(i),0.0)*l2c(i,k)/ptimestep
368	ENDIF
369	ENDDO
370
371
372	endif ! if precip_scheme=1
373
374	9999 continue
375
376	! Rain or snow on the ground
377	DO i = 1, ngridmx
378	if(zrfl(i).lt.0.0)then
379	print*,'Droplets of negative rain are falling...'
380	call abort
381	endif
382	IF (t(i,1) .LT. T_h2O_ice_liq) THEN
383	dqssnow(i) = zrfl(i)
384	dqsrain(i) = 0.0
385	ELSE
386	dqssnow(i) = 0.0
387	dqsrain(i) = zrfl(i) ! liquid water = ice for now
388	ENDIF
389	ENDDO
390
391	! now subroutine -----> GCM variables
392	DO k = 1, nlayermx
393	DO i = 1, ngridmx
394
395	if(evap_prec)then
396	dqrain(i,k,i_vap) = d_q(i,k)/ptimestep
397	d_t(i,k) = d_t(i,k)/ptimestep
398	else
399	dqrain(i,k,i_vap) = 0.0
400	d_t(i,k) = 0.0
401	endif
402	dqrain(i,k,i_ice) = d_ql(i,k)/ptimestep
403
404	ENDDO
405	ENDDO
406
407	RETURN
408	end subroutine rain

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

Download in other formats:

Original Format