10 |
! Unsaturated (precipitating) downdrafts |
! Unsaturated (precipitating) downdrafts |
11 |
|
|
12 |
use cv30_param_m, only: nl, sigd |
use cv30_param_m, only: nl, sigd |
13 |
use cv_thermo_m, only: cpd, ginv |
use cv_thermo_m, only: ginv |
14 |
use SUPHEC_M, only: rg |
use SUPHEC_M, only: rg, rcpd |
15 |
|
|
16 |
integer, intent(in):: icb(:) ! (ncum) |
integer, intent(in):: icb(:) ! (ncum) |
17 |
! {2 <= icb <= nl - 3} |
! {2 <= icb <= nl - 3} |
20 |
! first model level above the level of neutral buoyancy of the |
! first model level above the level of neutral buoyancy of the |
21 |
! parcel (1 <= inb <= nl - 1) |
! parcel (1 <= inb <= nl - 1) |
22 |
|
|
23 |
real, intent(in):: t(:, :), q(:, :), qs(:, :) ! (ncum, nl) |
real, intent(in):: t(:, :) ! (ncum, nl) temperature (K) |
24 |
|
real, intent(in):: q(:, :), qs(:, :) ! (ncum, nl) |
25 |
real, intent(in):: gz(:, :) ! (klon, klev) |
real, intent(in):: gz(:, :) ! (klon, klev) |
26 |
real, intent(in):: u(:, :), v(:, :) ! (ncum, nl) |
real, intent(in):: u(:, :), v(:, :) ! (ncum, nl) |
27 |
real, intent(in):: p(:, :) ! (klon, klev) pressure at full level, in hPa |
real, intent(in):: p(:, :) ! (klon, klev) pressure at full level, in hPa |
28 |
real, intent(in):: ph(:, :) ! (ncum, klev + 1) |
real, intent(in):: ph(:, :) ! (ncum, klev + 1) |
29 |
real, intent(in):: th(:, :) ! (ncum, nl - 1) |
real, intent(in):: th(:, :) ! (ncum, nl - 1) potential temperature, in K |
30 |
real, intent(in):: tv(:, :) ! (klon, klev) |
real, intent(in):: tv(:, :) ! (klon, klev) |
31 |
|
|
32 |
real, intent(in):: lv(:, :) ! (ncum, nl) |
real, intent(in):: lv(:, :) ! (ncum, nl) |
33 |
real, intent(in):: cpn(:, :) ! (klon, klev) |
! specific latent heat of vaporization of water, in J kg-1 |
34 |
|
|
35 |
|
real, intent(in):: cpn(:, :) ! (ncum, nl) |
36 |
|
! specific heat capacity at constant pressure of humid air, in J K-1 kg-1 |
37 |
|
|
38 |
real, intent(in):: ep(:, :) ! (ncum, klev) |
real, intent(in):: ep(:, :) ! (ncum, klev) |
39 |
real, intent(in):: clw(:, :) ! (ncum, klev) |
real, intent(in):: clw(:, :) ! (ncum, klev) |
40 |
real, intent(in):: m(:, :) ! (ncum, klev) |
real, intent(in):: m(:, :) ! (ncum, klev) |
43 |
real, intent(in):: delt |
real, intent(in):: delt |
44 |
real, intent(in):: plcl(:) ! (ncum) |
real, intent(in):: plcl(:) ! (ncum) |
45 |
|
|
46 |
real, intent(out):: mp(:, :) ! (klon, klev) |
real, intent(out):: mp(:, :) |
47 |
! mass flux of the unsaturated downdraft, defined positive downward |
! (ncum, nl) Mass flux of the unsaturated downdraft, defined |
48 |
! M_p in Emanuel (1991 928) |
! positive downward, in kg m-2 s-1. M_p in Emanuel (1991 928). |
49 |
|
|
50 |
real, intent(out):: qp(:, :), up(:, :), vp(:, :) ! (ncum, nl) |
real, intent(out):: qp(:, :), up(:, :), vp(:, :) ! (ncum, nl) |
51 |
real, intent(out):: wt(:, :) ! (ncum, nl) |
real, intent(out):: wt(:, :) ! (ncum, nl) |
67 |
|
|
68 |
integer ncum |
integer ncum |
69 |
integer i, il, imax |
integer i, il, imax |
70 |
real tinv, delti |
real, parameter:: tinv = 1. / 3. |
71 |
|
real delti |
72 |
real afac, afac1, afac2, bfac |
real afac, afac1, afac2, bfac |
73 |
real pr1, sigt, b6, c6, revap, tevap, delth |
real pr1, sigt, b6, c6, revap, tevap, delth |
74 |
real amfac, amp2, xf, tf, fac2, ur, sru, fac, d, af, bf |
real xf, tf, fac2, ur, sru, fac, d, af, bf |
75 |
real ampmax |
real ampmax |
76 |
real lvcp(size(icb), nl) ! (ncum, nl) |
real lvcp(size(icb), nl) ! (ncum, nl) L_v / C_p, in K |
77 |
real wdtrain(size(icb)) ! (ncum) |
real wdtrain(size(icb)) ! (ncum) |
78 |
logical lwork(size(icb)) ! (ncum) |
logical lwork(size(icb)) ! (ncum) |
79 |
|
|
81 |
|
|
82 |
ncum = size(icb) |
ncum = size(icb) |
83 |
delti = 1. / delt |
delti = 1. / delt |
|
tinv = 1. / 3. |
|
84 |
mp = 0. |
mp = 0. |
85 |
b = 0. |
b = 0. |
86 |
|
|
124 |
wt(il, i) = 45. |
wt(il, i) = 45. |
125 |
|
|
126 |
if (i < inb(il)) then |
if (i < inb(il)) then |
127 |
qp(il, i) = qp(il, i + 1) + (cpd * (t(il, i + 1) - t(il, i)) & |
qp(il, i) = qp(il, i + 1) + (rcpd * (t(il, i + 1) - t(il, i)) & |
128 |
+ gz(il, i + 1) - gz(il, i)) / lv(il, i) |
+ gz(il, i + 1) - gz(il, i)) / lv(il, i) |
129 |
qp(il, i) = 0.5 * (qp(il, i) + q(il, i)) |
qp(il, i) = 0.5 * (qp(il, i) + q(il, i)) |
130 |
endif |
endif |
137 |
afac = p(il, 1) * (qs(il, 1) - qp(il, 1)) & |
afac = p(il, 1) * (qs(il, 1) - qp(il, 1)) & |
138 |
/ (1e4 + 2000. * p(il, 1) * qs(il, 1)) |
/ (1e4 + 2000. * p(il, 1) * qs(il, 1)) |
139 |
else |
else |
140 |
qp(il, i - 1) = qp(il, i) + (cpd * (t(il, i) - t(il, i - 1)) & |
qp(il, i - 1) = qp(il, i) + (rcpd * (t(il, i) - t(il, i - 1)) & |
141 |
+ gz(il, i) - gz(il, i - 1)) / lv(il, i) |
+ gz(il, i) - gz(il, i - 1)) / lv(il, i) |
142 |
qp(il, i - 1) = 0.5 * (qp(il, i - 1) + q(il, i - 1)) |
qp(il, i - 1) = 0.5 * (qp(il, i - 1) + q(il, i - 1)) |
143 |
qp(il, i - 1) = min(qp(il, i - 1), qs(il, i - 1)) |
qp(il, i - 1) = min(qp(il, i - 1), qs(il, i - 1)) |
194 |
! If hydrostatic assumption fails, solve cubic |
! If hydrostatic assumption fails, solve cubic |
195 |
! difference equation for downdraft theta and mass |
! difference equation for downdraft theta and mass |
196 |
! flux from two simultaneous differential equations |
! flux from two simultaneous differential equations |
197 |
|
if (abs(mp(il, i + 1)**2 - mp(il, i)**2) > 0.1 * sigd**2 & |
198 |
amfac = sigd**2 * 70. * ph(il, i) * (p(il, i - 1) - p(il, i)) & |
* 70. * ph(il, i) * (p(il, i - 1) - p(il, i)) & |
199 |
* (th(il, i) - th(il, i - 1)) / (tv(il, i) * th(il, i)) |
* (th(il, i) - th(il, i - 1)) / (tv(il, i) * th(il, i))) & |
200 |
amp2 = abs(mp(il, i + 1)**2 - mp(il, i)**2) |
then |
|
|
|
|
if (amp2 > 0.1 * amfac) then |
|
201 |
xf = 100. * sigd**3 * (ph(il, i) - ph(il, i + 1)) |
xf = 100. * sigd**3 * (ph(il, i) - ph(il, i + 1)) |
202 |
tf = b(il, i) - 5. * (th(il, i) - th(il, i - 1)) & |
tf = b(il, i) - 5. * (th(il, i) - th(il, i - 1)) & |
203 |
* t(il, i) / (lvcp(il, i) * sigd * th(il, i)) |
* t(il, i) / (lvcp(il, i) * sigd * th(il, i)) |
218 |
+ sru)**tinv + fac * (abs(0.5 * bf - sru))**tinv |
+ sru)**tinv + fac * (abs(0.5 * bf - sru))**tinv |
219 |
else |
else |
220 |
d = atan(2. * sqrt(- ur) / (bf + 1e-28)) |
d = atan(2. * sqrt(- ur) / (bf + 1e-28)) |
221 |
if (fac2 < 0.)d = 3.14159 - d |
if (fac2 < 0.) d = 3.14159 - d |
222 |
mp(il, i) = mp(il, i + 1) * tinv + 2. * sqrt(af * tinv) & |
mp(il, i) = mp(il, i + 1) * tinv + 2. * sqrt(af * tinv) & |
223 |
* cos(d * tinv) |
* cos(d * tinv) |
224 |
endif |
endif |
229 |
! suivante : il faut diviser par (mp(il, i) * sigd |
! suivante : il faut diviser par (mp(il, i) * sigd |
230 |
! * rg) et non par (mp(il, i) + sigd * 0.1). Et il |
! * rg) et non par (mp(il, i) + sigd * 0.1). Et il |
231 |
! faut bien revoir les facteurs 100. |
! faut bien revoir les facteurs 100. |
232 |
b(il, i - 1) = b(il, i) + 100. * (p(il, i - 1) - p(il, i)) & |
b(il, i - 1) = max(b(il, i) + 100. * (p(il, i - 1) & |
233 |
* tevap / (mp(il, i) + sigd * 0.1) - 10. * (th(il, i) & |
- p(il, i)) * tevap / (mp(il, i) + sigd * 0.1) - 10. & |
234 |
- th(il, i - 1)) * t(il, i) / (lvcp(il, i) * sigd & |
* (th(il, i) - th(il, i - 1)) * t(il, i) & |
235 |
* th(il, i)) |
/ (lvcp(il, i) * sigd * th(il, i)), 0.) |
|
b(il, i - 1) = max(b(il, i - 1), 0.) |
|
236 |
endif |
endif |
237 |
|
|
238 |
! Limit magnitude of mp(i) to meet CFL condition: |
! Limit magnitude of mp to meet CFL condition: |
239 |
ampmax = 2. * (ph(il, i) - ph(il, i + 1)) * delti |
ampmax = 2. * (ph(il, i) - ph(il, i + 1)) * delti |
240 |
amp2 = 2. * (ph(il, i - 1) - ph(il, i)) * delti |
ampmax = min(ampmax, 2. * (ph(il, i - 1) - ph(il, i)) * delti) |
|
ampmax = min(ampmax, amp2) |
|
241 |
mp(il, i) = min(mp(il, i), ampmax) |
mp(il, i) = min(mp(il, i), ampmax) |
242 |
|
|
243 |
! Force mp to decrease linearly to zero between cloud |
! Force mp to decrease linearly to zero between cloud |