1 |
|
|
|
2 |
|
|
! $Id: cv_routines.F90 2311 2015-06-25 07:45:24Z emillour $ |
3 |
|
|
|
4 |
|
|
SUBROUTINE cv_param(nd) |
5 |
|
|
IMPLICIT NONE |
6 |
|
|
|
7 |
|
|
! ------------------------------------------------------------ |
8 |
|
|
! Set parameters for convectL |
9 |
|
|
! (includes microphysical parameters and parameters that |
10 |
|
|
! control the rate of approach to quasi-equilibrium) |
11 |
|
|
! ------------------------------------------------------------ |
12 |
|
|
|
13 |
|
|
! *** ELCRIT IS THE AUTOCONVERSION THERSHOLD WATER CONTENT (gm/gm) *** |
14 |
|
|
! *** TLCRIT IS CRITICAL TEMPERATURE BELOW WHICH THE AUTO- *** |
15 |
|
|
! *** CONVERSION THRESHOLD IS ASSUMED TO BE ZERO *** |
16 |
|
|
! *** (THE AUTOCONVERSION THRESHOLD VARIES LINEARLY *** |
17 |
|
|
! *** BETWEEN 0 C AND TLCRIT) *** |
18 |
|
|
! *** ENTP IS THE COEFFICIENT OF MIXING IN THE ENTRAINMENT *** |
19 |
|
|
! *** FORMULATION *** |
20 |
|
|
! *** SIGD IS THE FRACTIONAL AREA COVERED BY UNSATURATED DNDRAFT *** |
21 |
|
|
! *** SIGS IS THE FRACTION OF PRECIPITATION FALLING OUTSIDE *** |
22 |
|
|
! *** OF CLOUD *** |
23 |
|
|
! *** OMTRAIN IS THE ASSUMED FALL SPEED (P/s) OF RAIN *** |
24 |
|
|
! *** OMTSNOW IS THE ASSUMED FALL SPEED (P/s) OF SNOW *** |
25 |
|
|
! *** COEFFR IS A COEFFICIENT GOVERNING THE RATE OF EVAPORATION *** |
26 |
|
|
! *** OF RAIN *** |
27 |
|
|
! *** COEFFS IS A COEFFICIENT GOVERNING THE RATE OF EVAPORATION *** |
28 |
|
|
! *** OF SNOW *** |
29 |
|
|
! *** CU IS THE COEFFICIENT GOVERNING CONVECTIVE MOMENTUM *** |
30 |
|
|
! *** TRANSPORT *** |
31 |
|
|
! *** DTMAX IS THE MAXIMUM NEGATIVE TEMPERATURE PERTURBATION *** |
32 |
|
|
! *** A LIFTED PARCEL IS ALLOWED TO HAVE BELOW ITS LFC *** |
33 |
|
|
! *** ALPHA AND DAMP ARE PARAMETERS THAT CONTROL THE RATE OF *** |
34 |
|
|
! *** APPROACH TO QUASI-EQUILIBRIUM *** |
35 |
|
|
! *** (THEIR STANDARD VALUES ARE 0.20 AND 0.1, RESPECTIVELY) *** |
36 |
|
|
! *** (DAMP MUST BE LESS THAN 1) *** |
37 |
|
|
|
38 |
|
|
include "cvparam.h" |
39 |
|
|
INTEGER nd |
40 |
|
|
CHARACTER (LEN=20) :: modname = 'cv_routines' |
41 |
|
|
CHARACTER (LEN=80) :: abort_message |
42 |
|
|
|
43 |
|
|
! noff: integer limit for convection (nd-noff) |
44 |
|
|
! minorig: First level of convection |
45 |
|
|
|
46 |
|
|
noff = 2 |
47 |
|
|
minorig = 2 |
48 |
|
|
|
49 |
|
|
nl = nd - noff |
50 |
|
|
nlp = nl + 1 |
51 |
|
|
nlm = nl - 1 |
52 |
|
|
|
53 |
|
|
elcrit = 0.0011 |
54 |
|
|
tlcrit = -55.0 |
55 |
|
|
entp = 1.5 |
56 |
|
|
sigs = 0.12 |
57 |
|
|
sigd = 0.05 |
58 |
|
|
omtrain = 50.0 |
59 |
|
|
omtsnow = 5.5 |
60 |
|
|
coeffr = 1.0 |
61 |
|
|
coeffs = 0.8 |
62 |
|
|
dtmax = 0.9 |
63 |
|
|
|
64 |
|
|
cu = 0.70 |
65 |
|
|
|
66 |
|
|
betad = 10.0 |
67 |
|
|
|
68 |
|
|
damp = 0.1 |
69 |
|
|
alpha = 0.2 |
70 |
|
|
|
71 |
|
|
delta = 0.01 ! cld |
72 |
|
|
|
73 |
|
|
RETURN |
74 |
|
|
END SUBROUTINE cv_param |
75 |
|
|
|
76 |
|
|
SUBROUTINE cv_prelim(len, nd, ndp1, t, q, p, ph, lv, cpn, tv, gz, h, hm) |
77 |
|
|
IMPLICIT NONE |
78 |
|
|
|
79 |
|
|
! ===================================================================== |
80 |
|
|
! --- CALCULATE ARRAYS OF GEOPOTENTIAL, HEAT CAPACITY & STATIC ENERGY |
81 |
|
|
! ===================================================================== |
82 |
|
|
|
83 |
|
|
! inputs: |
84 |
|
|
INTEGER len, nd, ndp1 |
85 |
|
|
REAL t(len, nd), q(len, nd), p(len, nd), ph(len, ndp1) |
86 |
|
|
|
87 |
|
|
! outputs: |
88 |
|
|
REAL lv(len, nd), cpn(len, nd), tv(len, nd) |
89 |
|
|
REAL gz(len, nd), h(len, nd), hm(len, nd) |
90 |
|
|
|
91 |
|
|
! local variables: |
92 |
|
|
INTEGER k, i |
93 |
|
|
REAL cpx(len, nd) |
94 |
|
|
|
95 |
|
|
include "cvthermo.h" |
96 |
|
|
include "cvparam.h" |
97 |
|
|
|
98 |
|
|
|
99 |
|
|
DO k = 1, nlp |
100 |
|
|
DO i = 1, len |
101 |
|
|
lv(i, k) = lv0 - clmcpv*(t(i,k)-t0) |
102 |
|
|
cpn(i, k) = cpd*(1.0-q(i,k)) + cpv*q(i, k) |
103 |
|
|
cpx(i, k) = cpd*(1.0-q(i,k)) + cl*q(i, k) |
104 |
|
|
tv(i, k) = t(i, k)*(1.0+q(i,k)*epsim1) |
105 |
|
|
END DO |
106 |
|
|
END DO |
107 |
|
|
|
108 |
|
|
! gz = phi at the full levels (same as p). |
109 |
|
|
|
110 |
|
|
DO i = 1, len |
111 |
|
|
gz(i, 1) = 0.0 |
112 |
|
|
END DO |
113 |
|
|
DO k = 2, nlp |
114 |
|
|
DO i = 1, len |
115 |
|
|
gz(i, k) = gz(i, k-1) + hrd*(tv(i,k-1)+tv(i,k))*(p(i,k-1)-p(i,k))/ph(i, & |
116 |
|
|
k) |
117 |
|
|
END DO |
118 |
|
|
END DO |
119 |
|
|
|
120 |
|
|
! h = phi + cpT (dry static energy). |
121 |
|
|
! hm = phi + cp(T-Tbase)+Lq |
122 |
|
|
|
123 |
|
|
DO k = 1, nlp |
124 |
|
|
DO i = 1, len |
125 |
|
|
h(i, k) = gz(i, k) + cpn(i, k)*t(i, k) |
126 |
|
|
hm(i, k) = gz(i, k) + cpx(i, k)*(t(i,k)-t(i,1)) + lv(i, k)*q(i, k) |
127 |
|
|
END DO |
128 |
|
|
END DO |
129 |
|
|
|
130 |
|
|
RETURN |
131 |
|
|
END SUBROUTINE cv_prelim |
132 |
|
|
|
133 |
|
|
SUBROUTINE cv_feed(len, nd, t, q, qs, p, hm, gz, nk, icb, icbmax, iflag, tnk, & |
134 |
|
|
qnk, gznk, plcl) |
135 |
|
|
IMPLICIT NONE |
136 |
|
|
|
137 |
|
|
! ================================================================ |
138 |
|
|
! Purpose: CONVECTIVE FEED |
139 |
|
|
! ================================================================ |
140 |
|
|
|
141 |
|
|
include "cvparam.h" |
142 |
|
|
|
143 |
|
|
! inputs: |
144 |
|
|
INTEGER len, nd |
145 |
|
|
REAL t(len, nd), q(len, nd), qs(len, nd), p(len, nd) |
146 |
|
|
REAL hm(len, nd), gz(len, nd) |
147 |
|
|
|
148 |
|
|
! outputs: |
149 |
|
|
INTEGER iflag(len), nk(len), icb(len), icbmax |
150 |
|
|
REAL tnk(len), qnk(len), gznk(len), plcl(len) |
151 |
|
|
|
152 |
|
|
! local variables: |
153 |
|
|
INTEGER i, k |
154 |
|
|
INTEGER ihmin(len) |
155 |
|
|
REAL work(len) |
156 |
|
|
REAL pnk(len), qsnk(len), rh(len), chi(len) |
157 |
|
|
|
158 |
|
|
! ------------------------------------------------------------------- |
159 |
|
|
! --- Find level of minimum moist static energy |
160 |
|
|
! --- If level of minimum moist static energy coincides with |
161 |
|
|
! --- or is lower than minimum allowable parcel origin level, |
162 |
|
|
! --- set iflag to 6. |
163 |
|
|
! ------------------------------------------------------------------- |
164 |
|
|
|
165 |
|
|
DO i = 1, len |
166 |
|
|
work(i) = 1.0E12 |
167 |
|
|
ihmin(i) = nl |
168 |
|
|
END DO |
169 |
|
|
DO k = 2, nlp |
170 |
|
|
DO i = 1, len |
171 |
|
|
IF ((hm(i,k)<work(i)) .AND. (hm(i,k)<hm(i,k-1))) THEN |
172 |
|
|
work(i) = hm(i, k) |
173 |
|
|
ihmin(i) = k |
174 |
|
|
END IF |
175 |
|
|
END DO |
176 |
|
|
END DO |
177 |
|
|
DO i = 1, len |
178 |
|
|
ihmin(i) = min(ihmin(i), nlm) |
179 |
|
|
IF (ihmin(i)<=minorig) THEN |
180 |
|
|
iflag(i) = 6 |
181 |
|
|
END IF |
182 |
|
|
END DO |
183 |
|
|
|
184 |
|
|
! ------------------------------------------------------------------- |
185 |
|
|
! --- Find that model level below the level of minimum moist static |
186 |
|
|
! --- energy that has the maximum value of moist static energy |
187 |
|
|
! ------------------------------------------------------------------- |
188 |
|
|
|
189 |
|
|
DO i = 1, len |
190 |
|
|
work(i) = hm(i, minorig) |
191 |
|
|
nk(i) = minorig |
192 |
|
|
END DO |
193 |
|
|
DO k = minorig + 1, nl |
194 |
|
|
DO i = 1, len |
195 |
|
|
IF ((hm(i,k)>work(i)) .AND. (k<=ihmin(i))) THEN |
196 |
|
|
work(i) = hm(i, k) |
197 |
|
|
nk(i) = k |
198 |
|
|
END IF |
199 |
|
|
END DO |
200 |
|
|
END DO |
201 |
|
|
! ------------------------------------------------------------------- |
202 |
|
|
! --- Check whether parcel level temperature and specific humidity |
203 |
|
|
! --- are reasonable |
204 |
|
|
! ------------------------------------------------------------------- |
205 |
|
|
DO i = 1, len |
206 |
|
|
IF (((t(i,nk(i))<250.0) .OR. (q(i,nk(i))<=0.0) .OR. (p(i,ihmin(i))< & |
207 |
|
|
400.0)) .AND. (iflag(i)==0)) iflag(i) = 7 |
208 |
|
|
END DO |
209 |
|
|
! ------------------------------------------------------------------- |
210 |
|
|
! --- Calculate lifted condensation level of air at parcel origin level |
211 |
|
|
! --- (Within 0.2% of formula of Bolton, MON. WEA. REV.,1980) |
212 |
|
|
! ------------------------------------------------------------------- |
213 |
|
|
DO i = 1, len |
214 |
|
|
tnk(i) = t(i, nk(i)) |
215 |
|
|
qnk(i) = q(i, nk(i)) |
216 |
|
|
gznk(i) = gz(i, nk(i)) |
217 |
|
|
pnk(i) = p(i, nk(i)) |
218 |
|
|
qsnk(i) = qs(i, nk(i)) |
219 |
|
|
|
220 |
|
|
rh(i) = qnk(i)/qsnk(i) |
221 |
|
|
rh(i) = min(1.0, rh(i)) |
222 |
|
|
chi(i) = tnk(i)/(1669.0-122.0*rh(i)-tnk(i)) |
223 |
|
|
plcl(i) = pnk(i)*(rh(i)**chi(i)) |
224 |
|
|
IF (((plcl(i)<200.0) .OR. (plcl(i)>=2000.0)) .AND. (iflag(i)==0)) iflag(i & |
225 |
|
|
) = 8 |
226 |
|
|
END DO |
227 |
|
|
! ------------------------------------------------------------------- |
228 |
|
|
! --- Calculate first level above lcl (=icb) |
229 |
|
|
! ------------------------------------------------------------------- |
230 |
|
|
DO i = 1, len |
231 |
|
|
icb(i) = nlm |
232 |
|
|
END DO |
233 |
|
|
|
234 |
|
|
DO k = minorig, nl |
235 |
|
|
DO i = 1, len |
236 |
|
|
IF ((k>=(nk(i)+1)) .AND. (p(i,k)<plcl(i))) icb(i) = min(icb(i), k) |
237 |
|
|
END DO |
238 |
|
|
END DO |
239 |
|
|
|
240 |
|
|
DO i = 1, len |
241 |
|
|
IF ((icb(i)>=nlm) .AND. (iflag(i)==0)) iflag(i) = 9 |
242 |
|
|
END DO |
243 |
|
|
|
244 |
|
|
! Compute icbmax. |
245 |
|
|
|
246 |
|
|
icbmax = 2 |
247 |
|
|
DO i = 1, len |
248 |
|
|
icbmax = max(icbmax, icb(i)) |
249 |
|
|
END DO |
250 |
|
|
|
251 |
|
|
RETURN |
252 |
|
|
END SUBROUTINE cv_feed |
253 |
|
|
|
254 |
|
|
SUBROUTINE cv_undilute1(len, nd, t, q, qs, gz, p, nk, icb, icbmax, tp, tvp, & |
255 |
|
|
clw) |
256 |
|
|
IMPLICIT NONE |
257 |
|
|
|
258 |
|
|
include "cvthermo.h" |
259 |
|
|
include "cvparam.h" |
260 |
|
|
|
261 |
|
|
! inputs: |
262 |
|
|
INTEGER len, nd |
263 |
|
|
INTEGER nk(len), icb(len), icbmax |
264 |
|
|
REAL t(len, nd), q(len, nd), qs(len, nd), gz(len, nd) |
265 |
|
|
REAL p(len, nd) |
266 |
|
|
|
267 |
|
|
! outputs: |
268 |
|
|
REAL tp(len, nd), tvp(len, nd), clw(len, nd) |
269 |
|
|
|
270 |
|
|
! local variables: |
271 |
|
|
INTEGER i, k |
272 |
|
|
REAL tg, qg, alv, s, ahg, tc, denom, es, rg |
273 |
|
|
REAL ah0(len), cpp(len) |
274 |
|
|
REAL tnk(len), qnk(len), gznk(len), ticb(len), gzicb(len) |
275 |
|
|
|
276 |
|
|
! ------------------------------------------------------------------- |
277 |
|
|
! --- Calculates the lifted parcel virtual temperature at nk, |
278 |
|
|
! --- the actual temperature, and the adiabatic |
279 |
|
|
! --- liquid water content. The procedure is to solve the equation. |
280 |
|
|
! cp*tp+L*qp+phi=cp*tnk+L*qnk+gznk. |
281 |
|
|
! ------------------------------------------------------------------- |
282 |
|
|
|
283 |
|
|
DO i = 1, len |
284 |
|
|
tnk(i) = t(i, nk(i)) |
285 |
|
|
qnk(i) = q(i, nk(i)) |
286 |
|
|
gznk(i) = gz(i, nk(i)) |
287 |
|
|
ticb(i) = t(i, icb(i)) |
288 |
|
|
gzicb(i) = gz(i, icb(i)) |
289 |
|
|
END DO |
290 |
|
|
|
291 |
|
|
! *** Calculate certain parcel quantities, including static energy *** |
292 |
|
|
|
293 |
|
|
DO i = 1, len |
294 |
|
|
ah0(i) = (cpd*(1.-qnk(i))+cl*qnk(i))*tnk(i) + qnk(i)*(lv0-clmcpv*(tnk(i)- & |
295 |
|
|
273.15)) + gznk(i) |
296 |
|
|
cpp(i) = cpd*(1.-qnk(i)) + qnk(i)*cpv |
297 |
|
|
END DO |
298 |
|
|
|
299 |
|
|
! *** Calculate lifted parcel quantities below cloud base *** |
300 |
|
|
|
301 |
|
|
DO k = minorig, icbmax - 1 |
302 |
|
|
DO i = 1, len |
303 |
|
|
tp(i, k) = tnk(i) - (gz(i,k)-gznk(i))/cpp(i) |
304 |
|
|
tvp(i, k) = tp(i, k)*(1.+qnk(i)*epsi) |
305 |
|
|
END DO |
306 |
|
|
END DO |
307 |
|
|
|
308 |
|
|
! *** Find lifted parcel quantities above cloud base *** |
309 |
|
|
|
310 |
|
|
DO i = 1, len |
311 |
|
|
tg = ticb(i) |
312 |
|
|
qg = qs(i, icb(i)) |
313 |
|
|
alv = lv0 - clmcpv*(ticb(i)-t0) |
314 |
|
|
|
315 |
|
|
! First iteration. |
316 |
|
|
|
317 |
|
|
s = cpd + alv*alv*qg/(rrv*ticb(i)*ticb(i)) |
318 |
|
|
s = 1./s |
319 |
|
|
ahg = cpd*tg + (cl-cpd)*qnk(i)*ticb(i) + alv*qg + gzicb(i) |
320 |
|
|
tg = tg + s*(ah0(i)-ahg) |
321 |
|
|
tg = max(tg, 35.0) |
322 |
|
|
tc = tg - t0 |
323 |
|
|
denom = 243.5 + tc |
324 |
|
|
IF (tc>=0.0) THEN |
325 |
|
|
es = 6.112*exp(17.67*tc/denom) |
326 |
|
|
ELSE |
327 |
|
|
es = exp(23.33086-6111.72784/tg+0.15215*log(tg)) |
328 |
|
|
END IF |
329 |
|
|
qg = eps*es/(p(i,icb(i))-es*(1.-eps)) |
330 |
|
|
|
331 |
|
|
! Second iteration. |
332 |
|
|
|
333 |
|
|
s = cpd + alv*alv*qg/(rrv*ticb(i)*ticb(i)) |
334 |
|
|
s = 1./s |
335 |
|
|
ahg = cpd*tg + (cl-cpd)*qnk(i)*ticb(i) + alv*qg + gzicb(i) |
336 |
|
|
tg = tg + s*(ah0(i)-ahg) |
337 |
|
|
tg = max(tg, 35.0) |
338 |
|
|
tc = tg - t0 |
339 |
|
|
denom = 243.5 + tc |
340 |
|
|
IF (tc>=0.0) THEN |
341 |
|
|
es = 6.112*exp(17.67*tc/denom) |
342 |
|
|
ELSE |
343 |
|
|
es = exp(23.33086-6111.72784/tg+0.15215*log(tg)) |
344 |
|
|
END IF |
345 |
|
|
qg = eps*es/(p(i,icb(i))-es*(1.-eps)) |
346 |
|
|
|
347 |
|
|
alv = lv0 - clmcpv*(ticb(i)-273.15) |
348 |
|
|
tp(i, icb(i)) = (ah0(i)-(cl-cpd)*qnk(i)*ticb(i)-gz(i,icb(i))-alv*qg)/cpd |
349 |
|
|
clw(i, icb(i)) = qnk(i) - qg |
350 |
|
|
clw(i, icb(i)) = max(0.0, clw(i,icb(i))) |
351 |
|
|
rg = qg/(1.-qnk(i)) |
352 |
|
|
tvp(i, icb(i)) = tp(i, icb(i))*(1.+rg*epsi) |
353 |
|
|
END DO |
354 |
|
|
|
355 |
|
|
DO k = minorig, icbmax |
356 |
|
|
DO i = 1, len |
357 |
|
|
tvp(i, k) = tvp(i, k) - tp(i, k)*qnk(i) |
358 |
|
|
END DO |
359 |
|
|
END DO |
360 |
|
|
|
361 |
|
|
RETURN |
362 |
|
|
END SUBROUTINE cv_undilute1 |
363 |
|
|
|
364 |
|
|
SUBROUTINE cv_trigger(len, nd, icb, cbmf, tv, tvp, iflag) |
365 |
|
|
IMPLICIT NONE |
366 |
|
|
|
367 |
|
|
! ------------------------------------------------------------------- |
368 |
|
|
! --- Test for instability. |
369 |
|
|
! --- If there was no convection at last time step and parcel |
370 |
|
|
! --- is stable at icb, then set iflag to 4. |
371 |
|
|
! ------------------------------------------------------------------- |
372 |
|
|
|
373 |
|
|
include "cvparam.h" |
374 |
|
|
|
375 |
|
|
! inputs: |
376 |
|
|
INTEGER len, nd, icb(len) |
377 |
|
|
REAL cbmf(len), tv(len, nd), tvp(len, nd) |
378 |
|
|
|
379 |
|
|
! outputs: |
380 |
|
|
INTEGER iflag(len) ! also an input |
381 |
|
|
|
382 |
|
|
! local variables: |
383 |
|
|
INTEGER i |
384 |
|
|
|
385 |
|
|
|
386 |
|
|
DO i = 1, len |
387 |
|
|
IF ((cbmf(i)==0.0) .AND. (iflag(i)==0) .AND. (tvp(i, & |
388 |
|
|
icb(i))<=(tv(i,icb(i))-dtmax))) iflag(i) = 4 |
389 |
|
|
END DO |
390 |
|
|
|
391 |
|
|
RETURN |
392 |
|
|
END SUBROUTINE cv_trigger |
393 |
|
|
|
394 |
|
|
SUBROUTINE cv_compress(len, nloc, ncum, nd, iflag1, nk1, icb1, cbmf1, plcl1, & |
395 |
|
|
tnk1, qnk1, gznk1, t1, q1, qs1, u1, v1, gz1, h1, lv1, cpn1, p1, ph1, tv1, & |
396 |
|
|
tp1, tvp1, clw1, iflag, nk, icb, cbmf, plcl, tnk, qnk, gznk, t, q, qs, u, & |
397 |
|
|
v, gz, h, lv, cpn, p, ph, tv, tp, tvp, clw, dph) |
398 |
|
|
USE print_control_mod, ONLY: lunout |
399 |
|
|
IMPLICIT NONE |
400 |
|
|
|
401 |
|
|
include "cvparam.h" |
402 |
|
|
|
403 |
|
|
! inputs: |
404 |
|
|
INTEGER len, ncum, nd, nloc |
405 |
|
|
INTEGER iflag1(len), nk1(len), icb1(len) |
406 |
|
|
REAL cbmf1(len), plcl1(len), tnk1(len), qnk1(len), gznk1(len) |
407 |
|
|
REAL t1(len, nd), q1(len, nd), qs1(len, nd), u1(len, nd), v1(len, nd) |
408 |
|
|
REAL gz1(len, nd), h1(len, nd), lv1(len, nd), cpn1(len, nd) |
409 |
|
|
REAL p1(len, nd), ph1(len, nd+1), tv1(len, nd), tp1(len, nd) |
410 |
|
|
REAL tvp1(len, nd), clw1(len, nd) |
411 |
|
|
|
412 |
|
|
! outputs: |
413 |
|
|
INTEGER iflag(nloc), nk(nloc), icb(nloc) |
414 |
|
|
REAL cbmf(nloc), plcl(nloc), tnk(nloc), qnk(nloc), gznk(nloc) |
415 |
|
|
REAL t(nloc, nd), q(nloc, nd), qs(nloc, nd), u(nloc, nd), v(nloc, nd) |
416 |
|
|
REAL gz(nloc, nd), h(nloc, nd), lv(nloc, nd), cpn(nloc, nd) |
417 |
|
|
REAL p(nloc, nd), ph(nloc, nd+1), tv(nloc, nd), tp(nloc, nd) |
418 |
|
|
REAL tvp(nloc, nd), clw(nloc, nd) |
419 |
|
|
REAL dph(nloc, nd) |
420 |
|
|
|
421 |
|
|
! local variables: |
422 |
|
|
INTEGER i, k, nn |
423 |
|
|
CHARACTER (LEN=20) :: modname = 'cv_compress' |
424 |
|
|
CHARACTER (LEN=80) :: abort_message |
425 |
|
|
|
426 |
|
|
|
427 |
|
|
DO k = 1, nl + 1 |
428 |
|
|
nn = 0 |
429 |
|
|
DO i = 1, len |
430 |
|
|
IF (iflag1(i)==0) THEN |
431 |
|
|
nn = nn + 1 |
432 |
|
|
t(nn, k) = t1(i, k) |
433 |
|
|
q(nn, k) = q1(i, k) |
434 |
|
|
qs(nn, k) = qs1(i, k) |
435 |
|
|
u(nn, k) = u1(i, k) |
436 |
|
|
v(nn, k) = v1(i, k) |
437 |
|
|
gz(nn, k) = gz1(i, k) |
438 |
|
|
h(nn, k) = h1(i, k) |
439 |
|
|
lv(nn, k) = lv1(i, k) |
440 |
|
|
cpn(nn, k) = cpn1(i, k) |
441 |
|
|
p(nn, k) = p1(i, k) |
442 |
|
|
ph(nn, k) = ph1(i, k) |
443 |
|
|
tv(nn, k) = tv1(i, k) |
444 |
|
|
tp(nn, k) = tp1(i, k) |
445 |
|
|
tvp(nn, k) = tvp1(i, k) |
446 |
|
|
clw(nn, k) = clw1(i, k) |
447 |
|
|
END IF |
448 |
|
|
END DO |
449 |
|
|
END DO |
450 |
|
|
|
451 |
|
|
IF (nn/=ncum) THEN |
452 |
|
|
WRITE (lunout, *) 'strange! nn not equal to ncum: ', nn, ncum |
453 |
|
|
abort_message = '' |
454 |
|
|
CALL abort_physic(modname, abort_message, 1) |
455 |
|
|
END IF |
456 |
|
|
|
457 |
|
|
nn = 0 |
458 |
|
|
DO i = 1, len |
459 |
|
|
IF (iflag1(i)==0) THEN |
460 |
|
|
nn = nn + 1 |
461 |
|
|
cbmf(nn) = cbmf1(i) |
462 |
|
|
plcl(nn) = plcl1(i) |
463 |
|
|
tnk(nn) = tnk1(i) |
464 |
|
|
qnk(nn) = qnk1(i) |
465 |
|
|
gznk(nn) = gznk1(i) |
466 |
|
|
nk(nn) = nk1(i) |
467 |
|
|
icb(nn) = icb1(i) |
468 |
|
|
iflag(nn) = iflag1(i) |
469 |
|
|
END IF |
470 |
|
|
END DO |
471 |
|
|
|
472 |
|
|
DO k = 1, nl |
473 |
|
|
DO i = 1, ncum |
474 |
|
|
dph(i, k) = ph(i, k) - ph(i, k+1) |
475 |
|
|
END DO |
476 |
|
|
END DO |
477 |
|
|
|
478 |
|
|
RETURN |
479 |
|
|
END SUBROUTINE cv_compress |
480 |
|
|
|
481 |
|
|
SUBROUTINE cv_undilute2(nloc, ncum, nd, icb, nk, tnk, qnk, gznk, t, q, qs, & |
482 |
|
|
gz, p, dph, h, tv, lv, inb, inb1, tp, tvp, clw, hp, ep, sigp, frac) |
483 |
|
|
IMPLICIT NONE |
484 |
|
|
|
485 |
|
|
! --------------------------------------------------------------------- |
486 |
|
|
! Purpose: |
487 |
|
|
! FIND THE REST OF THE LIFTED PARCEL TEMPERATURES |
488 |
|
|
! & |
489 |
|
|
! COMPUTE THE PRECIPITATION EFFICIENCIES AND THE |
490 |
|
|
! FRACTION OF PRECIPITATION FALLING OUTSIDE OF CLOUD |
491 |
|
|
! & |
492 |
|
|
! FIND THE LEVEL OF NEUTRAL BUOYANCY |
493 |
|
|
! --------------------------------------------------------------------- |
494 |
|
|
|
495 |
|
|
include "cvthermo.h" |
496 |
|
|
include "cvparam.h" |
497 |
|
|
|
498 |
|
|
! inputs: |
499 |
|
|
INTEGER ncum, nd, nloc |
500 |
|
|
INTEGER icb(nloc), nk(nloc) |
501 |
|
|
REAL t(nloc, nd), q(nloc, nd), qs(nloc, nd), gz(nloc, nd) |
502 |
|
|
REAL p(nloc, nd), dph(nloc, nd) |
503 |
|
|
REAL tnk(nloc), qnk(nloc), gznk(nloc) |
504 |
|
|
REAL lv(nloc, nd), tv(nloc, nd), h(nloc, nd) |
505 |
|
|
|
506 |
|
|
! outputs: |
507 |
|
|
INTEGER inb(nloc), inb1(nloc) |
508 |
|
|
REAL tp(nloc, nd), tvp(nloc, nd), clw(nloc, nd) |
509 |
|
|
REAL ep(nloc, nd), sigp(nloc, nd), hp(nloc, nd) |
510 |
|
|
REAL frac(nloc) |
511 |
|
|
|
512 |
|
|
! local variables: |
513 |
|
|
INTEGER i, k |
514 |
|
|
REAL tg, qg, ahg, alv, s, tc, es, denom, rg, tca, elacrit |
515 |
|
|
REAL by, defrac |
516 |
|
|
REAL ah0(nloc), cape(nloc), capem(nloc), byp(nloc) |
517 |
|
|
LOGICAL lcape(nloc) |
518 |
|
|
|
519 |
|
|
! ===================================================================== |
520 |
|
|
! --- SOME INITIALIZATIONS |
521 |
|
|
! ===================================================================== |
522 |
|
|
|
523 |
|
|
DO k = 1, nl |
524 |
|
|
DO i = 1, ncum |
525 |
|
|
ep(i, k) = 0.0 |
526 |
|
|
sigp(i, k) = sigs |
527 |
|
|
END DO |
528 |
|
|
END DO |
529 |
|
|
|
530 |
|
|
! ===================================================================== |
531 |
|
|
! --- FIND THE REST OF THE LIFTED PARCEL TEMPERATURES |
532 |
|
|
! ===================================================================== |
533 |
|
|
|
534 |
|
|
! --- The procedure is to solve the equation. |
535 |
|
|
! cp*tp+L*qp+phi=cp*tnk+L*qnk+gznk. |
536 |
|
|
|
537 |
|
|
! *** Calculate certain parcel quantities, including static energy *** |
538 |
|
|
|
539 |
|
|
|
540 |
|
|
DO i = 1, ncum |
541 |
|
|
ah0(i) = (cpd*(1.-qnk(i))+cl*qnk(i))*tnk(i) + qnk(i)*(lv0-clmcpv*(tnk(i)- & |
542 |
|
|
t0)) + gznk(i) |
543 |
|
|
END DO |
544 |
|
|
|
545 |
|
|
|
546 |
|
|
! *** Find lifted parcel quantities above cloud base *** |
547 |
|
|
|
548 |
|
|
|
549 |
|
|
DO k = minorig + 1, nl |
550 |
|
|
DO i = 1, ncum |
551 |
|
|
IF (k>=(icb(i)+1)) THEN |
552 |
|
|
tg = t(i, k) |
553 |
|
|
qg = qs(i, k) |
554 |
|
|
alv = lv0 - clmcpv*(t(i,k)-t0) |
555 |
|
|
|
556 |
|
|
! First iteration. |
557 |
|
|
|
558 |
|
|
s = cpd + alv*alv*qg/(rrv*t(i,k)*t(i,k)) |
559 |
|
|
s = 1./s |
560 |
|
|
ahg = cpd*tg + (cl-cpd)*qnk(i)*t(i, k) + alv*qg + gz(i, k) |
561 |
|
|
tg = tg + s*(ah0(i)-ahg) |
562 |
|
|
tg = max(tg, 35.0) |
563 |
|
|
tc = tg - t0 |
564 |
|
|
denom = 243.5 + tc |
565 |
|
|
IF (tc>=0.0) THEN |
566 |
|
|
es = 6.112*exp(17.67*tc/denom) |
567 |
|
|
ELSE |
568 |
|
|
es = exp(23.33086-6111.72784/tg+0.15215*log(tg)) |
569 |
|
|
END IF |
570 |
|
|
qg = eps*es/(p(i,k)-es*(1.-eps)) |
571 |
|
|
|
572 |
|
|
! Second iteration. |
573 |
|
|
|
574 |
|
|
s = cpd + alv*alv*qg/(rrv*t(i,k)*t(i,k)) |
575 |
|
|
s = 1./s |
576 |
|
|
ahg = cpd*tg + (cl-cpd)*qnk(i)*t(i, k) + alv*qg + gz(i, k) |
577 |
|
|
tg = tg + s*(ah0(i)-ahg) |
578 |
|
|
tg = max(tg, 35.0) |
579 |
|
|
tc = tg - t0 |
580 |
|
|
denom = 243.5 + tc |
581 |
|
|
IF (tc>=0.0) THEN |
582 |
|
|
es = 6.112*exp(17.67*tc/denom) |
583 |
|
|
ELSE |
584 |
|
|
es = exp(23.33086-6111.72784/tg+0.15215*log(tg)) |
585 |
|
|
END IF |
586 |
|
|
qg = eps*es/(p(i,k)-es*(1.-eps)) |
587 |
|
|
|
588 |
|
|
alv = lv0 - clmcpv*(t(i,k)-t0) |
589 |
|
|
! print*,'cpd dans convect2 ',cpd |
590 |
|
|
! print*,'tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd' |
591 |
|
|
! print*,tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd |
592 |
|
|
tp(i, k) = (ah0(i)-(cl-cpd)*qnk(i)*t(i,k)-gz(i,k)-alv*qg)/cpd |
593 |
|
|
! if (.not.cpd.gt.1000.) then |
594 |
|
|
! print*,'CPD=',cpd |
595 |
|
|
! stop |
596 |
|
|
! endif |
597 |
|
|
clw(i, k) = qnk(i) - qg |
598 |
|
|
clw(i, k) = max(0.0, clw(i,k)) |
599 |
|
|
rg = qg/(1.-qnk(i)) |
600 |
|
|
tvp(i, k) = tp(i, k)*(1.+rg*epsi) |
601 |
|
|
END IF |
602 |
|
|
END DO |
603 |
|
|
END DO |
604 |
|
|
|
605 |
|
|
! ===================================================================== |
606 |
|
|
! --- SET THE PRECIPITATION EFFICIENCIES AND THE FRACTION OF |
607 |
|
|
! --- PRECIPITATION FALLING OUTSIDE OF CLOUD |
608 |
|
|
! --- THESE MAY BE FUNCTIONS OF TP(I), P(I) AND CLW(I) |
609 |
|
|
! ===================================================================== |
610 |
|
|
|
611 |
|
|
DO k = minorig + 1, nl |
612 |
|
|
DO i = 1, ncum |
613 |
|
|
IF (k>=(nk(i)+1)) THEN |
614 |
|
|
tca = tp(i, k) - t0 |
615 |
|
|
IF (tca>=0.0) THEN |
616 |
|
|
elacrit = elcrit |
617 |
|
|
ELSE |
618 |
|
|
elacrit = elcrit*(1.0-tca/tlcrit) |
619 |
|
|
END IF |
620 |
|
|
elacrit = max(elacrit, 0.0) |
621 |
|
|
ep(i, k) = 1.0 - elacrit/max(clw(i,k), 1.0E-8) |
622 |
|
|
ep(i, k) = max(ep(i,k), 0.0) |
623 |
|
|
ep(i, k) = min(ep(i,k), 1.0) |
624 |
|
|
sigp(i, k) = sigs |
625 |
|
|
END IF |
626 |
|
|
END DO |
627 |
|
|
END DO |
628 |
|
|
|
629 |
|
|
! ===================================================================== |
630 |
|
|
! --- CALCULATE VIRTUAL TEMPERATURE AND LIFTED PARCEL |
631 |
|
|
! --- VIRTUAL TEMPERATURE |
632 |
|
|
! ===================================================================== |
633 |
|
|
|
634 |
|
|
DO k = minorig + 1, nl |
635 |
|
|
DO i = 1, ncum |
636 |
|
|
IF (k>=(icb(i)+1)) THEN |
637 |
|
|
tvp(i, k) = tvp(i, k)*(1.0-qnk(i)+ep(i,k)*clw(i,k)) |
638 |
|
|
! print*,'i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)' |
639 |
|
|
! print*, i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k) |
640 |
|
|
END IF |
641 |
|
|
END DO |
642 |
|
|
END DO |
643 |
|
|
DO i = 1, ncum |
644 |
|
|
tvp(i, nlp) = tvp(i, nl) - (gz(i,nlp)-gz(i,nl))/cpd |
645 |
|
|
END DO |
646 |
|
|
|
647 |
|
|
! ===================================================================== |
648 |
|
|
! --- FIND THE FIRST MODEL LEVEL (INB1) ABOVE THE PARCEL'S |
649 |
|
|
! --- HIGHEST LEVEL OF NEUTRAL BUOYANCY |
650 |
|
|
! --- AND THE HIGHEST LEVEL OF POSITIVE CAPE (INB) |
651 |
|
|
! ===================================================================== |
652 |
|
|
|
653 |
|
|
DO i = 1, ncum |
654 |
|
|
cape(i) = 0.0 |
655 |
|
|
capem(i) = 0.0 |
656 |
|
|
inb(i) = icb(i) + 1 |
657 |
|
|
inb1(i) = inb(i) |
658 |
|
|
END DO |
659 |
|
|
|
660 |
|
|
! Originial Code |
661 |
|
|
|
662 |
|
|
! do 530 k=minorig+1,nl-1 |
663 |
|
|
! do 520 i=1,ncum |
664 |
|
|
! if(k.ge.(icb(i)+1))then |
665 |
|
|
! by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k) |
666 |
|
|
! byp=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1) |
667 |
|
|
! cape(i)=cape(i)+by |
668 |
|
|
! if(by.ge.0.0)inb1(i)=k+1 |
669 |
|
|
! if(cape(i).gt.0.0)then |
670 |
|
|
! inb(i)=k+1 |
671 |
|
|
! capem(i)=cape(i) |
672 |
|
|
! endif |
673 |
|
|
! endif |
674 |
|
|
! 520 continue |
675 |
|
|
! 530 continue |
676 |
|
|
! do 540 i=1,ncum |
677 |
|
|
! byp=(tvp(i,nl)-tv(i,nl))*dph(i,nl)/p(i,nl) |
678 |
|
|
! cape(i)=capem(i)+byp |
679 |
|
|
! defrac=capem(i)-cape(i) |
680 |
|
|
! defrac=max(defrac,0.001) |
681 |
|
|
! frac(i)=-cape(i)/defrac |
682 |
|
|
! frac(i)=min(frac(i),1.0) |
683 |
|
|
! frac(i)=max(frac(i),0.0) |
684 |
|
|
! 540 continue |
685 |
|
|
|
686 |
|
|
! K Emanuel fix |
687 |
|
|
|
688 |
|
|
! call zilch(byp,ncum) |
689 |
|
|
! do 530 k=minorig+1,nl-1 |
690 |
|
|
! do 520 i=1,ncum |
691 |
|
|
! if(k.ge.(icb(i)+1))then |
692 |
|
|
! by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k) |
693 |
|
|
! cape(i)=cape(i)+by |
694 |
|
|
! if(by.ge.0.0)inb1(i)=k+1 |
695 |
|
|
! if(cape(i).gt.0.0)then |
696 |
|
|
! inb(i)=k+1 |
697 |
|
|
! capem(i)=cape(i) |
698 |
|
|
! byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1) |
699 |
|
|
! endif |
700 |
|
|
! endif |
701 |
|
|
! 520 continue |
702 |
|
|
! 530 continue |
703 |
|
|
! do 540 i=1,ncum |
704 |
|
|
! inb(i)=max(inb(i),inb1(i)) |
705 |
|
|
! cape(i)=capem(i)+byp(i) |
706 |
|
|
! defrac=capem(i)-cape(i) |
707 |
|
|
! defrac=max(defrac,0.001) |
708 |
|
|
! frac(i)=-cape(i)/defrac |
709 |
|
|
! frac(i)=min(frac(i),1.0) |
710 |
|
|
! frac(i)=max(frac(i),0.0) |
711 |
|
|
! 540 continue |
712 |
|
|
|
713 |
|
|
! J Teixeira fix |
714 |
|
|
|
715 |
|
|
CALL zilch(byp, ncum) |
716 |
|
|
DO i = 1, ncum |
717 |
|
|
lcape(i) = .TRUE. |
718 |
|
|
END DO |
719 |
|
|
DO k = minorig + 1, nl - 1 |
720 |
|
|
DO i = 1, ncum |
721 |
|
|
IF (cape(i)<0.0) lcape(i) = .FALSE. |
722 |
|
|
IF ((k>=(icb(i)+1)) .AND. lcape(i)) THEN |
723 |
|
|
by = (tvp(i,k)-tv(i,k))*dph(i, k)/p(i, k) |
724 |
|
|
byp(i) = (tvp(i,k+1)-tv(i,k+1))*dph(i, k+1)/p(i, k+1) |
725 |
|
|
cape(i) = cape(i) + by |
726 |
|
|
IF (by>=0.0) inb1(i) = k + 1 |
727 |
|
|
IF (cape(i)>0.0) THEN |
728 |
|
|
inb(i) = k + 1 |
729 |
|
|
capem(i) = cape(i) |
730 |
|
|
END IF |
731 |
|
|
END IF |
732 |
|
|
END DO |
733 |
|
|
END DO |
734 |
|
|
DO i = 1, ncum |
735 |
|
|
cape(i) = capem(i) + byp(i) |
736 |
|
|
defrac = capem(i) - cape(i) |
737 |
|
|
defrac = max(defrac, 0.001) |
738 |
|
|
frac(i) = -cape(i)/defrac |
739 |
|
|
frac(i) = min(frac(i), 1.0) |
740 |
|
|
frac(i) = max(frac(i), 0.0) |
741 |
|
|
END DO |
742 |
|
|
|
743 |
|
|
! ===================================================================== |
744 |
|
|
! --- CALCULATE LIQUID WATER STATIC ENERGY OF LIFTED PARCEL |
745 |
|
|
! ===================================================================== |
746 |
|
|
|
747 |
|
|
! initialization: |
748 |
|
|
DO i = 1, ncum*nlp |
749 |
|
|
hp(i, 1) = h(i, 1) |
750 |
|
|
END DO |
751 |
|
|
|
752 |
|
|
DO k = minorig + 1, nl |
753 |
|
|
DO i = 1, ncum |
754 |
|
|
IF ((k>=icb(i)) .AND. (k<=inb(i))) THEN |
755 |
|
|
hp(i, k) = h(i, nk(i)) + (lv(i,k)+(cpd-cpv)*t(i,k))*ep(i, k)*clw(i, k & |
756 |
|
|
) |
757 |
|
|
END IF |
758 |
|
|
END DO |
759 |
|
|
END DO |
760 |
|
|
|
761 |
|
|
RETURN |
762 |
|
|
END SUBROUTINE cv_undilute2 |
763 |
|
|
|
764 |
|
|
SUBROUTINE cv_closure(nloc, ncum, nd, nk, icb, tv, tvp, p, ph, dph, plcl, & |
765 |
|
|
cpn, iflag, cbmf) |
766 |
|
|
IMPLICIT NONE |
767 |
|
|
|
768 |
|
|
! inputs: |
769 |
|
|
INTEGER ncum, nd, nloc |
770 |
|
|
INTEGER nk(nloc), icb(nloc) |
771 |
|
|
REAL tv(nloc, nd), tvp(nloc, nd), p(nloc, nd), dph(nloc, nd) |
772 |
|
|
REAL ph(nloc, nd+1) ! caution nd instead ndp1 to be consistent... |
773 |
|
|
REAL plcl(nloc), cpn(nloc, nd) |
774 |
|
|
|
775 |
|
|
! outputs: |
776 |
|
|
INTEGER iflag(nloc) |
777 |
|
|
REAL cbmf(nloc) ! also an input |
778 |
|
|
|
779 |
|
|
! local variables: |
780 |
|
|
INTEGER i, k, icbmax |
781 |
|
|
REAL dtpbl(nloc), dtmin(nloc), tvpplcl(nloc), tvaplcl(nloc) |
782 |
|
|
REAL work(nloc) |
783 |
|
|
|
784 |
|
|
include "cvthermo.h" |
785 |
|
|
include "cvparam.h" |
786 |
|
|
|
787 |
|
|
! ------------------------------------------------------------------- |
788 |
|
|
! Compute icbmax. |
789 |
|
|
! ------------------------------------------------------------------- |
790 |
|
|
|
791 |
|
|
icbmax = 2 |
792 |
|
|
DO i = 1, ncum |
793 |
|
|
icbmax = max(icbmax, icb(i)) |
794 |
|
|
END DO |
795 |
|
|
|
796 |
|
|
! ===================================================================== |
797 |
|
|
! --- CALCULATE CLOUD BASE MASS FLUX |
798 |
|
|
! ===================================================================== |
799 |
|
|
|
800 |
|
|
! tvpplcl = parcel temperature lifted adiabatically from level |
801 |
|
|
! icb-1 to the LCL. |
802 |
|
|
! tvaplcl = virtual temperature at the LCL. |
803 |
|
|
|
804 |
|
|
DO i = 1, ncum |
805 |
|
|
dtpbl(i) = 0.0 |
806 |
|
|
tvpplcl(i) = tvp(i, icb(i)-1) - rrd*tvp(i, icb(i)-1)*(p(i,icb(i)-1)-plcl( & |
807 |
|
|
i))/(cpn(i,icb(i)-1)*p(i,icb(i)-1)) |
808 |
|
|
tvaplcl(i) = tv(i, icb(i)) + (tvp(i,icb(i))-tvp(i,icb(i)+1))*(plcl(i)-p(i & |
809 |
|
|
,icb(i)))/(p(i,icb(i))-p(i,icb(i)+1)) |
810 |
|
|
END DO |
811 |
|
|
|
812 |
|
|
! ------------------------------------------------------------------- |
813 |
|
|
! --- Interpolate difference between lifted parcel and |
814 |
|
|
! --- environmental temperatures to lifted condensation level |
815 |
|
|
! ------------------------------------------------------------------- |
816 |
|
|
|
817 |
|
|
! dtpbl = average of tvp-tv in the PBL (k=nk to icb-1). |
818 |
|
|
|
819 |
|
|
DO k = minorig, icbmax |
820 |
|
|
DO i = 1, ncum |
821 |
|
|
IF ((k>=nk(i)) .AND. (k<=(icb(i)-1))) THEN |
822 |
|
|
dtpbl(i) = dtpbl(i) + (tvp(i,k)-tv(i,k))*dph(i, k) |
823 |
|
|
END IF |
824 |
|
|
END DO |
825 |
|
|
END DO |
826 |
|
|
DO i = 1, ncum |
827 |
|
|
dtpbl(i) = dtpbl(i)/(ph(i,nk(i))-ph(i,icb(i))) |
828 |
|
|
dtmin(i) = tvpplcl(i) - tvaplcl(i) + dtmax + dtpbl(i) |
829 |
|
|
END DO |
830 |
|
|
|
831 |
|
|
! ------------------------------------------------------------------- |
832 |
|
|
! --- Adjust cloud base mass flux |
833 |
|
|
! ------------------------------------------------------------------- |
834 |
|
|
|
835 |
|
|
DO i = 1, ncum |
836 |
|
|
work(i) = cbmf(i) |
837 |
|
|
cbmf(i) = max(0.0, (1.0-damp)*cbmf(i)+0.1*alpha*dtmin(i)) |
838 |
|
|
IF ((work(i)==0.0) .AND. (cbmf(i)==0.0)) THEN |
839 |
|
|
iflag(i) = 3 |
840 |
|
|
END IF |
841 |
|
|
END DO |
842 |
|
|
|
843 |
|
|
RETURN |
844 |
|
|
END SUBROUTINE cv_closure |
845 |
|
|
|
846 |
|
|
SUBROUTINE cv_mixing(nloc, ncum, nd, icb, nk, inb, inb1, ph, t, q, qs, u, v, & |
847 |
|
|
h, lv, qnk, hp, tv, tvp, ep, clw, cbmf, m, ment, qent, uent, vent, nent, & |
848 |
|
|
sij, elij) |
849 |
|
|
IMPLICIT NONE |
850 |
|
|
|
851 |
|
|
include "cvthermo.h" |
852 |
|
|
include "cvparam.h" |
853 |
|
|
|
854 |
|
|
! inputs: |
855 |
|
|
INTEGER ncum, nd, nloc |
856 |
|
|
INTEGER icb(nloc), inb(nloc), inb1(nloc), nk(nloc) |
857 |
|
|
REAL cbmf(nloc), qnk(nloc) |
858 |
|
|
REAL ph(nloc, nd+1) |
859 |
|
|
REAL t(nloc, nd), q(nloc, nd), qs(nloc, nd), lv(nloc, nd) |
860 |
|
|
REAL u(nloc, nd), v(nloc, nd), h(nloc, nd), hp(nloc, nd) |
861 |
|
|
REAL tv(nloc, nd), tvp(nloc, nd), ep(nloc, nd), clw(nloc, nd) |
862 |
|
|
|
863 |
|
|
! outputs: |
864 |
|
|
INTEGER nent(nloc, nd) |
865 |
|
|
REAL m(nloc, nd), ment(nloc, nd, nd), qent(nloc, nd, nd) |
866 |
|
|
REAL uent(nloc, nd, nd), vent(nloc, nd, nd) |
867 |
|
|
REAL sij(nloc, nd, nd), elij(nloc, nd, nd) |
868 |
|
|
|
869 |
|
|
! local variables: |
870 |
|
|
INTEGER i, j, k, ij |
871 |
|
|
INTEGER num1, num2 |
872 |
|
|
REAL dbo, qti, bf2, anum, denom, dei, altem, cwat, stemp |
873 |
|
|
REAL alt, qp1, smid, sjmin, sjmax, delp, delm |
874 |
|
|
REAL work(nloc), asij(nloc), smin(nloc), scrit(nloc) |
875 |
|
|
REAL bsum(nloc, nd) |
876 |
|
|
LOGICAL lwork(nloc) |
877 |
|
|
|
878 |
|
|
! ===================================================================== |
879 |
|
|
! --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS |
880 |
|
|
! ===================================================================== |
881 |
|
|
|
882 |
|
|
DO i = 1, ncum*nlp |
883 |
|
|
nent(i, 1) = 0 |
884 |
|
|
m(i, 1) = 0.0 |
885 |
|
|
END DO |
886 |
|
|
|
887 |
|
|
DO k = 1, nlp |
888 |
|
|
DO j = 1, nlp |
889 |
|
|
DO i = 1, ncum |
890 |
|
|
qent(i, k, j) = q(i, j) |
891 |
|
|
uent(i, k, j) = u(i, j) |
892 |
|
|
vent(i, k, j) = v(i, j) |
893 |
|
|
elij(i, k, j) = 0.0 |
894 |
|
|
ment(i, k, j) = 0.0 |
895 |
|
|
sij(i, k, j) = 0.0 |
896 |
|
|
END DO |
897 |
|
|
END DO |
898 |
|
|
END DO |
899 |
|
|
|
900 |
|
|
! ------------------------------------------------------------------- |
901 |
|
|
! --- Calculate rates of mixing, m(i) |
902 |
|
|
! ------------------------------------------------------------------- |
903 |
|
|
|
904 |
|
|
CALL zilch(work, ncum) |
905 |
|
|
|
906 |
|
|
DO j = minorig + 1, nl |
907 |
|
|
DO i = 1, ncum |
908 |
|
|
IF ((j>=(icb(i)+1)) .AND. (j<=inb(i))) THEN |
909 |
|
|
k = min(j, inb1(i)) |
910 |
|
|
dbo = abs(tv(i,k+1)-tvp(i,k+1)-tv(i,k-1)+tvp(i,k-1)) + & |
911 |
|
|
entp*0.04*(ph(i,k)-ph(i,k+1)) |
912 |
|
|
work(i) = work(i) + dbo |
913 |
|
|
m(i, j) = cbmf(i)*dbo |
914 |
|
|
END IF |
915 |
|
|
END DO |
916 |
|
|
END DO |
917 |
|
|
DO k = minorig + 1, nl |
918 |
|
|
DO i = 1, ncum |
919 |
|
|
IF ((k>=(icb(i)+1)) .AND. (k<=inb(i))) THEN |
920 |
|
|
m(i, k) = m(i, k)/work(i) |
921 |
|
|
END IF |
922 |
|
|
END DO |
923 |
|
|
END DO |
924 |
|
|
|
925 |
|
|
|
926 |
|
|
! ===================================================================== |
927 |
|
|
! --- CALCULATE ENTRAINED AIR MASS FLUX (ment), TOTAL WATER MIXING |
928 |
|
|
! --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING |
929 |
|
|
! --- FRACTION (sij) |
930 |
|
|
! ===================================================================== |
931 |
|
|
|
932 |
|
|
|
933 |
|
|
DO i = minorig + 1, nl |
934 |
|
|
DO j = minorig + 1, nl |
935 |
|
|
DO ij = 1, ncum |
936 |
|
|
IF ((i>=(icb(ij)+1)) .AND. (j>=icb(ij)) .AND. (i<=inb(ij)) .AND. (j<= & |
937 |
|
|
inb(ij))) THEN |
938 |
|
|
qti = qnk(ij) - ep(ij, i)*clw(ij, i) |
939 |
|
|
bf2 = 1. + lv(ij, j)*lv(ij, j)*qs(ij, j)/(rrv*t(ij,j)*t(ij,j)*cpd) |
940 |
|
|
anum = h(ij, j) - hp(ij, i) + (cpv-cpd)*t(ij, j)*(qti-q(ij,j)) |
941 |
|
|
denom = h(ij, i) - hp(ij, i) + (cpd-cpv)*(q(ij,i)-qti)*t(ij, j) |
942 |
|
|
dei = denom |
943 |
|
|
IF (abs(dei)<0.01) dei = 0.01 |
944 |
|
|
sij(ij, i, j) = anum/dei |
945 |
|
|
sij(ij, i, i) = 1.0 |
946 |
|
|
altem = sij(ij, i, j)*q(ij, i) + (1.-sij(ij,i,j))*qti - qs(ij, j) |
947 |
|
|
altem = altem/bf2 |
948 |
|
|
cwat = clw(ij, j)*(1.-ep(ij,j)) |
949 |
|
|
stemp = sij(ij, i, j) |
950 |
|
|
IF ((stemp<0.0 .OR. stemp>1.0 .OR. altem>cwat) .AND. j>i) THEN |
951 |
|
|
anum = anum - lv(ij, j)*(qti-qs(ij,j)-cwat*bf2) |
952 |
|
|
denom = denom + lv(ij, j)*(q(ij,i)-qti) |
953 |
|
|
IF (abs(denom)<0.01) denom = 0.01 |
954 |
|
|
sij(ij, i, j) = anum/denom |
955 |
|
|
altem = sij(ij, i, j)*q(ij, i) + (1.-sij(ij,i,j))*qti - qs(ij, j) |
956 |
|
|
altem = altem - (bf2-1.)*cwat |
957 |
|
|
END IF |
958 |
|
|
IF (sij(ij,i,j)>0.0 .AND. sij(ij,i,j)<0.9) THEN |
959 |
|
|
qent(ij, i, j) = sij(ij, i, j)*q(ij, i) + (1.-sij(ij,i,j))*qti |
960 |
|
|
uent(ij, i, j) = sij(ij, i, j)*u(ij, i) + & |
961 |
|
|
(1.-sij(ij,i,j))*u(ij, nk(ij)) |
962 |
|
|
vent(ij, i, j) = sij(ij, i, j)*v(ij, i) + & |
963 |
|
|
(1.-sij(ij,i,j))*v(ij, nk(ij)) |
964 |
|
|
elij(ij, i, j) = altem |
965 |
|
|
elij(ij, i, j) = max(0.0, elij(ij,i,j)) |
966 |
|
|
ment(ij, i, j) = m(ij, i)/(1.-sij(ij,i,j)) |
967 |
|
|
nent(ij, i) = nent(ij, i) + 1 |
968 |
|
|
END IF |
969 |
|
|
sij(ij, i, j) = max(0.0, sij(ij,i,j)) |
970 |
|
|
sij(ij, i, j) = min(1.0, sij(ij,i,j)) |
971 |
|
|
END IF |
972 |
|
|
END DO |
973 |
|
|
END DO |
974 |
|
|
|
975 |
|
|
! *** If no air can entrain at level i assume that updraft detrains |
976 |
|
|
! *** |
977 |
|
|
! *** at that level and calculate detrained air flux and properties |
978 |
|
|
! *** |
979 |
|
|
|
980 |
|
|
DO ij = 1, ncum |
981 |
|
|
IF ((i>=(icb(ij)+1)) .AND. (i<=inb(ij)) .AND. (nent(ij,i)==0)) THEN |
982 |
|
|
ment(ij, i, i) = m(ij, i) |
983 |
|
|
qent(ij, i, i) = q(ij, nk(ij)) - ep(ij, i)*clw(ij, i) |
984 |
|
|
uent(ij, i, i) = u(ij, nk(ij)) |
985 |
|
|
vent(ij, i, i) = v(ij, nk(ij)) |
986 |
|
|
elij(ij, i, i) = clw(ij, i) |
987 |
|
|
sij(ij, i, i) = 1.0 |
988 |
|
|
END IF |
989 |
|
|
END DO |
990 |
|
|
END DO |
991 |
|
|
|
992 |
|
|
DO i = 1, ncum |
993 |
|
|
sij(i, inb(i), inb(i)) = 1.0 |
994 |
|
|
END DO |
995 |
|
|
|
996 |
|
|
! ===================================================================== |
997 |
|
|
! --- NORMALIZE ENTRAINED AIR MASS FLUXES |
998 |
|
|
! --- TO REPRESENT EQUAL PROBABILITIES OF MIXING |
999 |
|
|
! ===================================================================== |
1000 |
|
|
|
1001 |
|
|
CALL zilch(bsum, ncum*nlp) |
1002 |
|
|
DO ij = 1, ncum |
1003 |
|
|
lwork(ij) = .FALSE. |
1004 |
|
|
END DO |
1005 |
|
|
DO i = minorig + 1, nl |
1006 |
|
|
|
1007 |
|
|
num1 = 0 |
1008 |
|
|
DO ij = 1, ncum |
1009 |
|
|
IF ((i>=icb(ij)+1) .AND. (i<=inb(ij))) num1 = num1 + 1 |
1010 |
|
|
END DO |
1011 |
|
|
IF (num1<=0) GO TO 789 |
1012 |
|
|
|
1013 |
|
|
DO ij = 1, ncum |
1014 |
|
|
IF ((i>=icb(ij)+1) .AND. (i<=inb(ij))) THEN |
1015 |
|
|
lwork(ij) = (nent(ij,i)/=0) |
1016 |
|
|
qp1 = q(ij, nk(ij)) - ep(ij, i)*clw(ij, i) |
1017 |
|
|
anum = h(ij, i) - hp(ij, i) - lv(ij, i)*(qp1-qs(ij,i)) |
1018 |
|
|
denom = h(ij, i) - hp(ij, i) + lv(ij, i)*(q(ij,i)-qp1) |
1019 |
|
|
IF (abs(denom)<0.01) denom = 0.01 |
1020 |
|
|
scrit(ij) = anum/denom |
1021 |
|
|
alt = qp1 - qs(ij, i) + scrit(ij)*(q(ij,i)-qp1) |
1022 |
|
|
IF (scrit(ij)<0.0 .OR. alt<0.0) scrit(ij) = 1.0 |
1023 |
|
|
asij(ij) = 0.0 |
1024 |
|
|
smin(ij) = 1.0 |
1025 |
|
|
END IF |
1026 |
|
|
END DO |
1027 |
|
|
DO j = minorig, nl |
1028 |
|
|
|
1029 |
|
|
num2 = 0 |
1030 |
|
|
DO ij = 1, ncum |
1031 |
|
|
IF ((i>=icb(ij)+1) .AND. (i<=inb(ij)) .AND. (j>=icb( & |
1032 |
|
|
ij)) .AND. (j<=inb(ij)) .AND. lwork(ij)) num2 = num2 + 1 |
1033 |
|
|
END DO |
1034 |
|
|
IF (num2<=0) GO TO 783 |
1035 |
|
|
|
1036 |
|
|
DO ij = 1, ncum |
1037 |
|
|
IF ((i>=icb(ij)+1) .AND. (i<=inb(ij)) .AND. (j>=icb( & |
1038 |
|
|
ij)) .AND. (j<=inb(ij)) .AND. lwork(ij)) THEN |
1039 |
|
|
IF (sij(ij,i,j)>0.0 .AND. sij(ij,i,j)<0.9) THEN |
1040 |
|
|
IF (j>i) THEN |
1041 |
|
|
smid = min(sij(ij,i,j), scrit(ij)) |
1042 |
|
|
sjmax = smid |
1043 |
|
|
sjmin = smid |
1044 |
|
|
IF (smid<smin(ij) .AND. sij(ij,i,j+1)<smid) THEN |
1045 |
|
|
smin(ij) = smid |
1046 |
|
|
sjmax = min(sij(ij,i,j+1), sij(ij,i,j), scrit(ij)) |
1047 |
|
|
sjmin = max(sij(ij,i,j-1), sij(ij,i,j)) |
1048 |
|
|
sjmin = min(sjmin, scrit(ij)) |
1049 |
|
|
END IF |
1050 |
|
|
ELSE |
1051 |
|
|
sjmax = max(sij(ij,i,j+1), scrit(ij)) |
1052 |
|
|
smid = max(sij(ij,i,j), scrit(ij)) |
1053 |
|
|
sjmin = 0.0 |
1054 |
|
|
IF (j>1) sjmin = sij(ij, i, j-1) |
1055 |
|
|
sjmin = max(sjmin, scrit(ij)) |
1056 |
|
|
END IF |
1057 |
|
|
delp = abs(sjmax-smid) |
1058 |
|
|
delm = abs(sjmin-smid) |
1059 |
|
|
asij(ij) = asij(ij) + (delp+delm)*(ph(ij,j)-ph(ij,j+1)) |
1060 |
|
|
ment(ij, i, j) = ment(ij, i, j)*(delp+delm)*(ph(ij,j)-ph(ij,j+1)) |
1061 |
|
|
END IF |
1062 |
|
|
END IF |
1063 |
|
|
END DO |
1064 |
|
|
783 END DO |
1065 |
|
|
DO ij = 1, ncum |
1066 |
|
|
IF ((i>=icb(ij)+1) .AND. (i<=inb(ij)) .AND. lwork(ij)) THEN |
1067 |
|
|
asij(ij) = max(1.0E-21, asij(ij)) |
1068 |
|
|
asij(ij) = 1.0/asij(ij) |
1069 |
|
|
bsum(ij, i) = 0.0 |
1070 |
|
|
END IF |
1071 |
|
|
END DO |
1072 |
|
|
DO j = minorig, nl + 1 |
1073 |
|
|
DO ij = 1, ncum |
1074 |
|
|
IF ((i>=icb(ij)+1) .AND. (i<=inb(ij)) .AND. (j>=icb( & |
1075 |
|
|
ij)) .AND. (j<=inb(ij)) .AND. lwork(ij)) THEN |
1076 |
|
|
ment(ij, i, j) = ment(ij, i, j)*asij(ij) |
1077 |
|
|
bsum(ij, i) = bsum(ij, i) + ment(ij, i, j) |
1078 |
|
|
END IF |
1079 |
|
|
END DO |
1080 |
|
|
END DO |
1081 |
|
|
DO ij = 1, ncum |
1082 |
|
|
IF ((i>=icb(ij)+1) .AND. (i<=inb(ij)) .AND. (bsum(ij, & |
1083 |
|
|
i)<1.0E-18) .AND. lwork(ij)) THEN |
1084 |
|
|
nent(ij, i) = 0 |
1085 |
|
|
ment(ij, i, i) = m(ij, i) |
1086 |
|
|
qent(ij, i, i) = q(ij, nk(ij)) - ep(ij, i)*clw(ij, i) |
1087 |
|
|
uent(ij, i, i) = u(ij, nk(ij)) |
1088 |
|
|
vent(ij, i, i) = v(ij, nk(ij)) |
1089 |
|
|
elij(ij, i, i) = clw(ij, i) |
1090 |
|
|
sij(ij, i, i) = 1.0 |
1091 |
|
|
END IF |
1092 |
|
|
END DO |
1093 |
|
|
789 END DO |
1094 |
|
|
|
1095 |
|
|
RETURN |
1096 |
|
|
END SUBROUTINE cv_mixing |
1097 |
|
|
|
1098 |
|
|
SUBROUTINE cv_unsat(nloc, ncum, nd, inb, t, q, qs, gz, u, v, p, ph, h, lv, & |
1099 |
|
|
ep, sigp, clw, m, ment, elij, iflag, mp, qp, up, vp, wt, water, evap) |
1100 |
|
|
IMPLICIT NONE |
1101 |
|
|
|
1102 |
|
|
|
1103 |
|
|
include "cvthermo.h" |
1104 |
|
|
include "cvparam.h" |
1105 |
|
|
|
1106 |
|
|
! inputs: |
1107 |
|
|
INTEGER ncum, nd, nloc |
1108 |
|
|
INTEGER inb(nloc) |
1109 |
|
|
REAL t(nloc, nd), q(nloc, nd), qs(nloc, nd) |
1110 |
|
|
REAL gz(nloc, nd), u(nloc, nd), v(nloc, nd) |
1111 |
|
|
REAL p(nloc, nd), ph(nloc, nd+1), h(nloc, nd) |
1112 |
|
|
REAL lv(nloc, nd), ep(nloc, nd), sigp(nloc, nd), clw(nloc, nd) |
1113 |
|
|
REAL m(nloc, nd), ment(nloc, nd, nd), elij(nloc, nd, nd) |
1114 |
|
|
|
1115 |
|
|
! outputs: |
1116 |
|
|
INTEGER iflag(nloc) ! also an input |
1117 |
|
|
REAL mp(nloc, nd), qp(nloc, nd), up(nloc, nd), vp(nloc, nd) |
1118 |
|
|
REAL water(nloc, nd), evap(nloc, nd), wt(nloc, nd) |
1119 |
|
|
|
1120 |
|
|
! local variables: |
1121 |
|
|
INTEGER i, j, k, ij, num1 |
1122 |
|
|
INTEGER jtt(nloc) |
1123 |
|
|
REAL awat, coeff, qsm, afac, sigt, b6, c6, revap |
1124 |
|
|
REAL dhdp, fac, qstm, rat |
1125 |
|
|
REAL wdtrain(nloc) |
1126 |
|
|
LOGICAL lwork(nloc) |
1127 |
|
|
|
1128 |
|
|
! ===================================================================== |
1129 |
|
|
! --- PRECIPITATING DOWNDRAFT CALCULATION |
1130 |
|
|
! ===================================================================== |
1131 |
|
|
|
1132 |
|
|
! Initializations: |
1133 |
|
|
|
1134 |
|
|
DO i = 1, ncum |
1135 |
|
|
DO k = 1, nl + 1 |
1136 |
|
|
wt(i, k) = omtsnow |
1137 |
|
|
mp(i, k) = 0.0 |
1138 |
|
|
evap(i, k) = 0.0 |
1139 |
|
|
water(i, k) = 0.0 |
1140 |
|
|
END DO |
1141 |
|
|
END DO |
1142 |
|
|
|
1143 |
|
|
DO i = 1, ncum |
1144 |
|
|
qp(i, 1) = q(i, 1) |
1145 |
|
|
up(i, 1) = u(i, 1) |
1146 |
|
|
vp(i, 1) = v(i, 1) |
1147 |
|
|
END DO |
1148 |
|
|
|
1149 |
|
|
DO k = 2, nl + 1 |
1150 |
|
|
DO i = 1, ncum |
1151 |
|
|
qp(i, k) = q(i, k-1) |
1152 |
|
|
up(i, k) = u(i, k-1) |
1153 |
|
|
vp(i, k) = v(i, k-1) |
1154 |
|
|
END DO |
1155 |
|
|
END DO |
1156 |
|
|
|
1157 |
|
|
|
1158 |
|
|
! *** Check whether ep(inb)=0, if so, skip precipitating *** |
1159 |
|
|
! *** downdraft calculation *** |
1160 |
|
|
|
1161 |
|
|
|
1162 |
|
|
! *** Integrate liquid water equation to find condensed water *** |
1163 |
|
|
! *** and condensed water flux *** |
1164 |
|
|
|
1165 |
|
|
|
1166 |
|
|
DO i = 1, ncum |
1167 |
|
|
jtt(i) = 2 |
1168 |
|
|
IF (ep(i,inb(i))<=0.0001) iflag(i) = 2 |
1169 |
|
|
IF (iflag(i)==0) THEN |
1170 |
|
|
lwork(i) = .TRUE. |
1171 |
|
|
ELSE |
1172 |
|
|
lwork(i) = .FALSE. |
1173 |
|
|
END IF |
1174 |
|
|
END DO |
1175 |
|
|
|
1176 |
|
|
! *** Begin downdraft loop *** |
1177 |
|
|
|
1178 |
|
|
|
1179 |
|
|
CALL zilch(wdtrain, ncum) |
1180 |
|
|
DO i = nl + 1, 1, -1 |
1181 |
|
|
|
1182 |
|
|
num1 = 0 |
1183 |
|
|
DO ij = 1, ncum |
1184 |
|
|
IF ((i<=inb(ij)) .AND. lwork(ij)) num1 = num1 + 1 |
1185 |
|
|
END DO |
1186 |
|
|
IF (num1<=0) GO TO 899 |
1187 |
|
|
|
1188 |
|
|
|
1189 |
|
|
! *** Calculate detrained precipitation *** |
1190 |
|
|
|
1191 |
|
|
DO ij = 1, ncum |
1192 |
|
|
IF ((i<=inb(ij)) .AND. (lwork(ij))) THEN |
1193 |
|
|
wdtrain(ij) = g*ep(ij, i)*m(ij, i)*clw(ij, i) |
1194 |
|
|
END IF |
1195 |
|
|
END DO |
1196 |
|
|
|
1197 |
|
|
IF (i>1) THEN |
1198 |
|
|
DO j = 1, i - 1 |
1199 |
|
|
DO ij = 1, ncum |
1200 |
|
|
IF ((i<=inb(ij)) .AND. (lwork(ij))) THEN |
1201 |
|
|
awat = elij(ij, j, i) - (1.-ep(ij,i))*clw(ij, i) |
1202 |
|
|
awat = max(0.0, awat) |
1203 |
|
|
wdtrain(ij) = wdtrain(ij) + g*awat*ment(ij, j, i) |
1204 |
|
|
END IF |
1205 |
|
|
END DO |
1206 |
|
|
END DO |
1207 |
|
|
END IF |
1208 |
|
|
|
1209 |
|
|
! *** Find rain water and evaporation using provisional *** |
1210 |
|
|
! *** estimates of qp(i)and qp(i-1) *** |
1211 |
|
|
|
1212 |
|
|
|
1213 |
|
|
! *** Value of terminal velocity and coeffecient of evaporation for snow |
1214 |
|
|
! *** |
1215 |
|
|
|
1216 |
|
|
DO ij = 1, ncum |
1217 |
|
|
IF ((i<=inb(ij)) .AND. (lwork(ij))) THEN |
1218 |
|
|
coeff = coeffs |
1219 |
|
|
wt(ij, i) = omtsnow |
1220 |
|
|
|
1221 |
|
|
! *** Value of terminal velocity and coeffecient of evaporation for |
1222 |
|
|
! rain *** |
1223 |
|
|
|
1224 |
|
|
IF (t(ij,i)>273.0) THEN |
1225 |
|
|
coeff = coeffr |
1226 |
|
|
wt(ij, i) = omtrain |
1227 |
|
|
END IF |
1228 |
|
|
qsm = 0.5*(q(ij,i)+qp(ij,i+1)) |
1229 |
|
|
afac = coeff*ph(ij, i)*(qs(ij,i)-qsm)/(1.0E4+2.0E3*ph(ij,i)*qs(ij,i)) |
1230 |
|
|
afac = max(afac, 0.0) |
1231 |
|
|
sigt = sigp(ij, i) |
1232 |
|
|
sigt = max(0.0, sigt) |
1233 |
|
|
sigt = min(1.0, sigt) |
1234 |
|
|
b6 = 100.*(ph(ij,i)-ph(ij,i+1))*sigt*afac/wt(ij, i) |
1235 |
|
|
c6 = (water(ij,i+1)*wt(ij,i+1)+wdtrain(ij)/sigd)/wt(ij, i) |
1236 |
|
|
revap = 0.5*(-b6+sqrt(b6*b6+4.*c6)) |
1237 |
|
|
evap(ij, i) = sigt*afac*revap |
1238 |
|
|
water(ij, i) = revap*revap |
1239 |
|
|
|
1240 |
|
|
! *** Calculate precipitating downdraft mass flux under *** |
1241 |
|
|
! *** hydrostatic approximation *** |
1242 |
|
|
|
1243 |
|
|
IF (i>1) THEN |
1244 |
|
|
dhdp = (h(ij,i)-h(ij,i-1))/(p(ij,i-1)-p(ij,i)) |
1245 |
|
|
dhdp = max(dhdp, 10.0) |
1246 |
|
|
mp(ij, i) = 100.*ginv*lv(ij, i)*sigd*evap(ij, i)/dhdp |
1247 |
|
|
mp(ij, i) = max(mp(ij,i), 0.0) |
1248 |
|
|
|
1249 |
|
|
! *** Add small amount of inertia to downdraft *** |
1250 |
|
|
|
1251 |
|
|
fac = 20.0/(ph(ij,i-1)-ph(ij,i)) |
1252 |
|
|
mp(ij, i) = (fac*mp(ij,i+1)+mp(ij,i))/(1.+fac) |
1253 |
|
|
|
1254 |
|
|
! *** Force mp to decrease linearly to zero |
1255 |
|
|
! *** |
1256 |
|
|
! *** between about 950 mb and the surface |
1257 |
|
|
! *** |
1258 |
|
|
|
1259 |
|
|
IF (p(ij,i)>(0.949*p(ij,1))) THEN |
1260 |
|
|
jtt(ij) = max(jtt(ij), i) |
1261 |
|
|
mp(ij, i) = mp(ij, jtt(ij))*(p(ij,1)-p(ij,i))/ & |
1262 |
|
|
(p(ij,1)-p(ij,jtt(ij))) |
1263 |
|
|
END IF |
1264 |
|
|
END IF |
1265 |
|
|
|
1266 |
|
|
! *** Find mixing ratio of precipitating downdraft *** |
1267 |
|
|
|
1268 |
|
|
IF (i/=inb(ij)) THEN |
1269 |
|
|
IF (i==1) THEN |
1270 |
|
|
qstm = qs(ij, 1) |
1271 |
|
|
ELSE |
1272 |
|
|
qstm = qs(ij, i-1) |
1273 |
|
|
END IF |
1274 |
|
|
IF (mp(ij,i)>mp(ij,i+1)) THEN |
1275 |
|
|
rat = mp(ij, i+1)/mp(ij, i) |
1276 |
|
|
qp(ij, i) = qp(ij, i+1)*rat + q(ij, i)*(1.0-rat) + & |
1277 |
|
|
100.*ginv*sigd*(ph(ij,i)-ph(ij,i+1))*(evap(ij,i)/mp(ij,i)) |
1278 |
|
|
up(ij, i) = up(ij, i+1)*rat + u(ij, i)*(1.-rat) |
1279 |
|
|
vp(ij, i) = vp(ij, i+1)*rat + v(ij, i)*(1.-rat) |
1280 |
|
|
ELSE |
1281 |
|
|
IF (mp(ij,i+1)>0.0) THEN |
1282 |
|
|
qp(ij, i) = (gz(ij,i+1)-gz(ij,i)+qp(ij,i+1)*(lv(ij,i+1)+t(ij, & |
1283 |
|
|
i+1)*(cl-cpd))+cpd*(t(ij,i+1)-t(ij, & |
1284 |
|
|
i)))/(lv(ij,i)+t(ij,i)*(cl-cpd)) |
1285 |
|
|
up(ij, i) = up(ij, i+1) |
1286 |
|
|
vp(ij, i) = vp(ij, i+1) |
1287 |
|
|
END IF |
1288 |
|
|
END IF |
1289 |
|
|
qp(ij, i) = min(qp(ij,i), qstm) |
1290 |
|
|
qp(ij, i) = max(qp(ij,i), 0.0) |
1291 |
|
|
END IF |
1292 |
|
|
END IF |
1293 |
|
|
END DO |
1294 |
|
|
899 END DO |
1295 |
|
|
|
1296 |
|
|
RETURN |
1297 |
|
|
END SUBROUTINE cv_unsat |
1298 |
|
|
|
1299 |
|
|
SUBROUTINE cv_yield(nloc, ncum, nd, nk, icb, inb, delt, t, q, u, v, gz, p, & |
1300 |
|
|
ph, h, hp, lv, cpn, ep, clw, frac, m, mp, qp, up, vp, wt, water, evap, & |
1301 |
|
|
ment, qent, uent, vent, nent, elij, tv, tvp, iflag, wd, qprime, tprime, & |
1302 |
|
|
precip, cbmf, ft, fq, fu, fv, ma, qcondc) |
1303 |
|
|
IMPLICIT NONE |
1304 |
|
|
|
1305 |
|
|
include "cvthermo.h" |
1306 |
|
|
include "cvparam.h" |
1307 |
|
|
|
1308 |
|
|
! inputs |
1309 |
|
|
INTEGER ncum, nd, nloc |
1310 |
|
|
INTEGER nk(nloc), icb(nloc), inb(nloc) |
1311 |
|
|
INTEGER nent(nloc, nd) |
1312 |
|
|
REAL delt |
1313 |
|
|
REAL t(nloc, nd), q(nloc, nd), u(nloc, nd), v(nloc, nd) |
1314 |
|
|
REAL gz(nloc, nd) |
1315 |
|
|
REAL p(nloc, nd), ph(nloc, nd+1), h(nloc, nd) |
1316 |
|
|
REAL hp(nloc, nd), lv(nloc, nd) |
1317 |
|
|
REAL cpn(nloc, nd), ep(nloc, nd), clw(nloc, nd), frac(nloc) |
1318 |
|
|
REAL m(nloc, nd), mp(nloc, nd), qp(nloc, nd) |
1319 |
|
|
REAL up(nloc, nd), vp(nloc, nd) |
1320 |
|
|
REAL wt(nloc, nd), water(nloc, nd), evap(nloc, nd) |
1321 |
|
|
REAL ment(nloc, nd, nd), qent(nloc, nd, nd), elij(nloc, nd, nd) |
1322 |
|
|
REAL uent(nloc, nd, nd), vent(nloc, nd, nd) |
1323 |
|
|
REAL tv(nloc, nd), tvp(nloc, nd) |
1324 |
|
|
|
1325 |
|
|
! outputs |
1326 |
|
|
INTEGER iflag(nloc) ! also an input |
1327 |
|
|
REAL cbmf(nloc) ! also an input |
1328 |
|
|
REAL wd(nloc), tprime(nloc), qprime(nloc) |
1329 |
|
|
REAL precip(nloc) |
1330 |
|
|
REAL ft(nloc, nd), fq(nloc, nd), fu(nloc, nd), fv(nloc, nd) |
1331 |
|
|
REAL ma(nloc, nd) |
1332 |
|
|
REAL qcondc(nloc, nd) |
1333 |
|
|
|
1334 |
|
|
! local variables |
1335 |
|
|
INTEGER i, j, ij, k, num1 |
1336 |
|
|
REAL dpinv, cpinv, awat, fqold, ftold, fuold, fvold, delti |
1337 |
|
|
REAL work(nloc), am(nloc), amp1(nloc), ad(nloc) |
1338 |
|
|
REAL ents(nloc), uav(nloc), vav(nloc), lvcp(nloc, nd) |
1339 |
|
|
REAL qcond(nloc, nd), nqcond(nloc, nd), wa(nloc, nd) ! cld |
1340 |
|
|
REAL siga(nloc, nd), ax(nloc, nd), mac(nloc, nd) ! cld |
1341 |
|
|
|
1342 |
|
|
|
1343 |
|
|
! -- initializations: |
1344 |
|
|
|
1345 |
|
|
delti = 1.0/delt |
1346 |
|
|
|
1347 |
|
|
DO i = 1, ncum |
1348 |
|
|
precip(i) = 0.0 |
1349 |
|
|
wd(i) = 0.0 |
1350 |
|
|
tprime(i) = 0.0 |
1351 |
|
|
qprime(i) = 0.0 |
1352 |
|
|
DO k = 1, nl + 1 |
1353 |
|
|
ft(i, k) = 0.0 |
1354 |
|
|
fu(i, k) = 0.0 |
1355 |
|
|
fv(i, k) = 0.0 |
1356 |
|
|
fq(i, k) = 0.0 |
1357 |
|
|
lvcp(i, k) = lv(i, k)/cpn(i, k) |
1358 |
|
|
qcondc(i, k) = 0.0 ! cld |
1359 |
|
|
qcond(i, k) = 0.0 ! cld |
1360 |
|
|
nqcond(i, k) = 0.0 ! cld |
1361 |
|
|
END DO |
1362 |
|
|
END DO |
1363 |
|
|
|
1364 |
|
|
|
1365 |
|
|
! *** Calculate surface precipitation in mm/day *** |
1366 |
|
|
|
1367 |
|
|
DO i = 1, ncum |
1368 |
|
|
IF (iflag(i)<=1) THEN |
1369 |
|
|
! c precip(i)=precip(i)+wt(i,1)*sigd*water(i,1)*3600.*24000. |
1370 |
|
|
! c & /(rowl*g) |
1371 |
|
|
! c precip(i)=precip(i)*delt/86400. |
1372 |
|
|
precip(i) = wt(i, 1)*sigd*water(i, 1)*86400/g |
1373 |
|
|
END IF |
1374 |
|
|
END DO |
1375 |
|
|
|
1376 |
|
|
|
1377 |
|
|
! *** Calculate downdraft velocity scale and surface temperature and *** |
1378 |
|
|
! *** water vapor fluctuations *** |
1379 |
|
|
|
1380 |
|
|
DO i = 1, ncum |
1381 |
|
|
wd(i) = betad*abs(mp(i,icb(i)))*0.01*rrd*t(i, icb(i))/(sigd*p(i,icb(i))) |
1382 |
|
|
qprime(i) = 0.5*(qp(i,1)-q(i,1)) |
1383 |
|
|
tprime(i) = lv0*qprime(i)/cpd |
1384 |
|
|
END DO |
1385 |
|
|
|
1386 |
|
|
! *** Calculate tendencies of lowest level potential temperature *** |
1387 |
|
|
! *** and mixing ratio *** |
1388 |
|
|
|
1389 |
|
|
DO i = 1, ncum |
1390 |
|
|
work(i) = 0.01/(ph(i,1)-ph(i,2)) |
1391 |
|
|
am(i) = 0.0 |
1392 |
|
|
END DO |
1393 |
|
|
DO k = 2, nl |
1394 |
|
|
DO i = 1, ncum |
1395 |
|
|
IF ((nk(i)==1) .AND. (k<=inb(i)) .AND. (nk(i)==1)) THEN |
1396 |
|
|
am(i) = am(i) + m(i, k) |
1397 |
|
|
END IF |
1398 |
|
|
END DO |
1399 |
|
|
END DO |
1400 |
|
|
DO i = 1, ncum |
1401 |
|
|
IF ((g*work(i)*am(i))>=delti) iflag(i) = 1 |
1402 |
|
|
ft(i, 1) = ft(i, 1) + g*work(i)*am(i)*(t(i,2)-t(i,1)+(gz(i,2)-gz(i, & |
1403 |
|
|
1))/cpn(i,1)) |
1404 |
|
|
ft(i, 1) = ft(i, 1) - lvcp(i, 1)*sigd*evap(i, 1) |
1405 |
|
|
ft(i, 1) = ft(i, 1) + sigd*wt(i, 2)*(cl-cpd)*water(i, 2)*(t(i,2)-t(i,1))* & |
1406 |
|
|
work(i)/cpn(i, 1) |
1407 |
|
|
fq(i, 1) = fq(i, 1) + g*mp(i, 2)*(qp(i,2)-q(i,1))*work(i) + & |
1408 |
|
|
sigd*evap(i, 1) |
1409 |
|
|
fq(i, 1) = fq(i, 1) + g*am(i)*(q(i,2)-q(i,1))*work(i) |
1410 |
|
|
fu(i, 1) = fu(i, 1) + g*work(i)*(mp(i,2)*(up(i,2)-u(i,1))+am(i)*(u(i, & |
1411 |
|
|
2)-u(i,1))) |
1412 |
|
|
fv(i, 1) = fv(i, 1) + g*work(i)*(mp(i,2)*(vp(i,2)-v(i,1))+am(i)*(v(i, & |
1413 |
|
|
2)-v(i,1))) |
1414 |
|
|
END DO |
1415 |
|
|
DO j = 2, nl |
1416 |
|
|
DO i = 1, ncum |
1417 |
|
|
IF (j<=inb(i)) THEN |
1418 |
|
|
fq(i, 1) = fq(i, 1) + g*work(i)*ment(i, j, 1)*(qent(i,j,1)-q(i,1)) |
1419 |
|
|
fu(i, 1) = fu(i, 1) + g*work(i)*ment(i, j, 1)*(uent(i,j,1)-u(i,1)) |
1420 |
|
|
fv(i, 1) = fv(i, 1) + g*work(i)*ment(i, j, 1)*(vent(i,j,1)-v(i,1)) |
1421 |
|
|
END IF |
1422 |
|
|
END DO |
1423 |
|
|
END DO |
1424 |
|
|
|
1425 |
|
|
! *** Calculate tendencies of potential temperature and mixing ratio *** |
1426 |
|
|
! *** at levels above the lowest level *** |
1427 |
|
|
|
1428 |
|
|
! *** First find the net saturated updraft and downdraft mass fluxes *** |
1429 |
|
|
! *** through each level *** |
1430 |
|
|
|
1431 |
|
|
DO i = 2, nl + 1 |
1432 |
|
|
|
1433 |
|
|
num1 = 0 |
1434 |
|
|
DO ij = 1, ncum |
1435 |
|
|
IF (i<=inb(ij)) num1 = num1 + 1 |
1436 |
|
|
END DO |
1437 |
|
|
IF (num1<=0) GO TO 1500 |
1438 |
|
|
|
1439 |
|
|
CALL zilch(amp1, ncum) |
1440 |
|
|
CALL zilch(ad, ncum) |
1441 |
|
|
|
1442 |
|
|
DO k = i + 1, nl + 1 |
1443 |
|
|
DO ij = 1, ncum |
1444 |
|
|
IF ((i>=nk(ij)) .AND. (i<=inb(ij)) .AND. (k<=(inb(ij)+1))) THEN |
1445 |
|
|
amp1(ij) = amp1(ij) + m(ij, k) |
1446 |
|
|
END IF |
1447 |
|
|
END DO |
1448 |
|
|
END DO |
1449 |
|
|
|
1450 |
|
|
DO k = 1, i |
1451 |
|
|
DO j = i + 1, nl + 1 |
1452 |
|
|
DO ij = 1, ncum |
1453 |
|
|
IF ((j<=(inb(ij)+1)) .AND. (i<=inb(ij))) THEN |
1454 |
|
|
amp1(ij) = amp1(ij) + ment(ij, k, j) |
1455 |
|
|
END IF |
1456 |
|
|
END DO |
1457 |
|
|
END DO |
1458 |
|
|
END DO |
1459 |
|
|
DO k = 1, i - 1 |
1460 |
|
|
DO j = i, nl + 1 |
1461 |
|
|
DO ij = 1, ncum |
1462 |
|
|
IF ((i<=inb(ij)) .AND. (j<=inb(ij))) THEN |
1463 |
|
|
ad(ij) = ad(ij) + ment(ij, j, k) |
1464 |
|
|
END IF |
1465 |
|
|
END DO |
1466 |
|
|
END DO |
1467 |
|
|
END DO |
1468 |
|
|
|
1469 |
|
|
DO ij = 1, ncum |
1470 |
|
|
IF (i<=inb(ij)) THEN |
1471 |
|
|
dpinv = 0.01/(ph(ij,i)-ph(ij,i+1)) |
1472 |
|
|
cpinv = 1.0/cpn(ij, i) |
1473 |
|
|
|
1474 |
|
|
ft(ij, i) = ft(ij, i) + g*dpinv*(amp1(ij)*(t(ij,i+1)-t(ij, & |
1475 |
|
|
i)+(gz(ij,i+1)-gz(ij,i))*cpinv)-ad(ij)*(t(ij,i)-t(ij, & |
1476 |
|
|
i-1)+(gz(ij,i)-gz(ij,i-1))*cpinv)) - sigd*lvcp(ij, i)*evap(ij, i) |
1477 |
|
|
ft(ij, i) = ft(ij, i) + g*dpinv*ment(ij, i, i)*(hp(ij,i)-h(ij,i)+t(ij & |
1478 |
|
|
,i)*(cpv-cpd)*(q(ij,i)-qent(ij,i,i)))*cpinv |
1479 |
|
|
ft(ij, i) = ft(ij, i) + sigd*wt(ij, i+1)*(cl-cpd)*water(ij, i+1)*(t( & |
1480 |
|
|
ij,i+1)-t(ij,i))*dpinv*cpinv |
1481 |
|
|
fq(ij, i) = fq(ij, i) + g*dpinv*(amp1(ij)*(q(ij,i+1)-q(ij, & |
1482 |
|
|
i))-ad(ij)*(q(ij,i)-q(ij,i-1))) |
1483 |
|
|
fu(ij, i) = fu(ij, i) + g*dpinv*(amp1(ij)*(u(ij,i+1)-u(ij, & |
1484 |
|
|
i))-ad(ij)*(u(ij,i)-u(ij,i-1))) |
1485 |
|
|
fv(ij, i) = fv(ij, i) + g*dpinv*(amp1(ij)*(v(ij,i+1)-v(ij, & |
1486 |
|
|
i))-ad(ij)*(v(ij,i)-v(ij,i-1))) |
1487 |
|
|
END IF |
1488 |
|
|
END DO |
1489 |
|
|
DO k = 1, i - 1 |
1490 |
|
|
DO ij = 1, ncum |
1491 |
|
|
IF (i<=inb(ij)) THEN |
1492 |
|
|
awat = elij(ij, k, i) - (1.-ep(ij,i))*clw(ij, i) |
1493 |
|
|
awat = max(awat, 0.0) |
1494 |
|
|
fq(ij, i) = fq(ij, i) + g*dpinv*ment(ij, k, i)*(qent(ij,k,i)-awat-q & |
1495 |
|
|
(ij,i)) |
1496 |
|
|
fu(ij, i) = fu(ij, i) + g*dpinv*ment(ij, k, i)*(uent(ij,k,i)-u(ij,i & |
1497 |
|
|
)) |
1498 |
|
|
fv(ij, i) = fv(ij, i) + g*dpinv*ment(ij, k, i)*(vent(ij,k,i)-v(ij,i & |
1499 |
|
|
)) |
1500 |
|
|
! (saturated updrafts resulting from mixing) ! cld |
1501 |
|
|
qcond(ij, i) = qcond(ij, i) + (elij(ij,k,i)-awat) ! cld |
1502 |
|
|
nqcond(ij, i) = nqcond(ij, i) + 1. ! cld |
1503 |
|
|
END IF |
1504 |
|
|
END DO |
1505 |
|
|
END DO |
1506 |
|
|
DO k = i, nl + 1 |
1507 |
|
|
DO ij = 1, ncum |
1508 |
|
|
IF ((i<=inb(ij)) .AND. (k<=inb(ij))) THEN |
1509 |
|
|
fq(ij, i) = fq(ij, i) + g*dpinv*ment(ij, k, i)*(qent(ij,k,i)-q(ij,i & |
1510 |
|
|
)) |
1511 |
|
|
fu(ij, i) = fu(ij, i) + g*dpinv*ment(ij, k, i)*(uent(ij,k,i)-u(ij,i & |
1512 |
|
|
)) |
1513 |
|
|
fv(ij, i) = fv(ij, i) + g*dpinv*ment(ij, k, i)*(vent(ij,k,i)-v(ij,i & |
1514 |
|
|
)) |
1515 |
|
|
END IF |
1516 |
|
|
END DO |
1517 |
|
|
END DO |
1518 |
|
|
DO ij = 1, ncum |
1519 |
|
|
IF (i<=inb(ij)) THEN |
1520 |
|
|
fq(ij, i) = fq(ij, i) + sigd*evap(ij, i) + g*(mp(ij,i+1)*(qp(ij, & |
1521 |
|
|
i+1)-q(ij,i))-mp(ij,i)*(qp(ij,i)-q(ij,i-1)))*dpinv |
1522 |
|
|
fu(ij, i) = fu(ij, i) + g*(mp(ij,i+1)*(up(ij,i+1)-u(ij, & |
1523 |
|
|
i))-mp(ij,i)*(up(ij,i)-u(ij,i-1)))*dpinv |
1524 |
|
|
fv(ij, i) = fv(ij, i) + g*(mp(ij,i+1)*(vp(ij,i+1)-v(ij, & |
1525 |
|
|
i))-mp(ij,i)*(vp(ij,i)-v(ij,i-1)))*dpinv |
1526 |
|
|
! (saturated downdrafts resulting from mixing) ! cld |
1527 |
|
|
DO k = i + 1, inb(ij) ! cld |
1528 |
|
|
qcond(ij, i) = qcond(ij, i) + elij(ij, k, i) ! cld |
1529 |
|
|
nqcond(ij, i) = nqcond(ij, i) + 1. ! cld |
1530 |
|
|
END DO ! cld |
1531 |
|
|
! (particular case: no detraining level is found) ! cld |
1532 |
|
|
IF (nent(ij,i)==0) THEN ! cld |
1533 |
|
|
qcond(ij, i) = qcond(ij, i) + (1.-ep(ij,i))*clw(ij, i) ! cld |
1534 |
|
|
nqcond(ij, i) = nqcond(ij, i) + 1. ! cld |
1535 |
|
|
END IF ! cld |
1536 |
|
|
IF (nqcond(ij,i)/=0.) THEN ! cld |
1537 |
|
|
qcond(ij, i) = qcond(ij, i)/nqcond(ij, i) ! cld |
1538 |
|
|
END IF ! cld |
1539 |
|
|
END IF |
1540 |
|
|
END DO |
1541 |
|
|
1500 END DO |
1542 |
|
|
|
1543 |
|
|
! *** Adjust tendencies at top of convection layer to reflect *** |
1544 |
|
|
! *** actual position of the level zero cape *** |
1545 |
|
|
|
1546 |
|
|
DO ij = 1, ncum |
1547 |
|
|
fqold = fq(ij, inb(ij)) |
1548 |
|
|
fq(ij, inb(ij)) = fq(ij, inb(ij))*(1.-frac(ij)) |
1549 |
|
|
fq(ij, inb(ij)-1) = fq(ij, inb(ij)-1) + frac(ij)*fqold*((ph(ij, & |
1550 |
|
|
inb(ij))-ph(ij,inb(ij)+1))/(ph(ij,inb(ij)-1)-ph(ij, & |
1551 |
|
|
inb(ij))))*lv(ij, inb(ij))/lv(ij, inb(ij)-1) |
1552 |
|
|
ftold = ft(ij, inb(ij)) |
1553 |
|
|
ft(ij, inb(ij)) = ft(ij, inb(ij))*(1.-frac(ij)) |
1554 |
|
|
ft(ij, inb(ij)-1) = ft(ij, inb(ij)-1) + frac(ij)*ftold*((ph(ij, & |
1555 |
|
|
inb(ij))-ph(ij,inb(ij)+1))/(ph(ij,inb(ij)-1)-ph(ij, & |
1556 |
|
|
inb(ij))))*cpn(ij, inb(ij))/cpn(ij, inb(ij)-1) |
1557 |
|
|
fuold = fu(ij, inb(ij)) |
1558 |
|
|
fu(ij, inb(ij)) = fu(ij, inb(ij))*(1.-frac(ij)) |
1559 |
|
|
fu(ij, inb(ij)-1) = fu(ij, inb(ij)-1) + frac(ij)*fuold*((ph(ij, & |
1560 |
|
|
inb(ij))-ph(ij,inb(ij)+1))/(ph(ij,inb(ij)-1)-ph(ij,inb(ij)))) |
1561 |
|
|
fvold = fv(ij, inb(ij)) |
1562 |
|
|
fv(ij, inb(ij)) = fv(ij, inb(ij))*(1.-frac(ij)) |
1563 |
|
|
fv(ij, inb(ij)-1) = fv(ij, inb(ij)-1) + frac(ij)*fvold*((ph(ij, & |
1564 |
|
|
inb(ij))-ph(ij,inb(ij)+1))/(ph(ij,inb(ij)-1)-ph(ij,inb(ij)))) |
1565 |
|
|
END DO |
1566 |
|
|
|
1567 |
|
|
! *** Very slightly adjust tendencies to force exact *** |
1568 |
|
|
! *** enthalpy, momentum and tracer conservation *** |
1569 |
|
|
|
1570 |
|
|
DO ij = 1, ncum |
1571 |
|
|
ents(ij) = 0.0 |
1572 |
|
|
uav(ij) = 0.0 |
1573 |
|
|
vav(ij) = 0.0 |
1574 |
|
|
DO i = 1, inb(ij) |
1575 |
|
|
ents(ij) = ents(ij) + (cpn(ij,i)*ft(ij,i)+lv(ij,i)*fq(ij,i))*(ph(ij,i)- & |
1576 |
|
|
ph(ij,i+1)) |
1577 |
|
|
uav(ij) = uav(ij) + fu(ij, i)*(ph(ij,i)-ph(ij,i+1)) |
1578 |
|
|
vav(ij) = vav(ij) + fv(ij, i)*(ph(ij,i)-ph(ij,i+1)) |
1579 |
|
|
END DO |
1580 |
|
|
END DO |
1581 |
|
|
DO ij = 1, ncum |
1582 |
|
|
ents(ij) = ents(ij)/(ph(ij,1)-ph(ij,inb(ij)+1)) |
1583 |
|
|
uav(ij) = uav(ij)/(ph(ij,1)-ph(ij,inb(ij)+1)) |
1584 |
|
|
vav(ij) = vav(ij)/(ph(ij,1)-ph(ij,inb(ij)+1)) |
1585 |
|
|
END DO |
1586 |
|
|
DO ij = 1, ncum |
1587 |
|
|
DO i = 1, inb(ij) |
1588 |
|
|
ft(ij, i) = ft(ij, i) - ents(ij)/cpn(ij, i) |
1589 |
|
|
fu(ij, i) = (1.-cu)*(fu(ij,i)-uav(ij)) |
1590 |
|
|
fv(ij, i) = (1.-cu)*(fv(ij,i)-vav(ij)) |
1591 |
|
|
END DO |
1592 |
|
|
END DO |
1593 |
|
|
|
1594 |
|
|
DO k = 1, nl + 1 |
1595 |
|
|
DO i = 1, ncum |
1596 |
|
|
IF ((q(i,k)+delt*fq(i,k))<0.0) iflag(i) = 10 |
1597 |
|
|
END DO |
1598 |
|
|
END DO |
1599 |
|
|
|
1600 |
|
|
|
1601 |
|
|
DO i = 1, ncum |
1602 |
|
|
IF (iflag(i)>2) THEN |
1603 |
|
|
precip(i) = 0.0 |
1604 |
|
|
cbmf(i) = 0.0 |
1605 |
|
|
END IF |
1606 |
|
|
END DO |
1607 |
|
|
DO k = 1, nl |
1608 |
|
|
DO i = 1, ncum |
1609 |
|
|
IF (iflag(i)>2) THEN |
1610 |
|
|
ft(i, k) = 0.0 |
1611 |
|
|
fq(i, k) = 0.0 |
1612 |
|
|
fu(i, k) = 0.0 |
1613 |
|
|
fv(i, k) = 0.0 |
1614 |
|
|
qcondc(i, k) = 0.0 ! cld |
1615 |
|
|
END IF |
1616 |
|
|
END DO |
1617 |
|
|
END DO |
1618 |
|
|
|
1619 |
|
|
DO k = 1, nl + 1 |
1620 |
|
|
DO i = 1, ncum |
1621 |
|
|
ma(i, k) = 0. |
1622 |
|
|
END DO |
1623 |
|
|
END DO |
1624 |
|
|
DO k = nl, 1, -1 |
1625 |
|
|
DO i = 1, ncum |
1626 |
|
|
ma(i, k) = ma(i, k+1) + m(i, k) |
1627 |
|
|
END DO |
1628 |
|
|
END DO |
1629 |
|
|
|
1630 |
|
|
|
1631 |
|
|
! *** diagnose the in-cloud mixing ratio *** ! cld |
1632 |
|
|
! *** of condensed water *** ! cld |
1633 |
|
|
! ! cld |
1634 |
|
|
DO ij = 1, ncum ! cld |
1635 |
|
|
DO i = 1, nd ! cld |
1636 |
|
|
mac(ij, i) = 0.0 ! cld |
1637 |
|
|
wa(ij, i) = 0.0 ! cld |
1638 |
|
|
siga(ij, i) = 0.0 ! cld |
1639 |
|
|
END DO ! cld |
1640 |
|
|
DO i = nk(ij), inb(ij) ! cld |
1641 |
|
|
DO k = i + 1, inb(ij) + 1 ! cld |
1642 |
|
|
mac(ij, i) = mac(ij, i) + m(ij, k) ! cld |
1643 |
|
|
END DO ! cld |
1644 |
|
|
END DO ! cld |
1645 |
|
|
DO i = icb(ij), inb(ij) - 1 ! cld |
1646 |
|
|
ax(ij, i) = 0. ! cld |
1647 |
|
|
DO j = icb(ij), i ! cld |
1648 |
|
|
ax(ij, i) = ax(ij, i) + rrd*(tvp(ij,j)-tv(ij,j)) & ! cld |
1649 |
|
|
*(ph(ij,j)-ph(ij,j+1))/p(ij, j) ! cld |
1650 |
|
|
END DO ! cld |
1651 |
|
|
IF (ax(ij,i)>0.0) THEN ! cld |
1652 |
|
|
wa(ij, i) = sqrt(2.*ax(ij,i)) ! cld |
1653 |
|
|
END IF ! cld |
1654 |
|
|
END DO ! cld |
1655 |
|
|
DO i = 1, nl ! cld |
1656 |
|
|
IF (wa(ij,i)>0.0) & ! cld |
1657 |
|
|
siga(ij, i) = mac(ij, i)/wa(ij, i) & ! cld |
1658 |
|
|
*rrd*tvp(ij, i)/p(ij, i)/100./delta ! cld |
1659 |
|
|
siga(ij, i) = min(siga(ij,i), 1.0) ! cld |
1660 |
|
|
qcondc(ij, i) = siga(ij, i)*clw(ij, i)*(1.-ep(ij,i)) & ! cld |
1661 |
|
|
+(1.-siga(ij,i))*qcond(ij, i) ! cld |
1662 |
|
|
END DO ! cld |
1663 |
|
|
END DO ! cld |
1664 |
|
|
|
1665 |
|
|
RETURN |
1666 |
|
|
END SUBROUTINE cv_yield |
1667 |
|
|
|
1668 |
|
|
SUBROUTINE cv_uncompress(nloc, len, ncum, nd, idcum, iflag, precip, cbmf, ft, & |
1669 |
|
|
fq, fu, fv, ma, qcondc, iflag1, precip1, cbmf1, ft1, fq1, fu1, fv1, ma1, & |
1670 |
|
|
qcondc1) |
1671 |
|
|
IMPLICIT NONE |
1672 |
|
|
|
1673 |
|
|
include "cvparam.h" |
1674 |
|
|
|
1675 |
|
|
! inputs: |
1676 |
|
|
INTEGER len, ncum, nd, nloc |
1677 |
|
|
INTEGER idcum(nloc) |
1678 |
|
|
INTEGER iflag(nloc) |
1679 |
|
|
REAL precip(nloc), cbmf(nloc) |
1680 |
|
|
REAL ft(nloc, nd), fq(nloc, nd), fu(nloc, nd), fv(nloc, nd) |
1681 |
|
|
REAL ma(nloc, nd) |
1682 |
|
|
REAL qcondc(nloc, nd) !cld |
1683 |
|
|
|
1684 |
|
|
! outputs: |
1685 |
|
|
INTEGER iflag1(len) |
1686 |
|
|
REAL precip1(len), cbmf1(len) |
1687 |
|
|
REAL ft1(len, nd), fq1(len, nd), fu1(len, nd), fv1(len, nd) |
1688 |
|
|
REAL ma1(len, nd) |
1689 |
|
|
REAL qcondc1(len, nd) !cld |
1690 |
|
|
|
1691 |
|
|
! local variables: |
1692 |
|
|
INTEGER i, k |
1693 |
|
|
|
1694 |
|
|
DO i = 1, ncum |
1695 |
|
|
precip1(idcum(i)) = precip(i) |
1696 |
|
|
cbmf1(idcum(i)) = cbmf(i) |
1697 |
|
|
iflag1(idcum(i)) = iflag(i) |
1698 |
|
|
END DO |
1699 |
|
|
|
1700 |
|
|
DO k = 1, nl |
1701 |
|
|
DO i = 1, ncum |
1702 |
|
|
ft1(idcum(i), k) = ft(i, k) |
1703 |
|
|
fq1(idcum(i), k) = fq(i, k) |
1704 |
|
|
fu1(idcum(i), k) = fu(i, k) |
1705 |
|
|
fv1(idcum(i), k) = fv(i, k) |
1706 |
|
|
ma1(idcum(i), k) = ma(i, k) |
1707 |
|
|
qcondc1(idcum(i), k) = qcondc(i, k) |
1708 |
|
|
END DO |
1709 |
|
|
END DO |
1710 |
|
|
|
1711 |
|
|
RETURN |
1712 |
|
|
END SUBROUTINE cv_uncompress |
1713 |
|
|
|