/[lmdze]/trunk/Sources/phylmd/CV30_routines/cv30_yield.f
ViewVC logotype

Contents of /trunk/Sources/phylmd/CV30_routines/cv30_yield.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 200 - (show annotations)
Thu Jun 2 15:40:30 2016 UTC (7 years, 11 months ago) by guez
File size: 19036 byte(s)
Changes results.
1 module cv30_yield_m
2
3 implicit none
4
5 contains
6
7 SUBROUTINE cv30_yield(icb, inb, delt, t, rr, u, v, gz, p, ph, h, hp, lv, &
8 cpn, th, ep, clw, m, tp, mp, rp, up, vp, wt, water, evap, b, ment, &
9 qent, uent, vent, nent, elij, sig, tv, tvp, iflag, precip, VPrecip, &
10 ft, fr, fu, fv, upwd, dnwd, dnwd0, ma, mike, tls, tps, qcondc)
11
12 ! Tendencies, precipitation, variables of interface with other
13 ! processes, etc.
14
15 use conema3_m, only: iflag_clw
16 use cv30_param_m, only: minorig, nl, sigd
17 use cv_thermo_m, only: cl, cpd, cpv, rowl, rrd, rrv
18 USE dimphy, ONLY: klev, klon
19 use SUPHEC_M, only: rg
20
21 ! inputs:
22 integer, intent(in):: icb(:), inb(:) ! (ncum)
23 real, intent(in):: delt
24 real, intent(in):: t(klon, klev), rr(klon, klev)
25 real, intent(in):: u(klon, klev), v(klon, klev)
26 real gz(klon, klev)
27 real p(klon, klev)
28 real ph(klon, klev + 1), h(klon, klev), hp(klon, klev)
29 real lv(klon, klev), cpn(klon, klev)
30 real th(klon, klev)
31 real ep(klon, klev), clw(klon, klev)
32 real m(klon, klev)
33 real tp(klon, klev)
34 real mp(klon, klev), rp(klon, klev), up(klon, klev)
35 real, intent(in):: vp(:, 2:) ! (ncum, 2:nl)
36 real, intent(in):: wt(:, :) ! (ncum, nl - 1)
37 real, intent(in):: water(:, :), evap(:, :) ! (ncum, nl)
38 real, intent(in):: b(:, :) ! (ncum, nl - 1)
39 real ment(klon, klev, klev), qent(klon, klev, klev), uent(klon, klev, klev)
40 real vent(klon, klev, klev)
41 integer nent(klon, klev)
42 real elij(klon, klev, klev)
43 real sig(klon, klev)
44 real tv(klon, klev), tvp(klon, klev)
45
46 ! outputs:
47 integer, intent(out):: iflag(:) ! (ncum)
48 real precip(klon)
49 real VPrecip(klon, klev + 1)
50 real ft(klon, klev), fr(klon, klev), fu(klon, klev), fv(klon, klev)
51 real upwd(klon, klev), dnwd(klon, klev)
52 real dnwd0(klon, klev)
53 real ma(klon, klev)
54 real mike(klon, klev)
55 real tls(klon, klev), tps(klon, klev)
56 real qcondc(klon, klev)
57
58 ! Local:
59 real, parameter:: delta = 0.01 ! interface cloud parameterization
60 integer ncum
61 integer i, k, il, n, j
62 real awat, delti
63 real ax, bx, cx, dx
64 real cpinv, rdcp, dpinv
65 real lvcp(klon, klev)
66 real am(klon), work(klon), ad(klon), amp1(klon)
67 real up1(klon, klev, klev), dn1(klon, klev, klev)
68 real asum(klon), bsum(klon), csum(klon), dsum(klon)
69 real qcond(klon, klev), nqcond(klon, klev), wa(klon, klev)
70 real siga(klon, klev), sax(klon, klev), mac(klon, klev)
71
72 !-------------------------------------------------------------
73
74 ncum = size(icb)
75 iflag = 0
76
77 ! initialization:
78
79 delti = 1.0 / delt
80
81 do il = 1, ncum
82 precip(il) = 0.0
83 VPrecip(il, klev + 1) = 0.
84 enddo
85
86 do i = 1, klev
87 do il = 1, ncum
88 VPrecip(il, i) = 0.0
89 ft(il, i) = 0.0
90 fr(il, i) = 0.0
91 fu(il, i) = 0.0
92 fv(il, i) = 0.0
93 qcondc(il, i) = 0.0
94 qcond(il, i) = 0.0
95 nqcond(il, i) = 0.0
96 enddo
97 enddo
98
99 do i = 1, nl
100 do il = 1, ncum
101 lvcp(il, i) = lv(il, i) / cpn(il, i)
102 enddo
103 enddo
104
105 ! calculate surface precipitation in mm / day
106
107 do il = 1, ncum
108 if (ep(il, inb(il)) >= 1e-4) precip(il) = wt(il, 1) * sigd &
109 * water(il, 1) * 86400. * 1000. / (rowl * rg)
110 enddo
111
112 ! CALCULATE VERTICAL PROFILE OF PRECIPITATIONs IN kg / m2 / s ===
113
114 ! MAF rajout pour lessivage
115 do k = 1, nl - 1
116 do il = 1, ncum
117 if (k <= inb(il)) VPrecip(il, k) = wt(il, k) * sigd * water(il, k) &
118 / rg
119 end do
120 end do
121
122 ! calculate tendencies of lowest level potential temperature
123 ! and mixing ratio
124
125 do il = 1, ncum
126 work(il) = 1.0 / (ph(il, 1) - ph(il, 2))
127 am(il) = 0.0
128 enddo
129
130 do k = 2, nl
131 do il = 1, ncum
132 if (k <= inb(il)) am(il) = am(il) + m(il, k)
133 enddo
134 enddo
135
136 do il = 1, ncum
137 if (0.01 * rg * work(il) * am(il) >= delti) iflag(il) = 1
138
139 ft(il, 1) = 0.01 * rg * work(il) * am(il) * (t(il, 2) - t(il, 1) &
140 + (gz(il, 2) - gz(il, 1)) / cpn(il, 1)) - 0.5 * lvcp(il, 1) &
141 * sigd * (evap(il, 1) + evap(il, 2)) - 0.009 * rg * sigd &
142 * mp(il, 2) * t(il, 1) * b(il, 1) * work(il) + 0.01 * sigd &
143 * wt(il, 1) * (cl - cpd) * water(il, 2) * (t(il, 2) - t(il, 1)) &
144 * work(il) / cpn(il, 1)
145
146 !jyg1 Correction pour mieux conserver l'eau (conformite avec CONVECT4.3)
147 ! (sb: pour l'instant, on ne fait que le chgt concernant rg, pas evap)
148 fr(il, 1) = 0.01 * rg * mp(il, 2) * (rp(il, 2) - rr(il, 1)) &
149 * work(il) + sigd * 0.5 * (evap(il, 1) + evap(il, 2))
150 ! + tard : + sigd * evap(il, 1)
151
152 fr(il, 1) = fr(il, 1) + 0.01 * rg * am(il) * (rr(il, 2) - rr(il, 1)) &
153 * work(il)
154
155 fu(il, 1) = fu(il, 1) + 0.01 * rg * work(il) * (mp(il, 2) &
156 * (up(il, 2) - u(il, 1)) + am(il) * (u(il, 2) - u(il, 1)))
157 fv(il, 1) = fv(il, 1) + 0.01 * rg * work(il) * (mp(il, 2) &
158 * (vp(il, 2) - v(il, 1)) + am(il) * (v(il, 2) - v(il, 1)))
159 enddo
160
161 do j = 2, nl
162 do il = 1, ncum
163 if (j <= inb(il)) then
164 fr(il, 1) = fr(il, 1) + 0.01 * rg * work(il) * ment(il, j, 1) &
165 * (qent(il, j, 1) - rr(il, 1))
166 fu(il, 1) = fu(il, 1) + 0.01 * rg * work(il) * ment(il, j, 1) &
167 * (uent(il, j, 1) - u(il, 1))
168 fv(il, 1) = fv(il, 1) + 0.01 * rg * work(il) * ment(il, j, 1) &
169 * (vent(il, j, 1) - v(il, 1))
170 endif
171 enddo
172 enddo
173
174 ! calculate tendencies of potential temperature and mixing ratio
175 ! at levels above the lowest level
176
177 ! first find the net saturated updraft and downdraft mass fluxes
178 ! through each level
179
180 loop_i: do i = 2, nl - 1
181 if (any(inb >= i)) then
182 amp1(:ncum) = 0.
183 ad(:ncum) = 0.
184
185 do k = i + 1, nl + 1
186 do il = 1, ncum
187 if (i <= inb(il) .and. k <= (inb(il) + 1)) then
188 amp1(il) = amp1(il) + m(il, k)
189 endif
190 end do
191 end do
192
193 do k = 1, i
194 do j = i + 1, nl + 1
195 do il = 1, ncum
196 if (i <= inb(il) .and. j <= (inb(il) + 1)) then
197 amp1(il) = amp1(il) + ment(il, k, j)
198 endif
199 end do
200 end do
201 end do
202
203 do k = 1, i - 1
204 do j = i, nl + 1 ! newvecto: nl au lieu nl + 1?
205 do il = 1, ncum
206 if (i <= inb(il) .and. j <= inb(il)) then
207 ad(il) = ad(il) + ment(il, j, k)
208 endif
209 end do
210 end do
211 end do
212
213 do il = 1, ncum
214 if (i <= inb(il)) then
215 dpinv = 1.0 / (ph(il, i) - ph(il, i + 1))
216 cpinv = 1.0 / cpn(il, i)
217
218 if (0.01 * rg * dpinv * amp1(il) >= delti) iflag(il) = 1
219
220 ft(il, i) = 0.01 * rg * dpinv * (amp1(il) * (t(il, i + 1) &
221 - t(il, i) + (gz(il, i + 1) - gz(il, i)) * cpinv) &
222 - ad(il) * (t(il, i) - t(il, i - 1) + (gz(il, i) &
223 - gz(il, i - 1)) * cpinv)) - 0.5 * sigd * lvcp(il, i) &
224 * (evap(il, i) + evap(il, i + 1)) - 0.009 * rg * sigd &
225 * (mp(il, i + 1) * t(il, i) * b(il, i) - mp(il, i) &
226 * t(il, i - 1) * cpn(il, i - 1) * cpinv * b(il, i - 1)) &
227 * dpinv + 0.01 * rg * dpinv * ment(il, i, i) &
228 * (hp(il, i) - h(il, i) + t(il, i) * (cpv - cpd) &
229 * (rr(il, i) - qent(il, i, i))) * cpinv + 0.01 * sigd &
230 * wt(il, i) * (cl - cpd) * water(il, i + 1) &
231 * (t(il, i + 1) - t(il, i)) * dpinv * cpinv
232 fr(il, i) = 0.01 * rg * dpinv * (amp1(il) * (rr(il, i + 1) &
233 - rr(il, i)) - ad(il) * (rr(il, i) - rr(il, i - 1)))
234 fu(il, i) = fu(il, i) + 0.01 * rg * dpinv * (amp1(il) &
235 * (u(il, i + 1) - u(il, i)) - ad(il) * (u(il, i) &
236 - u(il, i - 1)))
237 fv(il, i) = fv(il, i) + 0.01 * rg * dpinv * (amp1(il) &
238 * (v(il, i + 1) - v(il, i)) - ad(il) * (v(il, i) &
239 - v(il, i - 1)))
240 endif
241 end do
242
243 do k = 1, i - 1
244 do il = 1, ncum
245 if (i <= inb(il)) then
246 dpinv = 1.0 / (ph(il, i) - ph(il, i + 1))
247 cpinv = 1.0 / cpn(il, i)
248
249 awat = elij(il, k, i) - (1. - ep(il, i)) * clw(il, i)
250 awat = amax1(awat, 0.0)
251
252 fr(il, i) = fr(il, i) + 0.01 * rg * dpinv &
253 * ment(il, k, i) * (qent(il, k, i) - awat - rr(il, i))
254 fu(il, i) = fu(il, i) + 0.01 * rg * dpinv &
255 * ment(il, k, i) * (uent(il, k, i) - u(il, i))
256 fv(il, i) = fv(il, i) + 0.01 * rg * dpinv &
257 * ment(il, k, i) * (vent(il, k, i) - v(il, i))
258
259 ! (saturated updrafts resulting from mixing)
260 qcond(il, i) = qcond(il, i) + (elij(il, k, i) - awat)
261 nqcond(il, i) = nqcond(il, i) + 1.
262 endif ! i
263 end do
264 end do
265
266 do k = i, nl + 1
267 do il = 1, ncum
268 if (i <= inb(il) .and. k <= inb(il)) then
269 dpinv = 1.0 / (ph(il, i) - ph(il, i + 1))
270 cpinv = 1.0 / cpn(il, i)
271
272 fr(il, i) = fr(il, i) + 0.01 * rg * dpinv &
273 * ment(il, k, i) * (qent(il, k, i) - rr(il, i))
274 fu(il, i) = fu(il, i) + 0.01 * rg * dpinv &
275 * ment(il, k, i) * (uent(il, k, i) - u(il, i))
276 fv(il, i) = fv(il, i) + 0.01 * rg * dpinv &
277 * ment(il, k, i) * (vent(il, k, i) - v(il, i))
278 endif
279 end do
280 end do
281
282 do il = 1, ncum
283 if (i <= inb(il)) then
284 dpinv = 1.0 / (ph(il, i) - ph(il, i + 1))
285 cpinv = 1.0 / cpn(il, i)
286
287 ! sb: on ne fait pas encore la correction permettant de mieux
288 ! conserver l'eau:
289 fr(il, i) = fr(il, i) + 0.5 * sigd * (evap(il, i) &
290 + evap(il, i + 1)) + 0.01 * rg * (mp(il, i + 1) &
291 * (rp(il, i + 1) - rr(il, i)) - mp(il, i) * (rp(il, i) &
292 - rr(il, i - 1))) * dpinv
293
294 fu(il, i) = fu(il, i) + 0.01 * rg * (mp(il, i + 1) &
295 * (up(il, i + 1) - u(il, i)) - mp(il, i) * (up(il, i) &
296 - u(il, i - 1))) * dpinv
297 fv(il, i) = fv(il, i) + 0.01 * rg * (mp(il, i + 1) &
298 * (vp(il, i + 1) - v(il, i)) - mp(il, i) * (vp(il, i) &
299 - v(il, i - 1))) * dpinv
300 endif
301 end do
302
303 ! sb: interface with the cloud parameterization:
304
305 do k = i + 1, nl
306 do il = 1, ncum
307 if (k <= inb(il) .and. i <= inb(il)) then
308 ! (saturated downdrafts resulting from mixing)
309 qcond(il, i) = qcond(il, i) + elij(il, k, i)
310 nqcond(il, i) = nqcond(il, i) + 1.
311 endif
312 enddo
313 enddo
314
315 ! (particular case: no detraining level is found)
316 do il = 1, ncum
317 if (i <= inb(il) .and. nent(il, i) == 0) then
318 qcond(il, i) = qcond(il, i) + (1. - ep(il, i)) * clw(il, i)
319 nqcond(il, i) = nqcond(il, i) + 1.
320 endif
321 enddo
322
323 do il = 1, ncum
324 if (i <= inb(il) .and. nqcond(il, i) /= 0.) then
325 qcond(il, i) = qcond(il, i) / nqcond(il, i)
326 endif
327 enddo
328 end if
329 end do loop_i
330
331 ! move the detrainment at level inb down to level inb - 1
332 ! in such a way as to preserve the vertically
333 ! integrated enthalpy and water tendencies
334
335 do il = 1, ncum
336 ax = 0.1 * ment(il, inb(il), inb(il)) * (hp(il, inb(il)) &
337 - h(il, inb(il)) + t(il, inb(il)) * (cpv - cpd) &
338 * (rr(il, inb(il)) - qent(il, inb(il), inb(il)))) &
339 / (cpn(il, inb(il)) * (ph(il, inb(il)) - ph(il, inb(il) + 1)))
340 ft(il, inb(il)) = ft(il, inb(il)) - ax
341 ft(il, inb(il) - 1) = ft(il, inb(il) - 1) + ax * cpn(il, inb(il)) &
342 * (ph(il, inb(il)) - ph(il, inb(il) + 1)) / (cpn(il, inb(il) - 1) &
343 * (ph(il, inb(il) - 1) - ph(il, inb(il))))
344
345 bx = 0.1 * ment(il, inb(il), inb(il)) * (qent(il, inb(il), inb(il)) &
346 - rr(il, inb(il))) / (ph(il, inb(il)) - ph(il, inb(il) + 1))
347 fr(il, inb(il)) = fr(il, inb(il)) - bx
348 fr(il, inb(il) - 1) = fr(il, inb(il) - 1) &
349 + bx * (ph(il, inb(il)) - ph(il, inb(il) + 1)) &
350 / (ph(il, inb(il) - 1) - ph(il, inb(il)))
351
352 cx = 0.1 * ment(il, inb(il), inb(il)) * (uent(il, inb(il), inb(il)) &
353 - u(il, inb(il))) / (ph(il, inb(il)) - ph(il, inb(il) + 1))
354 fu(il, inb(il)) = fu(il, inb(il)) - cx
355 fu(il, inb(il) - 1) = fu(il, inb(il) - 1) &
356 + cx * (ph(il, inb(il)) - ph(il, inb(il) + 1)) &
357 / (ph(il, inb(il) - 1) - ph(il, inb(il)))
358
359 dx = 0.1 * ment(il, inb(il), inb(il)) * (vent(il, inb(il), inb(il)) &
360 - v(il, inb(il))) / (ph(il, inb(il)) - ph(il, inb(il) + 1))
361 fv(il, inb(il)) = fv(il, inb(il)) - dx
362 fv(il, inb(il) - 1) = fv(il, inb(il) - 1) &
363 + dx * (ph(il, inb(il)) - ph(il, inb(il) + 1)) &
364 / (ph(il, inb(il) - 1) - ph(il, inb(il)))
365
366 end do
367
368 ! homoginize tendencies below cloud base
369
370 do il = 1, ncum
371 asum(il) = 0.0
372 bsum(il) = 0.0
373 csum(il) = 0.0
374 dsum(il) = 0.0
375 enddo
376
377 do i = 1, nl
378 do il = 1, ncum
379 if (i <= (icb(il) - 1)) then
380 asum(il) = asum(il) + ft(il, i) * (ph(il, i) - ph(il, i + 1))
381 bsum(il) = bsum(il) + fr(il, i) * (lv(il, i) + (cl - cpd) &
382 * (t(il, i) - t(il, 1))) * (ph(il, i) - ph(il, i + 1))
383 csum(il) = csum(il) + (lv(il, i) + (cl - cpd) * (t(il, i) &
384 - t(il, 1))) * (ph(il, i) - ph(il, i + 1))
385 dsum(il) = dsum(il) + t(il, i) * (ph(il, i) - ph(il, i + 1)) &
386 / th(il, i)
387 endif
388 enddo
389 enddo
390
391 do i = 1, nl
392 do il = 1, ncum
393 if (i <= (icb(il) - 1)) then
394 ft(il, i) = asum(il) * t(il, i) / (th(il, i) * dsum(il))
395 fr(il, i) = bsum(il) / csum(il)
396 endif
397 enddo
398 enddo
399
400 ! reset counter and return
401
402 do il = 1, ncum
403 sig(il, klev) = 2.0
404 enddo
405
406 do i = 1, klev
407 do il = 1, ncum
408 upwd(il, i) = 0.0
409 dnwd(il, i) = 0.0
410 enddo
411 enddo
412
413 do i = 1, nl
414 do il = 1, ncum
415 dnwd0(il, i) = - mp(il, i)
416 enddo
417 enddo
418 do i = nl + 1, klev
419 do il = 1, ncum
420 dnwd0(il, i) = 0.
421 enddo
422 enddo
423
424 do i = 1, nl
425 do il = 1, ncum
426 if (i >= icb(il) .and. i <= inb(il)) then
427 upwd(il, i) = 0.0
428 dnwd(il, i) = 0.0
429 endif
430 enddo
431 enddo
432
433 do i = 1, nl
434 do k = 1, nl
435 do il = 1, ncum
436 up1(il, k, i) = 0.0
437 dn1(il, k, i) = 0.0
438 enddo
439 enddo
440 enddo
441
442 do i = 1, nl
443 do k = i, nl
444 do n = 1, i - 1
445 do il = 1, ncum
446 if (i >= icb(il).and.i <= inb(il).and.k <= inb(il)) then
447 up1(il, k, i) = up1(il, k, i) + ment(il, n, k)
448 dn1(il, k, i) = dn1(il, k, i) - ment(il, k, n)
449 endif
450 enddo
451 enddo
452 enddo
453 enddo
454
455 do i = 2, nl
456 do k = i, nl
457 do il = 1, ncum
458 if (i <= inb(il).and.k <= inb(il)) then
459 upwd(il, i) = upwd(il, i) + m(il, k) + up1(il, k, i)
460 dnwd(il, i) = dnwd(il, i) + dn1(il, k, i)
461 endif
462 enddo
463 enddo
464 enddo
465
466 ! D\'etermination de la variation de flux ascendant entre
467 ! deux niveaux non dilu\'es mike
468
469 do i = 1, nl
470 do il = 1, ncum
471 mike(il, i) = m(il, i)
472 enddo
473 enddo
474
475 do i = nl + 1, klev
476 do il = 1, ncum
477 mike(il, i) = 0.
478 enddo
479 enddo
480
481 do i = 1, klev
482 do il = 1, ncum
483 ma(il, i) = 0
484 enddo
485 enddo
486
487 do i = 1, nl
488 do j = i, nl
489 do il = 1, ncum
490 ma(il, i) = ma(il, i) + m(il, j)
491 enddo
492 enddo
493 enddo
494
495 do i = nl + 1, klev
496 do il = 1, ncum
497 ma(il, i) = 0.
498 enddo
499 enddo
500
501 do i = 1, nl
502 do il = 1, ncum
503 if (i <= (icb(il) - 1)) then
504 ma(il, i) = 0
505 endif
506 enddo
507 enddo
508
509 ! icb repr\'esente le niveau o\`u se trouve la base du nuage, et
510 ! inb le sommet du nuage
511
512 do i = 1, klev
513 DO il = 1, ncum
514 rdcp = (rrd * (1. - rr(il, i)) - rr(il, i) * rrv) &
515 / (cpd * (1. - rr(il, i)) + rr(il, i) * cpv)
516 tls(il, i) = t(il, i) * (1000.0 / p(il, i))**rdcp
517 tps(il, i) = tp(il, i)
518 end DO
519 enddo
520
521 ! Diagnose the in-cloud mixing ratio of condensed water
522
523 do i = 1, klev
524 do il = 1, ncum
525 mac(il, i) = 0.0
526 wa(il, i) = 0.0
527 siga(il, i) = 0.0
528 sax(il, i) = 0.0
529 enddo
530 enddo
531
532 do i = minorig, nl
533 do k = i + 1, nl + 1
534 do il = 1, ncum
535 if (i <= inb(il) .and. k <= (inb(il) + 1)) then
536 mac(il, i) = mac(il, i) + m(il, k)
537 endif
538 enddo
539 enddo
540 enddo
541
542 do i = 1, nl
543 do j = 1, i
544 do il = 1, ncum
545 if (i >= icb(il) .and. i <= (inb(il) - 1) &
546 .and. j >= icb(il)) then
547 sax(il, i) = sax(il, i) + rrd * (tvp(il, j) - tv(il, j)) &
548 * (ph(il, j) - ph(il, j + 1)) / p(il, j)
549 endif
550 enddo
551 enddo
552 enddo
553
554 do i = 1, nl
555 do il = 1, ncum
556 if (i >= icb(il) .and. i <= (inb(il) - 1) &
557 .and. sax(il, i) > 0.0) then
558 wa(il, i) = sqrt(2. * sax(il, i))
559 endif
560 enddo
561 enddo
562
563 do i = 1, nl
564 do il = 1, ncum
565 if (wa(il, i) > 0.0) siga(il, i) = mac(il, i) / wa(il, i) * rrd &
566 * tvp(il, i) / p(il, i) / 100. / delta
567 siga(il, i) = min(siga(il, i), 1.0)
568
569 if (iflag_clw == 0) then
570 qcondc(il, i) = siga(il, i) * clw(il, i) * (1. - ep(il, i)) &
571 + (1. - siga(il, i)) * qcond(il, i)
572 else if (iflag_clw == 1) then
573 qcondc(il, i) = qcond(il, i)
574 endif
575 enddo
576 enddo
577
578 end SUBROUTINE cv30_yield
579
580 end module cv30_yield_m

  ViewVC Help
Powered by ViewVC 1.1.21