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