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

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

Parent Directory Parent Directory | Revision Log Revision Log

Revision 186 - (hide annotations)
Mon Mar 21 15:36:26 2016 UTC (8 years, 1 month ago) by guez
File size: 13265 byte(s) 
Removed variables nlm and nlp of module cv30_param_m. We do not
believe much in the benefit of these intermediary variables so we go
for clarity.

Removed variable noff of module cv30_param_m. Never used anywhere
else. Just set the value of nl explicitly in cv30_param.

Removed argument nd of cv30_param. Only called with nd = klev.

Replaced calls to zilch by array assignments. There was a strange
double call to zilch with the same arguments in cv30_mixing.

Removed procedure cv_flag. Just set the value of variable cvflag_grav
of module cvflag at declaration.
1 guez 185 module cv30_unsat_m
2 guez 47
3 guez 97 implicit none
4 guez 47
5 guez 97 contains
6 guez 47
7 guez 186 SUBROUTINE cv30_unsat(nloc, ncum, nd, na, icb, inb, t, rr, rs, gz, u, v, p, &
8     ph, th, tv, lv, cpn, ep, sigp, clw, m, ment, elij, delt, plcl, mp, rp, &
9     up, vp, wt, water, evap, b)
10 guez 47
11 guez 186 use cv30_param_m, only: nl, sigd
12     use cvflag, only: cvflag_grav
13     use cvthermo, only: cpd, ginv, grav
14 guez 47
15 guez 97 ! inputs:
16 guez 186 integer, intent(in):: nloc, ncum, nd, na
17     integer, intent(in):: icb(:), inb(:) ! (ncum)
18     real t(nloc, nd), rr(nloc, nd), rs(nloc, nd)
19     real gz(nloc, na)
20     real u(nloc, nd), v(nloc, nd)
21     real p(nloc, nd), ph(nloc, nd + 1)
22     real th(nloc, na)
23     real tv(nloc, na)
24     real lv(nloc, na)
25     real cpn(nloc, na)
26     real ep(nloc, na), sigp(nloc, na), clw(nloc, na)
27     real m(nloc, na), ment(nloc, na, na), elij(nloc, na, na)
28 guez 97 real, intent(in):: delt
29     real plcl(nloc)
30 guez 47
31 guez 97 ! outputs:
32 guez 186 real mp(nloc, na), rp(nloc, na), up(nloc, na), vp(nloc, na)
33     real wt(nloc, na), water(nloc, na), evap(nloc, na)
34     real b(:, :) ! (nloc, na)
35 guez 47
36 guez 186 ! Local:
37     integer i, j, il, num1
38 guez 97 real tinv, delti
39     real awat, afac, afac1, afac2, bfac
40     real pr1, pr2, sigt, b6, c6, revap, tevap, delth
41     real amfac, amp2, xf, tf, fac2, ur, sru, fac, d, af, bf
42     real ampmax
43 guez 186 real lvcp(nloc, na)
44 guez 97 real wdtrain(nloc)
45     logical lwork(nloc)
46 guez 47
47 guez 97 !------------------------------------------------------
48 guez 47
49 guez 186 delti = 1. / delt
50     tinv = 1. / 3.
51     mp = 0.
52 guez 47
53 guez 186 do i = 1, nl
54     do il = 1, ncum
55     mp(il, i) = 0.
56     rp(il, i) = rr(il, i)
57     up(il, i) = u(il, i)
58     vp(il, i) = v(il, i)
59     wt(il, i) = 0.001
60     water(il, i) = 0.
61     evap(il, i) = 0.
62     b(il, i) = 0.
63     lvcp(il, i) = lv(il, i) / cpn(il, i)
64 guez 97 enddo
65     enddo
66 guez 47
67 guez 186 ! check whether ep(inb) = 0, if so, skip precipitating
68     ! downdraft calculation
69 guez 47
70 guez 186 do il = 1, ncum
71     lwork(il) = .TRUE.
72     if (ep(il, inb(il)) < 0.0001) lwork(il) = .FALSE.
73 guez 97 enddo
74 guez 47
75 guez 186 wdtrain(:ncum) = 0.
76 guez 47
77 guez 186 downdraft_loop: DO i = nl + 1, 1, - 1
78     num1 = 0
79 guez 47
80 guez 186 do il = 1, ncum
81     if (i <= inb(il) .and. lwork(il)) num1 = num1 + 1
82 guez 97 enddo
83 guez 47
84 guez 186 if (num1 > 0) then
85     ! integrate liquid water equation to find condensed water
86     ! and condensed water flux
87 guez 47
88 guez 186 ! calculate detrained precipitation
89 guez 47
90 guez 186 do il = 1, ncum
91     if (i <= inb(il) .and. lwork(il)) then
92     if (cvflag_grav) then
93     wdtrain(il) = grav * ep(il, i) * m(il, i) * clw(il, i)
94     else
95     wdtrain(il) = 10. * ep(il, i) * m(il, i) * clw(il, i)
96     endif
97 guez 97 endif
98 guez 186 enddo
99 guez 47
100 guez 186 if (i > 1) then
101     do j = 1, i - 1
102     do il = 1, ncum
103     if (i <= inb(il) .and. lwork(il)) then
104     awat = elij(il, j, i) - (1. - ep(il, i)) * clw(il, i)
105     awat = amax1(awat, 0.)
106     if (cvflag_grav) then
107     wdtrain(il) = wdtrain(il) + grav * awat &
108     * ment(il, j, i)
109     else
110     wdtrain(il) = wdtrain(il) + 10. * awat * ment(il, j, i)
111     endif
112 guez 97 endif
113 guez 186 enddo
114     end do
115     endif
116 guez 47
117 guez 186 ! find rain water and evaporation using provisional
118     ! estimates of rp(i)and rp(i - 1)
119 guez 47
120 guez 186 do il = 1, ncum
121     if (i <= inb(il) .and. lwork(il)) then
122     wt(il, i) = 45.
123 guez 47
124 guez 186 if (i < inb(il)) then
125     rp(il, i) = rp(il, i + 1) + (cpd * (t(il, i + 1) &
126     - t(il, i)) + gz(il, i + 1) - gz(il, i)) / lv(il, i)
127     rp(il, i) = 0.5 * (rp(il, i) + rr(il, i))
128     endif
129     rp(il, i) = amax1(rp(il, i), 0.)
130     rp(il, i) = amin1(rp(il, i), rs(il, i))
131     rp(il, inb(il)) = rr(il, inb(il))
132 guez 47
133 guez 186 if (i == 1) then
134     afac = p(il, 1) * (rs(il, 1) - rp(il, 1)) &
135     / (1e4 + 2000. * p(il, 1) * rs(il, 1))
136     else
137     rp(il, i - 1) = rp(il, i) + (cpd * (t(il, i) &
138     - t(il, i - 1)) + gz(il, i) - gz(il, i - 1)) / lv(il, i)
139     rp(il, i - 1) = 0.5 * (rp(il, i - 1) + rr(il, i - 1))
140     rp(il, i - 1) = amin1(rp(il, i - 1), rs(il, i - 1))
141     rp(il, i - 1) = amax1(rp(il, i - 1), 0.)
142     afac1 = p(il, i) * (rs(il, i) - rp(il, i)) &
143     / (1e4 + 2000. * p(il, i) * rs(il, i))
144     afac2 = p(il, i - 1) * (rs(il, i - 1) - rp(il, i - 1)) &
145     / (1e4 + 2000. * p(il, i - 1) * rs(il, i - 1))
146     afac = 0.5 * (afac1 + afac2)
147     endif
148     if (i == inb(il))afac = 0.
149     afac = amax1(afac, 0.)
150     bfac = 1. / (sigd * wt(il, i))
151 guez 47
152 guez 186 ! prise en compte de la variation progressive de sigt dans
153     ! les couches icb et icb - 1:
154     ! pour plcl < ph(i + 1), pr1 = 0 & pr2 = 1
155     ! pour plcl > ph(i), pr1 = 1 & pr2 = 0
156     ! pour ph(i + 1) < plcl < ph(i), pr1 est la proportion a cheval
157     ! sur le nuage, et pr2 est la proportion sous la base du
158     ! nuage.
159     pr1 = (plcl(il) - ph(il, i + 1)) / (ph(il, i) - ph(il, i + 1))
160     pr1 = max(0., min(1., pr1))
161     pr2 = (ph(il, i) - plcl(il)) / (ph(il, i) - ph(il, i + 1))
162     pr2 = max(0., min(1., pr2))
163     sigt = sigp(il, i) * pr1 + pr2
164 guez 47
165 guez 186 b6 = bfac * 50. * sigd * (ph(il, i) - ph(il, i + 1)) * sigt &
166     * afac
167     c6 = water(il, i + 1) + bfac * wdtrain(il) - 50. * sigd * bfac &
168     * (ph(il, i) - ph(il, i + 1)) * evap(il, i + 1)
169     if (c6 > 0.) then
170     revap = 0.5 * (- b6 + sqrt(b6 * b6 + 4. * c6))
171     evap(il, i) = sigt * afac * revap
172     water(il, i) = revap * revap
173 guez 97 else
174 guez 186 evap(il, i) = - evap(il, i + 1) + 0.02 * (wdtrain(il) &
175     + sigd * wt(il, i) * water(il, i + 1)) &
176     / (sigd * (ph(il, i) - ph(il, i + 1)))
177     end if
178 guez 47
179 guez 186 ! calculate precipitating downdraft mass flux under
180     ! hydrostatic approximation
181    
182     if (i /= 1) then
183     tevap = amax1(0., evap(il, i))
184     delth = amax1(0.001, (th(il, i) - th(il, i - 1)))
185 guez 97 if (cvflag_grav) then
186 guez 186 mp(il, i) = 100. * ginv * lvcp(il, i) * sigd * tevap &
187     * (p(il, i - 1) - p(il, i)) / delth
188 guez 97 else
189 guez 186 mp(il, i) = 10. * lvcp(il, i) * sigd * tevap &
190     * (p(il, i - 1) - p(il, i)) / delth
191 guez 97 endif
192 guez 47
193 guez 186 ! if hydrostatic assumption fails,
194     ! solve cubic difference equation for downdraft theta
195     ! and mass flux from two simultaneous differential eqns
196 guez 47
197 guez 186 amfac = sigd * sigd * 70. * ph(il, i) &
198     * (p(il, i - 1) - p(il, i)) &
199     * (th(il, i) - th(il, i - 1)) / (tv(il, i) * th(il, i))
200     amp2 = abs(mp(il, i + 1) * mp(il, i + 1) - mp(il, i) &
201     * mp(il, i))
202     if (amp2 > (0.1 * amfac)) then
203     xf = 100. * sigd * sigd * sigd * (ph(il, i) &
204     - ph(il, i + 1))
205     tf = b(il, i) - 5. * (th(il, i) - th(il, i - 1)) &
206     * t(il, i) / (lvcp(il, i) * sigd * th(il, i))
207     af = xf * tf + mp(il, i + 1) * mp(il, i + 1) * tinv
208     bf = 2. * (tinv * mp(il, i + 1))**3 + tinv &
209     * mp(il, i + 1) * xf * tf + 50. * (p(il, i - 1) &
210     - p(il, i)) * xf * tevap
211     fac2 = 1.
212     if (bf < 0.)fac2 = - 1.
213     bf = abs(bf)
214     ur = 0.25 * bf * bf - af * af * af * tinv * tinv * tinv
215     if (ur >= 0.) then
216     sru = sqrt(ur)
217     fac = 1.
218     if ((0.5 * bf - sru) < 0.)fac = - 1.
219     mp(il, i) = mp(il, i + 1) * tinv &
220     + (0.5 * bf + sru)**tinv &
221     + fac * (abs(0.5 * bf - sru))**tinv
222     else
223     d = atan(2. * sqrt(- ur) / (bf + 1e-28))
224     if (fac2 < 0.)d = 3.14159 - d
225     mp(il, i) = mp(il, i + 1) * tinv + 2. &
226     * sqrt(af * tinv) * cos(d * tinv)
227     endif
228     mp(il, i) = amax1(0., mp(il, i))
229 guez 47
230 guez 186 if (cvflag_grav) then
231     ! Il y a vraisemblablement une erreur dans la
232     ! ligne 2 suivante: il faut diviser par (mp(il,
233     ! i) * sigd * grav) et non par (mp(il, i) + sigd
234     ! * 0.1). Et il faut bien revoir les facteurs
235     ! 100.
236     b(il, i - 1) = b(il, i) + 100. * (p(il, i - 1) &
237     - p(il, i)) * tevap / (mp(il, i) + sigd * 0.1) &
238     - 10. * (th(il, i) - th(il, i - 1)) * t(il, i) &
239     / (lvcp(il, i) * sigd * th(il, i))
240     else
241     b(il, i - 1) = b(il, i) + 100. * (p(il, i - 1) &
242     - p(il, i)) * tevap / (mp(il, i) + sigd * 0.1) &
243     - 10. * (th(il, i) - th(il, i - 1)) * t(il, i) &
244     / (lvcp(il, i) * sigd * th(il, i))
245     endif
246     b(il, i - 1) = amax1(b(il, i - 1), 0.)
247     endif
248 guez 47
249 guez 186 ! limit magnitude of mp(i) to meet cfl condition
250 guez 47
251 guez 186 ampmax = 2. * (ph(il, i) - ph(il, i + 1)) * delti
252     amp2 = 2. * (ph(il, i - 1) - ph(il, i)) * delti
253     ampmax = amin1(ampmax, amp2)
254     mp(il, i) = amin1(mp(il, i), ampmax)
255    
256     ! force mp to decrease linearly to zero
257     ! between cloud base and the surface
258    
259     if (p(il, i) > p(il, icb(il))) then
260     mp(il, i) = mp(il, icb(il)) * (p(il, 1) - p(il, i)) &
261     / (p(il, 1) - p(il, icb(il)))
262 guez 97 endif
263 guez 186 endif ! i == 1
264 guez 47
265 guez 186 ! find mixing ratio of precipitating downdraft
266 guez 47
267 guez 186 if (i /= inb(il)) then
268     rp(il, i) = rr(il, i)
269    
270     if (mp(il, i) > mp(il, i + 1)) then
271 guez 97 if (cvflag_grav) then
272 guez 186 rp(il, i) = rp(il, i + 1) * mp(il, i + 1) + rr(il, i) &
273     * (mp(il, i) - mp(il, i + 1)) + 100. * ginv &
274     * 0.5 * sigd * (ph(il, i) - ph(il, i + 1)) &
275     * (evap(il, i + 1) + evap(il, i))
276 guez 97 else
277 guez 186 rp(il, i) = rp(il, i + 1) * mp(il, i + 1) + rr(il, i) &
278     * (mp(il, i) - mp(il, i + 1)) + 5. * sigd &
279     * (ph(il, i) - ph(il, i + 1)) &
280     * (evap(il, i + 1) + evap(il, i))
281 guez 97 endif
282 guez 186 rp(il, i) = rp(il, i) / mp(il, i)
283     up(il, i) = up(il, i + 1) * mp(il, i + 1) + u(il, i) &
284     * (mp(il, i) - mp(il, i + 1))
285     up(il, i) = up(il, i) / mp(il, i)
286     vp(il, i) = vp(il, i + 1) * mp(il, i + 1) + v(il, i) &
287     * (mp(il, i) - mp(il, i + 1))
288     vp(il, i) = vp(il, i) / mp(il, i)
289     else
290     if (mp(il, i + 1) > 1e-16) then
291     if (cvflag_grav) then
292     rp(il, i) = rp(il, i + 1) &
293     + 100. * ginv * 0.5 * sigd * (ph(il, i) &
294     - ph(il, i + 1)) &
295     * (evap(il, i + 1) + evap(il, i)) &
296     / mp(il, i + 1)
297     else
298     rp(il, i) = rp(il, i + 1) &
299     + 5. * sigd * (ph(il, i) - ph(il, i + 1)) &
300     * (evap(il, i + 1) + evap(il, i)) &
301     / mp(il, i + 1)
302     endif
303     up(il, i) = up(il, i + 1)
304     vp(il, i) = vp(il, i + 1)
305     endif
306 guez 97 endif
307 guez 186 rp(il, i) = amin1(rp(il, i), rs(il, i))
308     rp(il, i) = amax1(rp(il, i), 0.)
309 guez 97 endif
310     endif
311 guez 186 end do
312     end if
313     end DO downdraft_loop
314 guez 47
315 guez 185 end SUBROUTINE cv30_unsat
316 guez 97
317 guez 185 end module cv30_unsat_m
  ViewVC Help
Powered by ViewVC 1.1.21