/[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 190 - (hide annotations)
Thu Apr 14 15:15:56 2016 UTC (8 years, 1 month ago) by guez
File size: 19487 byte(s) 
Created module cv_thermo_m around procedure cv_thermo. Moved variables
from module cvthermo to module cv_thermo_m, where they are defined.

In ini_histins and initphysto, using part of rlon and rlat from
phyetat0_m is pretending that we do not know about the dynamical grid,
while the way we extract zx_lon(:, 1) and zx_lat(1, :) depends on
ordering inside rlon and rlat. So we might as well simplify and
clarify by using directly rlonv and rlatu.

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