1 |
guez |
185 |
module cv30_undilute2_m |
2 |
guez |
47 |
|
3 |
guez |
183 |
implicit none |
4 |
guez |
47 |
|
5 |
guez |
183 |
contains |
6 |
guez |
47 |
|
7 |
guez |
187 |
SUBROUTINE cv30_undilute2(ncum, icb, icbs, nk, tnk, qnk, gznk, t, qs, gz, & |
8 |
|
|
p, h, tv, lv, pbase, buoybase, plcl, inb, tp, tvp, clw, hp, ep, sigp, & |
9 |
|
|
buoy) |
10 |
guez |
47 |
|
11 |
guez |
187 |
! Undilute (adiabatic) updraft, second part. Purpose: find the |
12 |
|
|
! rest of the lifted parcel temperatures; compute the |
13 |
|
|
! precipitation efficiencies and the fraction of precipitation |
14 |
|
|
! falling outside of cloud; find the level of neutral buoyancy. |
15 |
guez |
47 |
|
16 |
guez |
183 |
! Vertical profile of buoyancy computed here (use of buoybase). |
17 |
guez |
47 |
|
18 |
guez |
183 |
use conema3_m, only: epmax |
19 |
guez |
186 |
use cv30_param_m, only: dtovsh, minorig, nl, pbcrit, ptcrit, spfac |
20 |
guez |
190 |
use cv_thermo_m, only: cl, clmcpv, cpd, cpv, eps, lv0, rrv |
21 |
guez |
187 |
USE dimphy, ONLY: klon, klev |
22 |
guez |
47 |
|
23 |
guez |
187 |
integer, intent(in):: ncum |
24 |
guez |
47 |
|
25 |
guez |
187 |
integer, intent(in):: icb(klon), icbs(klon) |
26 |
|
|
! icbs is the first level above LCL (may differ from icb) |
27 |
|
|
|
28 |
|
|
integer, intent(in):: nk(klon) |
29 |
|
|
real, intent(in):: tnk(klon), qnk(klon), gznk(klon) |
30 |
|
|
real, intent(in):: t(klon, klev), qs(klon, klev), gz(klon, klev) |
31 |
|
|
real, intent(in):: p(klon, klev), h(klon, klev) |
32 |
|
|
real, intent(in):: tv(klon, klev), lv(klon, klev) |
33 |
|
|
real, intent(in):: pbase(klon), buoybase(klon), plcl(klon) |
34 |
|
|
|
35 |
guez |
183 |
! outputs: |
36 |
guez |
187 |
integer, intent(out):: inb(:) ! (ncum) |
37 |
|
|
! first model level above the level of neutral buoyancy of the |
38 |
|
|
! parcel (<= nl - 1) |
39 |
|
|
|
40 |
|
|
real tp(klon, klev), tvp(klon, klev), clw(klon, klev) |
41 |
guez |
183 |
! condensed water not removed from tvp |
42 |
guez |
187 |
real hp(klon, klev), ep(klon, klev), sigp(klon, klev) |
43 |
|
|
real buoy(klon, klev) |
44 |
guez |
47 |
|
45 |
guez |
183 |
! Local: |
46 |
|
|
integer i, k |
47 |
|
|
real tg, qg, ahg, alv, s, tc, es, denom |
48 |
|
|
real pden |
49 |
guez |
187 |
real ah0(klon) |
50 |
guez |
47 |
|
51 |
guez |
183 |
!--------------------------------------------------------------------- |
52 |
guez |
47 |
|
53 |
guez |
183 |
! SOME INITIALIZATIONS |
54 |
guez |
47 |
|
55 |
guez |
183 |
do k = 1, nl |
56 |
|
|
do i = 1, ncum |
57 |
|
|
ep(i, k) = 0.0 |
58 |
|
|
sigp(i, k) = spfac |
59 |
|
|
end do |
60 |
|
|
end do |
61 |
guez |
47 |
|
62 |
guez |
183 |
! FIND THE REST OF THE LIFTED PARCEL TEMPERATURES |
63 |
guez |
47 |
|
64 |
guez |
183 |
! The procedure is to solve the equation. |
65 |
|
|
! cp * tp + L * qp + phi = cp * tnk + L * qnk + gznk. |
66 |
guez |
47 |
|
67 |
guez |
183 |
! Calculate certain parcel quantities, including static energy |
68 |
guez |
47 |
|
69 |
guez |
183 |
do i = 1, ncum |
70 |
|
|
ah0(i) = (cpd * (1. - qnk(i)) + cl * qnk(i)) * tnk(i) & |
71 |
|
|
+ qnk(i) * (lv0 - clmcpv * (tnk(i) - 273.15)) + gznk(i) |
72 |
|
|
end do |
73 |
guez |
47 |
|
74 |
guez |
183 |
! Find lifted parcel quantities above cloud base |
75 |
guez |
47 |
|
76 |
guez |
183 |
do k = minorig + 1, nl |
77 |
|
|
do i = 1, ncum |
78 |
|
|
if (k >= (icbs(i) + 1)) then |
79 |
|
|
tg = t(i, k) |
80 |
|
|
qg = qs(i, k) |
81 |
|
|
alv = lv0 - clmcpv * (t(i, k) - 273.15) |
82 |
guez |
47 |
|
83 |
guez |
183 |
! First iteration. |
84 |
guez |
47 |
|
85 |
guez |
183 |
s = cpd * (1. - qnk(i)) + cl * qnk(i) & |
86 |
|
|
+ alv * alv * qg / (rrv * t(i, k) * t(i, k)) |
87 |
|
|
s = 1. / s |
88 |
guez |
47 |
|
89 |
guez |
183 |
ahg = cpd * tg + (cl - cpd) * qnk(i) * tg + alv * qg + gz(i, k) |
90 |
|
|
tg = tg + s * (ah0(i) - ahg) |
91 |
guez |
47 |
|
92 |
guez |
183 |
tc = tg - 273.15 |
93 |
|
|
denom = 243.5 + tc |
94 |
|
|
denom = MAX(denom, 1.0) |
95 |
guez |
47 |
|
96 |
guez |
183 |
es = 6.112 * exp(17.67 * tc / denom) |
97 |
guez |
47 |
|
98 |
guez |
183 |
qg = eps * es / (p(i, k) - es * (1. - eps)) |
99 |
guez |
47 |
|
100 |
guez |
183 |
! Second iteration. |
101 |
guez |
47 |
|
102 |
guez |
183 |
ahg = cpd * tg + (cl - cpd) * qnk(i) * tg + alv * qg + gz(i, k) |
103 |
|
|
tg = tg + s * (ah0(i) - ahg) |
104 |
guez |
47 |
|
105 |
guez |
183 |
tc = tg - 273.15 |
106 |
|
|
denom = 243.5 + tc |
107 |
|
|
denom = MAX(denom, 1.0) |
108 |
guez |
47 |
|
109 |
guez |
183 |
es = 6.112 * exp(17.67 * tc / denom) |
110 |
|
|
|
111 |
|
|
qg = eps * es / (p(i, k) - es * (1. - eps)) |
112 |
|
|
|
113 |
|
|
alv = lv0 - clmcpv * (t(i, k) - 273.15) |
114 |
|
|
|
115 |
|
|
! no approximation: |
116 |
|
|
tp(i, k) = (ah0(i) - gz(i, k) - alv * qg) & |
117 |
|
|
/ (cpd + (cl - cpd) * qnk(i)) |
118 |
|
|
|
119 |
|
|
clw(i, k) = qnk(i) - qg |
120 |
|
|
clw(i, k) = max(0.0, clw(i, k)) |
121 |
|
|
! qg utilise au lieu du vrai mixing ratio rg: |
122 |
|
|
tvp(i, k) = tp(i, k) * (1. + qg / eps - qnk(i)) ! whole thing |
123 |
|
|
endif |
124 |
|
|
end do |
125 |
|
|
end do |
126 |
|
|
|
127 |
|
|
! SET THE PRECIPITATION EFFICIENCIES AND THE FRACTION OF |
128 |
|
|
! PRECIPITATION FALLING OUTSIDE OF CLOUD |
129 |
|
|
! THESE MAY BE FUNCTIONS OF TP(I), P(I) AND CLW(I) |
130 |
|
|
do k = 1, nl |
131 |
|
|
do i = 1, ncum |
132 |
|
|
pden = ptcrit - pbcrit |
133 |
|
|
ep(i, k) = (plcl(i) - p(i, k) - pbcrit) / pden * epmax |
134 |
|
|
ep(i, k) = amax1(ep(i, k), 0.0) |
135 |
|
|
ep(i, k) = amin1(ep(i, k), epmax) |
136 |
|
|
sigp(i, k) = spfac |
137 |
|
|
end do |
138 |
|
|
end do |
139 |
|
|
|
140 |
|
|
! CALCULATE VIRTUAL TEMPERATURE AND LIFTED PARCEL |
141 |
|
|
! VIRTUAL TEMPERATURE |
142 |
|
|
|
143 |
|
|
! tvp est calcule en une seule fois, et sans retirer |
144 |
|
|
! l'eau condensee (~> reversible CAPE) |
145 |
|
|
do i = 1, ncum |
146 |
guez |
186 |
tp(i, nl + 1) = tp(i, nl) |
147 |
guez |
183 |
end do |
148 |
|
|
|
149 |
|
|
! EFFECTIVE VERTICAL PROFILE OF BUOYANCY: |
150 |
|
|
|
151 |
|
|
! first estimate of buoyancy: |
152 |
|
|
do i = 1, ncum |
153 |
|
|
do k = 1, nl |
154 |
|
|
buoy(i, k) = tvp(i, k) - tv(i, k) |
155 |
|
|
end do |
156 |
|
|
end do |
157 |
|
|
|
158 |
|
|
! set buoyancy = buoybase for all levels below base |
159 |
|
|
! for safety, set buoy(icb) = buoybase |
160 |
|
|
do i = 1, ncum |
161 |
|
|
do k = 1, nl |
162 |
|
|
if ((k >= icb(i)) .and. (k <= nl) .and. (p(i, k) >= pbase(i))) then |
163 |
|
|
buoy(i, k) = buoybase(i) |
164 |
|
|
endif |
165 |
|
|
end do |
166 |
|
|
buoy(icb(i), k) = buoybase(i) |
167 |
|
|
end do |
168 |
|
|
|
169 |
guez |
187 |
! Compute inb: |
170 |
guez |
183 |
|
171 |
guez |
187 |
inb = nl - 1 |
172 |
guez |
183 |
|
173 |
|
|
do i = 1, ncum |
174 |
|
|
do k = 1, nl - 1 |
175 |
|
|
if ((k >= icb(i)) .and. (buoy(i, k) < dtovsh)) then |
176 |
|
|
inb(i) = MIN(inb(i), k) |
177 |
|
|
endif |
178 |
|
|
end do |
179 |
|
|
end do |
180 |
|
|
|
181 |
|
|
! CALCULATE LIQUID WATER STATIC ENERGY OF LIFTED PARCEL |
182 |
|
|
|
183 |
guez |
186 |
do k = 1, nl + 1 |
184 |
guez |
183 |
do i = 1, ncum |
185 |
|
|
hp(i, k) = h(i, k) |
186 |
|
|
enddo |
187 |
|
|
enddo |
188 |
|
|
|
189 |
|
|
do k = minorig + 1, nl |
190 |
|
|
do i = 1, ncum |
191 |
|
|
if (k >= icb(i) .and. k <= inb(i)) hp(i, k) = h(i, nk(i)) & |
192 |
|
|
+ (lv(i, k) + (cpd - cpv) * t(i, k)) * ep(i, k) * clw(i, k) |
193 |
|
|
end do |
194 |
|
|
end do |
195 |
|
|
|
196 |
guez |
185 |
end SUBROUTINE cv30_undilute2 |
197 |
guez |
183 |
|
198 |
guez |
185 |
end module cv30_undilute2_m |