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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 188 - (hide annotations)
Tue Mar 22 16:31:39 2016 UTC (8 years, 2 months ago) by guez
File size: 19501 byte(s)
Removed argument ncum of cv30_unsat, arguments nloc, ncum, nd, na of cv30_yield.

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

  ViewVC Help
Powered by ViewVC 1.1.21