Context Navigation

rain.F90 @ 601

Last change on this file since 601 was 526, checked in by jleconte, 13 years ago

13/02/2012 == JL + AS

All outputs are now in netCDF format. Even in 1D (No more G1D)
Clean up of the call to callcorrk when CLFvarying=true
Corrects a bug in writediagspecIR/VI. Output are now in W/m2/cm-1 as a function of the wavenumber in cm-1
Enable writediagspecIR/V to work in the CLFvarying=true case (output now done in Physiq after writediagfi)
Add a simple treatment for the supersaturation of CO2 (see forget et al 2012)
corrects a small bug when no clouds are present in aeropacity

File size: 10.5 KB

Line
1	subroutine rain(ptimestep,pplev,pplay,t,pdt,pq,pdq,d_t,dqrain,dqsrain,dqssnow,rneb)
2
3
4	use watercommon_h, only: To, RLVTT, RCPD, RCPV, RV, RVTMP2
5
6	implicit none
7
8	!==================================================================
9	!
10	! Purpose
11	! -------
12	! Calculates H2O precipitation using simplified microphysics.
13	!
14	! Authors
15	! -------
16	! Adapted from the LMDTERRE code by R. Wordsworth (2009)
17	! Original author Z. X. Li (1993)
18	!
19	!==================================================================
20
21	#include "dimensions.h"
22	#include "dimphys.h"
23	#include "tracer.h"
24	#include "comcstfi.h"
25	#include "callkeys.h"
26
27	! Pre-arguments (for universal model)
28	real pq(ngridmx,nlayermx,nqmx) ! tracer (kg/kg)
29	real qsurf(ngridmx,nqmx) ! tracer at the surface (kg.m-2)
30	REAL pdt(ngridmx,nlayermx),pdq(ngridmx,nlayermx,nqmx)
31
32	real dqrain(ngridmx,nlayermx,nqmx) ! tendency of H2O precipitation (kg/kg.s-1)
33	real dqsrain(ngridmx) ! rain flux at the surface (kg.m-2.s-1)
34	real dqssnow(ngridmx) ! snow flux at the surface (kg.m-2.s-1)
35	REAL d_t(ngridmx,nlayermx) ! temperature increment
36
37	! Arguments
38	REAL ptimestep ! time interval
39	REAL pplev(ngridmx,nlayermx+1) ! inter-layer pressure
40	REAL pplay(ngridmx,nlayermx) ! mid-layer pressure
41	REAL t(ngridmx,nlayermx) ! input temperature (K)
42	REAL zt(ngridmx,nlayermx) ! working temperature (K)
43	REAL ql(ngridmx,nlayermx) ! liquid water (Kg/Kg)
44	REAL q(ngridmx,nlayermx) ! specific humidity (Kg/Kg)
45	REAL rneb(ngridmx,nlayermx) ! cloud fraction
46	REAL d_q(ngridmx,nlayermx) ! water vapor increment
47	REAL d_ql(ngridmx,nlayermx) ! liquid water / ice increment
48
49	! Subroutine options
50	REAL seuil_neb ! Nebulosity threshold
51	PARAMETER (seuil_neb=0.001)
52
53	! REAL ct ! Inverse of cloud precipitation time
54	! PARAMETER (ct=1./1800.)
55	! PARAMETER (ct=1./1849.479)
56
57	REAL cl ! Precipitation threshold
58	PARAMETER (cl=2.0e-4)
59
60	INTEGER ninter
61	PARAMETER (ninter=5)
62
63	logical simple ! Use very simple Emanuel scheme?
64	parameter(simple=.false.)
65
66	logical evap_prec ! Does the rain evaporate?
67	parameter(evap_prec=.true.)
68
69	! for simple scheme
70	real t_crit
71	PARAMETER (t_crit=218.0)
72	real lconvert
73
74	! for precipitation evaporation (old scheme)
75	real eeff1
76	real eeff2
77	! parameter (eeff1=0.95)
78	! parameter (eeff2=0.98)
79	parameter (eeff1=0.5)
80	parameter (eeff2=0.8)
81
82	! Local variables
83	INTEGER i, k, n
84	REAL zqs(ngridmx,nlayermx), zdelta, zcor
85	REAL zrfl(ngridmx), zrfln(ngridmx), zqev, zqevt
86
87	REAL zoliq(ngridmx)
88	REAL ztglace
89	REAL zdz(ngridmx),zrho(ngridmx),ztot(ngridmx), zrhol(ngridmx)
90	REAL zchau(ngridmx),zfroi(ngridmx),zfrac(ngridmx),zneb(ngridmx)
91
92	real ttemp, ptemp
93	real tnext(ngridmx,nlayermx)
94
95	real l2c(ngridmx,nlayermx)
96	real dWtot
97
98
99	! Indices of water vapour and water ice tracers
100	INTEGER, SAVE :: i_vap=0 ! water vapour
101	INTEGER, SAVE :: i_ice=0 ! water ice
102
103	LOGICAL firstcall
104	SAVE firstcall
105
106	! Online functions
107	REAL fallv, zzz ! falling speed of ice crystals
108	fallv (zzz) = 3.29 * ((zzz)**0.16)
109
110	DATA firstcall /.true./
111
112	IF (firstcall) THEN
113
114	i_vap=igcm_h2o_vap
115	i_ice=igcm_h2o_ice
116
117	write(,) "rain: i_ice=",i_ice
118	write(,) " i_vap=",i_vap
119
120	PRINT*, 'in rain.F, ninter=', ninter
121	PRINT*, 'in rain.F, evap_prec=', evap_prec
122
123	!print*,ptimestep
124	!print*,1./ct
125	!if(.not.simple)then
126	! IF (ABS(ptimestep-1./ct).GT.0.001) THEN
127	! PRINT*, 'Must talk to Laurent Li!!!'
128	! call abort
129	! ENDIF
130	!endif
131
132	firstcall = .false.
133	ENDIF
134
135	! GCM -----> subroutine variables
136	DO k = 1, nlayermx
137	DO i = 1, ngridmx
138
139	zt(i,k) = t(i,k)+pdt(i,k)*ptimestep ! a big fat bug was here
140	q(i,k) = pq(i,k,i_vap)+pdq(i,k,i_vap)*ptimestep
141	ql(i,k) = pq(i,k,i_ice)+pdq(i,k,i_ice)*ptimestep
142
143	!q(i,k) = pq(i,k,i_vap)!+pdq(i,k,i_vap)
144	!ql(i,k) = pq(i,k,i_ice)!+pdq(i,k,i_ice)
145
146	if(q(i,k).lt.0.)then ! if this is not done, we don't conserve water
147	q(i,k)=0.
148	endif
149	if(ql(i,k).lt.0.)then
150	ql(i,k)=0.
151	endif
152
153	ENDDO
154	ENDDO
155
156	! Determine the cold clouds by their temperature
157	ztglace = To - 15.0
158
159	! Initialise the outputs
160	DO k = 1, nlayermx
161	DO i = 1, ngridmx
162	d_t(i,k) = 0.0
163	d_q(i,k) = 0.0
164	d_ql(i,k) = 0.0
165	ENDDO
166	ENDDO
167	DO i = 1, ngridmx
168	zrfl(i) = 0.0
169	zrfln(i) = 0.0
170	ENDDO
171
172	! calculate saturation mixing ratio
173	DO k = 1, nlayermx
174	DO i = 1, ngridmx
175	ttemp = zt(i,k)
176	ptemp = pplay(i,k)
177	call watersat(ttemp,ptemp,zqs(i,k))
178	ENDDO
179	ENDDO
180
181	! get column / layer conversion factor
182	DO k = 1, nlayermx
183	DO i = 1, ngridmx
184	!l2c(i,k)=(pplev(i,k)-pplev(i,k+1))/(g*ptimestep)
185	l2c(i,k)=(pplev(i,k)-pplev(i,k+1))/g
186	ENDDO
187	ENDDO
188
189
190	! Vertical loop (from top to bottom)
191	! We carry the rain with us and calculate that added by warm/cold precipitation
192	! processes and that subtracted by evaporation at each level.
193	DO 9999 k = nlayermx, 1, -1
194
195	IF (evap_prec) THEN ! note no rneb dependence!
196	DO i = 1, ngridmx
197	IF (zrfl(i) .GT.0.) THEN
198
199	zqev = MAX (0.0, (zqs(i,k)-q(i,k)))/ptimestep! BC modif here
200	zqevt = 2.0e-5*(1.0-q(i,k)/zqs(i,k)) &
201	sqrt(zrfl(i))l2c(i,k)/pplay(i,k)zt(i,k)R ! BC modif here
202	zqevt = MAX (zqevt, 0.0)
203	zqev = MIN (zqev, zqevt)
204	zqev = MAX (zqev, 0.0)
205	zrfln(i) = zrfl(i) - zqev
206	zrfln(i) = max(zrfln(i),0.0)
207
208	!zqev = MAX (0.0, (zqs(i,k)-q(i,k))*eeff1 )
209	!zqevt = (zrfl(i)/l2c(i,k))*eeff2
210	!zqev = MIN (zqev, zqevt)
211	!zrfln(i) = zrfl(i) - zqev*l2c(i,k)
212	!zrfln(i) = zrfl(i) - 1.5e-5(1.0-q(i,k)/zqs(i,k))sqrt(zrfl(i))
213	!zrfln(i) = min(zrfln(i),zrfl(i))
214	! this is what is actually written in the manual
215
216	d_q(i,k) = - (zrfln(i)-zrfl(i))/l2c(i,k)*ptimestep
217	!d_t(i,k) = d_q(i,k) * RLVTT/RCPD!/(1.0+RVTMP2*q(i,k)) ! double BC modif here
218	d_t(i,k) = - d_q(i,k) * RLVTT/RCPD ! was bugged!
219
220	zrfl(i) = zrfln(i)
221	ENDIF
222	ENDDO
223	ENDIF
224
225	DO i = 1, ngridmx
226	zoliq(i) = 0.0
227	ENDDO
228
229
230	if(simple)then
231
232	DO i = 1, ngridmx
233	ttemp = zt(i,k)
234	IF (ttemp .ge. To) THEN
235	lconvert=rainthreshold
236	ELSEIF (ttemp .gt. t_crit) THEN
237	lconvert=rainthreshold*(1.- t_crit/ttemp)
238	lconvert=MAX(0.0,lconvert)
239	ELSE
240	lconvert=0.
241	ENDIF
242
243
244	IF (ql(i,k).gt.1.e-9) then
245	zneb(i) = MAX(rneb(i,k), seuil_neb)
246	IF ((ql(i,k)/zneb(i)).gt.lconvert)THEN ! precipitate!
247	d_ql(i,k) = -MAX((ql(i,k)-lconvert),0.0)
248	zrfl(i) = zrfl(i) - d_ql(i,k)*l2c(i,k)/ptimestep
249	ENDIF
250	ENDIF
251	ENDDO
252
253	else
254
255	DO i = 1, ngridmx
256	IF (rneb(i,k).GT.0.0) THEN
257	zoliq(i) = ql(i,k)
258	zrho(i) = pplay(i,k) / ( zt(i,k) * R )
259	zdz(i) = (pplev(i,k)-pplev(i,k+1)) / (zrho(i)*g)
260	zfrac(i) = (zt(i,k)-ztglace) / (To-ztglace)
261	zfrac(i) = MAX(zfrac(i), 0.0)
262	zfrac(i) = MIN(zfrac(i), 1.0)
263	zneb(i) = MAX(rneb(i,k), seuil_neb)
264	ENDIF
265	ENDDO
266
267	DO n = 1, ninter
268	DO i = 1, ngridmx
269	IF (rneb(i,k).GT.0.0) THEN
270	zchau(i) = (1./FLOAT(ninter)) * zoliq(i) &
271	* (1.0-EXP(-(zoliq(i)/zneb(i)/cl)*2)) zfrac(i)
272	! warning: this may give dodgy results for physics calls .ne. 48 per day...
273
274	! this is the ONLY place where zneb, ct and cl are used
275	zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
276	zfroi(i) = ptimestep/FLOAT(ninter)/zdz(i)*zoliq(i) &
277	fallv(zrhol(i)) (1.0-zfrac(i)) ! zfroi behaves oddly...
278	! * 0.1 * (1.0-zfrac(i))
279	ztot(i) = zchau(i) + zfroi(i)
280
281	IF (zneb(i).EQ.seuil_neb) ztot(i) = 0.0
282	ztot(i) = MIN(MAX(ztot(i),0.0),zoliq(i))
283	zoliq(i) = MAX(zoliq(i)-ztot(i), 0.0)
284
285	ENDIF
286	ENDDO
287	ENDDO
288
289	! Change in cloud density and surface H2O values
290	DO i = 1, ngridmx
291	IF (rneb(i,k).GT.0.0) THEN
292	d_ql(i,k) = (zoliq(i) - ql(i,k))!/ptimestep
293	zrfl(i) = zrfl(i)+ MAX(ql(i,k)-zoliq(i),0.0)*l2c(i,k)/ptimestep
294	ENDIF
295	ENDDO
296
297	endif ! if simple
298
299	9999 continue
300
301	! Rain or snow on the ground
302	DO i = 1, ngridmx
303	if(zrfl(i).lt.0.0)then
304	print*,'Droplets of negative rain are falling...'
305	call abort
306	endif
307	IF (t(i,1) .LT. To) THEN
308	dqssnow(i) = zrfl(i)
309	dqsrain(i) = 0.0
310	ELSE
311	dqssnow(i) = 0.0
312	dqsrain(i) = zrfl(i) ! liquid water = ice for now
313	ENDIF
314	ENDDO
315
316	! now subroutine -----> GCM variables
317	DO k = 1, nlayermx
318	DO i = 1, ngridmx
319
320	if(evap_prec)then
321	dqrain(i,k,i_vap) = d_q(i,k)/ptimestep
322	d_t(i,k) = d_t(i,k)/ptimestep
323	else
324	dqrain(i,k,i_vap) = 0.0
325	d_t(i,k) = 0.0
326	endif
327	dqrain(i,k,i_ice) = d_ql(i,k)/ptimestep
328
329	ENDDO
330	ENDDO
331
332	RETURN
333	end subroutine rain

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

Download in other formats:

Original Format