/[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 200 - (hide annotations)
Thu Jun 2 15:40:30 2016 UTC (7 years, 11 months ago) by guez
File size: 19036 byte(s)
Changes results.
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 195 ! Tendencies, precipitation, variables of interface with other
13     ! processes, etc.
14    
15 guez 187 use conema3_m, only: iflag_clw
16 guez 195 use cv30_param_m, only: minorig, nl, sigd
17 guez 198 use cv_thermo_m, only: cl, cpd, cpv, rowl, rrd, rrv
18 guez 188 USE dimphy, ONLY: klev, klon
19 guez 198 use SUPHEC_M, only: rg
20 guez 187
21 guez 97 ! inputs:
22 guez 188 integer, intent(in):: icb(:), inb(:) ! (ncum)
23 guez 97 real, intent(in):: delt
24 guez 200 real, intent(in):: t(klon, klev), rr(klon, klev)
25     real, intent(in):: u(klon, klev), v(klon, klev)
26 guez 188 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 guez 195 real, intent(in):: vp(:, 2:) ! (ncum, 2:nl)
36     real, intent(in):: wt(:, :) ! (ncum, nl - 1)
37 guez 189 real, intent(in):: water(:, :), evap(:, :) ! (ncum, nl)
38     real, intent(in):: b(:, :) ! (ncum, nl - 1)
39 guez 188 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 guez 47
46 guez 200 ! outputs:
47 guez 196 integer, intent(out):: iflag(:) ! (ncum)
48 guez 188 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 guez 97
58 guez 188 ! Local:
59 guez 195 real, parameter:: delta = 0.01 ! interface cloud parameterization
60 guez 188 integer ncum
61 guez 200 integer i, k, il, n, j
62     real awat, delti
63 guez 105 real ax, bx, cx, dx
64 guez 97 real cpinv, rdcp, dpinv
65 guez 188 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 guez 47
72 guez 97 !-------------------------------------------------------------
73 guez 47
74 guez 188 ncum = size(icb)
75 guez 196 iflag = 0
76 guez 188
77 guez 97 ! initialization:
78 guez 47
79 guez 188 delti = 1.0 / delt
80 guez 47
81 guez 188 do il = 1, ncum
82     precip(il) = 0.0
83     VPrecip(il, klev + 1) = 0.
84 guez 97 enddo
85 guez 47
86 guez 188 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 guez 47 enddo
97 guez 97 enddo
98 guez 47
99 guez 188 do i = 1, nl
100     do il = 1, ncum
101     lvcp(il, i) = lv(il, i) / cpn(il, i)
102 guez 47 enddo
103 guez 97 enddo
104 guez 47
105 guez 200 ! calculate surface precipitation in mm / day
106 guez 47
107 guez 188 do il = 1, ncum
108     if (ep(il, inb(il)) >= 1e-4) precip(il) = wt(il, 1) * sigd &
109 guez 198 * water(il, 1) * 86400. * 1000. / (rowl * rg)
110 guez 97 enddo
111 guez 47
112 guez 188 ! CALCULATE VERTICAL PROFILE OF PRECIPITATIONs IN kg / m2 / s ===
113 guez 187
114 guez 97 ! MAF rajout pour lessivage
115 guez 188 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 guez 198 / rg
119 guez 47 end do
120 guez 97 end do
121 guez 187
122 guez 200 ! calculate tendencies of lowest level potential temperature
123     ! and mixing ratio
124 guez 187
125 guez 188 do il = 1, ncum
126     work(il) = 1.0 / (ph(il, 1) - ph(il, 2))
127     am(il) = 0.0
128 guez 97 enddo
129 guez 47
130 guez 188 do k = 2, nl
131     do il = 1, ncum
132     if (k <= inb(il)) am(il) = am(il) + m(il, k)
133 guez 47 enddo
134 guez 97 enddo
135 guez 47
136 guez 188 do il = 1, ncum
137 guez 198 if (0.01 * rg * work(il) * am(il) >= delti) iflag(il) = 1
138 guez 47
139 guez 198 ft(il, 1) = 0.01 * rg * work(il) * am(il) * (t(il, 2) - t(il, 1) &
140 guez 200 + (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 guez 47
146 guez 187 !jyg1 Correction pour mieux conserver l'eau (conformite avec CONVECT4.3)
147 guez 198 ! (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 guez 188 * work(il) + sigd * 0.5 * (evap(il, 1) + evap(il, 2))
150     ! + tard : + sigd * evap(il, 1)
151 guez 47
152 guez 198 fr(il, 1) = fr(il, 1) + 0.01 * rg * am(il) * (rr(il, 2) - rr(il, 1)) &
153 guez 188 * work(il)
154 guez 47
155 guez 198 fu(il, 1) = fu(il, 1) + 0.01 * rg * work(il) * (mp(il, 2) &
156 guez 188 * (up(il, 2) - u(il, 1)) + am(il) * (u(il, 2) - u(il, 1)))
157 guez 198 fv(il, 1) = fv(il, 1) + 0.01 * rg * work(il) * (mp(il, 2) &
158 guez 188 * (vp(il, 2) - v(il, 1)) + am(il) * (v(il, 2) - v(il, 1)))
159 guez 200 enddo
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 198 fr(il, 1) = fr(il, 1) + 0.01 * rg * work(il) * ment(il, j, 1) &
165 guez 188 * (qent(il, j, 1) - rr(il, 1))
166 guez 198 fu(il, 1) = fu(il, 1) + 0.01 * rg * work(il) * ment(il, j, 1) &
167 guez 188 * (uent(il, j, 1) - u(il, 1))
168 guez 198 fv(il, 1) = fv(il, 1) + 0.01 * rg * work(il) * ment(il, j, 1) &
169 guez 188 * (vent(il, j, 1) - v(il, 1))
170     endif
171 guez 47 enddo
172 guez 97 enddo
173 guez 47
174 guez 200 ! calculate tendencies of potential temperature and mixing ratio
175     ! at levels above the lowest level
176 guez 47
177 guez 200 ! 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 guez 200 if (any(inb >= i)) then
182 guez 188 amp1(:ncum) = 0.
183     ad(:ncum) = 0.
184 guez 47
185 guez 188 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 guez 97 end do
192 guez 47
193 guez 188 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 guez 97 end do
201     end do
202 guez 47
203 guez 188 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 guez 97 end do
211     end do
212 guez 47
213 guez 188 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 guez 47
218 guez 198 if (0.01 * rg * dpinv * amp1(il) >= delti) iflag(il) = 1
219 guez 47
220 guez 198 ft(il, i) = 0.01 * rg * dpinv * (amp1(il) * (t(il, i + 1) &
221 guez 188 - 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 guez 200 * (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 guez 188 * (hp(il, i) - h(il, i) + t(il, i) * (cpv - cpd) &
229 guez 200 * (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 guez 198 fr(il, i) = 0.01 * rg * dpinv * (amp1(il) * (rr(il, i + 1) &
233 guez 188 - rr(il, i)) - ad(il) * (rr(il, i) - rr(il, i - 1)))
234 guez 198 fu(il, i) = fu(il, i) + 0.01 * rg * dpinv * (amp1(il) &
235 guez 188 * (u(il, i + 1) - u(il, i)) - ad(il) * (u(il, i) &
236     - u(il, i - 1)))
237 guez 198 fv(il, i) = fv(il, i) + 0.01 * rg * dpinv * (amp1(il) &
238 guez 188 * (v(il, i + 1) - v(il, i)) - ad(il) * (v(il, i) &
239     - v(il, i - 1)))
240     endif
241     end do
242 guez 47
243 guez 188 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 guez 47
249 guez 188 awat = elij(il, k, i) - (1. - ep(il, i)) * clw(il, i)
250     awat = amax1(awat, 0.0)
251 guez 47
252 guez 198 fr(il, i) = fr(il, i) + 0.01 * rg * dpinv &
253 guez 188 * ment(il, k, i) * (qent(il, k, i) - awat - rr(il, i))
254 guez 198 fu(il, i) = fu(il, i) + 0.01 * rg * dpinv &
255 guez 188 * ment(il, k, i) * (uent(il, k, i) - u(il, i))
256 guez 198 fv(il, i) = fv(il, i) + 0.01 * rg * dpinv &
257 guez 188 * ment(il, k, i) * (vent(il, k, i) - v(il, i))
258 guez 47
259 guez 188 ! (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 guez 97 end do
265 guez 47
266 guez 188 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 guez 47
272 guez 198 fr(il, i) = fr(il, i) + 0.01 * rg * dpinv &
273 guez 188 * ment(il, k, i) * (qent(il, k, i) - rr(il, i))
274 guez 198 fu(il, i) = fu(il, i) + 0.01 * rg * dpinv &
275 guez 188 * ment(il, k, i) * (uent(il, k, i) - u(il, i))
276 guez 198 fv(il, i) = fv(il, i) + 0.01 * rg * dpinv &
277 guez 188 * ment(il, k, i) * (vent(il, k, i) - v(il, i))
278     endif
279     end do
280 guez 97 end do
281 guez 47
282 guez 188 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 guez 47
287 guez 188 ! 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 guez 198 + evap(il, i + 1)) + 0.01 * rg * (mp(il, i + 1) &
291 guez 188 * (rp(il, i + 1) - rr(il, i)) - mp(il, i) * (rp(il, i) &
292     - rr(il, i - 1))) * dpinv
293 guez 47
294 guez 198 fu(il, i) = fu(il, i) + 0.01 * rg * (mp(il, i + 1) &
295 guez 188 * (up(il, i + 1) - u(il, i)) - mp(il, i) * (up(il, i) &
296     - u(il, i - 1))) * dpinv
297 guez 198 fv(il, i) = fv(il, i) + 0.01 * rg * (mp(il, i + 1) &
298 guez 188 * (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 guez 47
303 guez 188 ! sb: interface with the cloud parameterization:
304 guez 47
305 guez 188 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 guez 47
315 guez 188 ! (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 guez 47
323 guez 188 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 guez 47
331 guez 200 ! 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 guez 47
335 guez 188 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 guez 47
345 guez 188 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 guez 47
352 guez 188 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 guez 47
359 guez 188 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 guez 47
366 guez 97 end do
367 guez 47
368 guez 200 ! homoginize tendencies below cloud base
369 guez 187
370 guez 188 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 guez 97 enddo
376 guez 47
377 guez 188 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 guez 97 endif
388 guez 47 enddo
389 guez 97 enddo
390 guez 47
391 guez 188 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 guez 97 endif
397 guez 47 enddo
398 guez 97 enddo
399 guez 47
400 guez 200 ! reset counter and return
401 guez 187
402 guez 188 do il = 1, ncum
403     sig(il, klev) = 2.0
404 guez 97 enddo
405 guez 47
406 guez 188 do i = 1, klev
407     do il = 1, ncum
408     upwd(il, i) = 0.0
409     dnwd(il, i) = 0.0
410 guez 47 enddo
411 guez 97 enddo
412 guez 47
413 guez 188 do i = 1, nl
414     do il = 1, ncum
415     dnwd0(il, i) = - mp(il, i)
416 guez 47 enddo
417 guez 97 enddo
418 guez 188 do i = nl + 1, klev
419     do il = 1, ncum
420     dnwd0(il, i) = 0.
421 guez 47 enddo
422 guez 97 enddo
423 guez 47
424 guez 188 do i = 1, nl
425     do il = 1, ncum
426 guez 187 if (i >= icb(il) .and. i <= inb(il)) then
427 guez 188 upwd(il, i) = 0.0
428     dnwd(il, i) = 0.0
429 guez 97 endif
430 guez 47 enddo
431 guez 97 enddo
432 guez 47
433 guez 188 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 guez 97 enddo
439 guez 47 enddo
440 guez 97 enddo
441 guez 47
442 guez 188 do i = 1, nl
443     do k = i, nl
444     do n = 1, i - 1
445     do il = 1, ncum
446 guez 187 if (i >= icb(il).and.i <= inb(il).and.k <= inb(il)) then
447 guez 188 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 guez 97 endif
450     enddo
451     enddo
452 guez 47 enddo
453 guez 97 enddo
454 guez 47
455 guez 188 do i = 2, nl
456     do k = i, nl
457     do il = 1, ncum
458 guez 187 if (i <= inb(il).and.k <= inb(il)) then
459 guez 188 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 guez 97 endif
462     enddo
463 guez 47 enddo
464 guez 97 enddo
465 guez 47
466 guez 200 ! D\'etermination de la variation de flux ascendant entre
467     ! deux niveaux non dilu\'es mike
468 guez 47
469 guez 188 do i = 1, nl
470     do il = 1, ncum
471     mike(il, i) = m(il, i)
472 guez 47 enddo
473 guez 97 enddo
474 guez 47
475 guez 188 do i = nl + 1, klev
476     do il = 1, ncum
477     mike(il, i) = 0.
478 guez 47 enddo
479 guez 97 enddo
480 guez 47
481 guez 188 do i = 1, klev
482     do il = 1, ncum
483     ma(il, i) = 0
484 guez 47 enddo
485 guez 97 enddo
486 guez 47
487 guez 188 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 guez 97 enddo
492 guez 47 enddo
493 guez 97 enddo
494 guez 47
495 guez 188 do i = nl + 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 il = 1, ncum
503     if (i <= (icb(il) - 1)) then
504     ma(il, i) = 0
505 guez 97 endif
506 guez 47 enddo
507 guez 97 enddo
508 guez 47
509 guez 200 ! icb repr\'esente le niveau o\`u se trouve la base du nuage, et
510     ! inb le sommet du nuage
511 guez 47
512 guez 188 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 guez 97 end DO
519     enddo
520 guez 47
521 guez 200 ! Diagnose the in-cloud mixing ratio of condensed water
522 guez 47
523 guez 188 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 guez 47
532 guez 188 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 guez 47
542 guez 188 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 guez 97
554 guez 188 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 guez 47
563 guez 188 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 guez 187 if (iflag_clw == 0) then
570 guez 188 qcondc(il, i) = siga(il, i) * clw(il, i) * (1. - ep(il, i)) &
571     + (1. - siga(il, i)) * qcond(il, i)
572 guez 187 else if (iflag_clw == 1) then
573 guez 188 qcondc(il, i) = qcond(il, i)
574 guez 97 endif
575 guez 188 enddo
576     enddo
577 guez 47
578 guez 185 end SUBROUTINE cv30_yield
579 guez 97
580 guez 185 end module cv30_yield_m

  ViewVC Help
Powered by ViewVC 1.1.21