Context Navigation

rain.F90 @ 135

Last change on this file since 135 was 135, checked in by aslmd, 14 years ago
CHANGEMENT ARBORESCENCE ETAPE 2 -- NON COMPLET
File size: 9.5 KB

Line
1
2	subroutine rain(ptimestep,pplev,pplay,t,pdt,pq,pdq,d_t,dqrain,dqsrain,dqssnow,rneb)
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) ! temperature (K)
42	REAL ql(ngridmx,nlayermx) ! liquid water (Kg/Kg)
43	REAL q(ngridmx,nlayermx) ! specific humidity (Kg/Kg)
44	REAL rneb(ngridmx,nlayermx) ! cloud fraction
45	REAL d_q(ngridmx,nlayermx) ! water vapor increment
46	REAL d_ql(ngridmx,nlayermx) ! liquid water / ice increment
47
48	! Subroutine options
49	REAL seuil_neb ! Nebulosity threshold
50	PARAMETER (seuil_neb=0.001)
51
52	REAL ct ! Inverse of cloud precipitation time
53	! PARAMETER (ct=1./1800.)
54	PARAMETER (ct=1./1849.479)
55
56	REAL cl ! Precipitation threshold
57	PARAMETER (cl=2.0e-4)
58
59	INTEGER ninter
60	PARAMETER (ninter=5)
61
62	logical simple ! Use very simple Emanuel scheme?
63	parameter(simple=.true.)
64
65	logical evap_prec ! Does the rain evaporate?
66	parameter(evap_prec=.true.)
67
68	! for simple scheme
69	real t_crit
70	PARAMETER (t_crit=218.0)
71	real lconvert
72
73	! for precipitation evaporation (old scheme)
74	real eeff1
75	real eeff2
76	! parameter (eeff1=0.95)
77	! parameter (eeff2=0.98)
78	parameter (eeff1=0.5)
79	parameter (eeff2=0.8)
80
81	! Local variables
82	INTEGER i, k, n
83	REAL zqs(ngridmx,nlayermx), zdelta, zcor
84	REAL zrfl(ngridmx), zrfln(ngridmx), zqev, zqevt
85
86	REAL zoliq(ngridmx)
87	REAL ztglace
88	REAL zdz(ngridmx),zrho(ngridmx),ztot(ngridmx), zrhol(ngridmx)
89	REAL zchau(ngridmx),zfroi(ngridmx),zfrac(ngridmx),zneb(ngridmx)
90
91	real ttemp, ptemp
92	real tnext(ngridmx,nlayermx)
93
94	REAL l2c(ngridmx,nlayermx)
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	print*,ptimestep
121	print*,1./ct
122
123	if(.not.simple)then
124	IF (ABS(ptimestep-1./ct).GT.0.001) THEN
125	PRINT*, 'Must talk to Laurent Li!!!'
126	call abort
127	ENDIF
128	endif
129
130	firstcall = .false.
131	ENDIF
132
133	! GCM -----> subroutine variables
134	DO k = 1, nlayermx
135	DO i = 1, ngridmx
136
137	q(i,k) = pq(i,k,i_vap)!+pdq(i,k,i_vap)
138	ql(i,k) = pq(i,k,i_ice)!+pdq(i,k,i_ice)
139
140	if(q(i,k).lt.0.)then ! if this is not done, we don't conserve water
141	q(i,k)=0.
142	endif
143	if(ql(i,k).lt.0.)then
144	ql(i,k)=0.
145	endif
146
147	ENDDO
148	ENDDO
149
150	! Determine the cold clouds by their temperature
151	ztglace = To - 15.0
152
153	! Initialise the outputs
154	DO k = 1, nlayermx
155	DO i = 1, ngridmx
156	d_t(i,k) = 0.0
157	d_q(i,k) = 0.0
158	d_ql(i,k) = 0.0
159	ENDDO
160	ENDDO
161	DO i = 1, ngridmx
162	zrfl(i) = 0.0
163	zrfln(i) = 0.0
164	ENDDO
165
166	! calculate saturation mixing ratio
167	DO k = 1, nlayermx
168	DO i = 1, ngridmx
169	ttemp = t(i,k)
170	ptemp = pplay(i,k)
171	call watersat_2(ttemp,ptemp,zqs(i,k))
172	ENDDO
173	ENDDO
174
175	! get column / layer conversion factor (+ptimstep)
176	DO k = 1, nlayermx
177	DO i = 1, ngridmx
178	l2c(i,k)=(pplev(i,k)-pplev(i,k+1))/(g*ptimestep)
179	ENDDO
180	ENDDO
181
182
183	! Vertical loop (from top to bottom)
184	! We carry the rain with us and calculate that added by warm/cold precipitation
185	! processes and that subtracted by evaporation at each level.
186	DO 9999 k = nlayermx, 1, -1
187
188	IF (evap_prec) THEN ! note no rneb dependence!
189	DO i = 1, ngridmx
190	IF (zrfl(i) .GT.0.) THEN
191	!zqev = MAX (0.0, (zqs(i,k)-q(i,k))*eeff1 )
192	!zqevt = (zrfl(i)/l2c(i,k))*eeff2
193	!zqev = MIN (zqev, zqevt)
194	!zrfln(i) = zrfl(i) - zqev*l2c(i,k)
195
196	zrfln(i) = zrfl(i) - 1.5e-5(1.0-q(i,k)/zqs(i,k))sqrt(zrfl(i))
197	zrfln(i) = min(zrfln(i),zrfl(i))
198	! this is what is actually written in the manual
199
200	d_q(i,k) = - (zrfln(i)-zrfl(i))/l2c(i,k)
201	d_t(i,k) = d_q(i,k) * RLVTT/RCPD/(1.0+RVTMP2*q(i,k))
202
203	zrfl(i) = zrfln(i)
204	ENDIF
205	ENDDO
206	ENDIF
207
208	DO i = 1, ngridmx
209	zoliq(i) = 0.0
210	ENDDO
211
212
213	if(simple)then
214
215	DO i = 1, ngridmx
216	ttemp = t(i,k)
217	IF (ttemp .ge. To) THEN
218	lconvert=rainthreshold
219	ELSEIF (ttemp .gt. t_crit) THEN
220	lconvert=rainthreshold*(1.- t_crit/ttemp)
221	ELSE
222	lconvert=0.
223	ENDIF
224
225	IF (ql(i,k).gt.lconvert)THEN ! precipitate!
226	d_ql(i,k) = (lconvert-ql(i,k))/ptimestep
227	zrfl(i) = zrfl(i) + max(ql(i,k) - lconvert,0.0)*l2c(i,k)
228	ENDIF
229	ENDDO
230
231	else
232
233	DO i = 1, ngridmx
234	IF (rneb(i,k).GT.0.0) THEN
235	zoliq(i) = ql(i,k)
236	zrho(i) = pplay(i,k) / ( t(i,k) * R )
237	zdz(i) = (pplev(i,k)-pplev(i,k+1)) / (zrho(i)*g)
238	zfrac(i) = (t(i,k)-ztglace) / (To-ztglace)
239	zfrac(i) = MAX(zfrac(i), 0.0)
240	zfrac(i) = MIN(zfrac(i), 1.0)
241	zneb(i) = MAX(rneb(i,k), seuil_neb)
242	ENDIF
243	ENDDO
244
245	DO n = 1, ninter
246	DO i = 1, ngridmx
247	IF (rneb(i,k).GT.0.0) THEN
248	zchau(i) = ctptimestep/FLOAT(ninter) zoliq(i) &
249	* (1.0-EXP(-(zoliq(i)/zneb(i)/cl)*2)) zfrac(i)
250	! this is the ONLY place where zneb, ct and cl are used
251	zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
252	zfroi(i) = ptimestep/FLOAT(ninter)/zdz(i)*zoliq(i) &
253	fallv(zrhol(i)) (1.0-zfrac(i)) ! zfroi behaves oddly...
254	! * 0.1 * (1.0-zfrac(i))
255	ztot(i) = zchau(i) + zfroi(i)
256
257	IF (zneb(i).EQ.seuil_neb) ztot(i) = 0.0
258	ztot(i) = MIN(MAX(ztot(i),0.0),zoliq(i))
259	zoliq(i) = MAX(zoliq(i)-ztot(i), 0.0)
260
261	ENDIF
262	ENDDO
263	ENDDO
264
265	! Change in cloud density and surface H2O values
266	DO i = 1, ngridmx
267	IF (rneb(i,k).GT.0.0) THEN
268	d_ql(i,k) = (zoliq(i) - ql(i,k))/ptimestep
269	zrfl(i) = zrfl(i)+ MAX(ql(i,k)-zoliq(i),0.0)*l2c(i,k)
270	ENDIF
271	ENDDO
272
273	endif ! if simple
274
275	9999 continue
276
277	! Rain or snow on the ground
278	DO i = 1, ngridmx
279	IF (t(i,1) .LT. To) THEN
280	dqssnow(i) = zrfl(i)
281	ELSE
282	dqsrain(i) = zrfl(i) ! liquid water = ice for now
283	ENDIF
284	ENDDO
285
286	! now subroutine -----> GCM variables
287	DO k = 1, nlayermx
288	DO i = 1, ngridmx
289
290	if(evap_prec)then
291	dqrain(i,k,i_vap) = d_q(i,k)/ptimestep
292	d_t(i,k) = d_t(i,k)/ptimestep
293	else
294	dqrain(i,k,i_vap) = 0.0
295	d_t(i,k) = 0.0
296	endif
297	dqrain(i,k,i_ice) = d_ql(i,k)
298
299	ENDDO
300	ENDDO
301
302	! debugging
303	! print*,'zrfl=',zrfl
304	! print*,'dqrain=',dqrain
305	! print*,'q=',q
306	! print*,'ql=',ql
307	! print*,'dql=',d_ql
308	! DO k = 1, nlayermx
309	! DO i = 1, ngridmx
310	! if(ql(i,k).lt.0.0) then
311	! print*,'below zero!!!!'
312	! call abort
313	! endif
314	! ENDDO
315	! ENDDO
316
317
318	RETURN
319	end subroutine rain

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

Download in other formats:

Original Format