1 |
SUBROUTINE cv_yield(nloc,ncum,nd,nk,icb,inb,delt & |
2 |
,t,q,u,v,gz,p,ph,h,hp,lv,cpn & |
3 |
,ep,clw,frac,m,mp,qp,up,vp & |
4 |
,wt,water,evap & |
5 |
,ment,qent,uent,vent,nent,elij & |
6 |
,tv,tvp & |
7 |
,iflag,wd,qprime,tprime & |
8 |
,precip,cbmf,ft,fq,fu,fv,Ma,qcondc) |
9 |
use cvthermo |
10 |
use cv_param |
11 |
implicit none |
12 |
|
13 |
|
14 |
! inputs |
15 |
integer, intent(in):: ncum, nd, nloc |
16 |
integer nk(nloc), icb(nloc), inb(nloc) |
17 |
integer nent(nloc,nd) |
18 |
real, intent(in):: delt |
19 |
real t(nloc,nd), q(nloc,nd), u(nloc,nd), v(nloc,nd) |
20 |
real gz(nloc,nd) |
21 |
real p(nloc,nd), ph(nloc,nd+1), h(nloc,nd) |
22 |
real hp(nloc,nd), lv(nloc,nd) |
23 |
real cpn(nloc,nd), ep(nloc,nd), clw(nloc,nd), frac(nloc) |
24 |
real m(nloc,nd), mp(nloc,nd), qp(nloc,nd) |
25 |
real up(nloc,nd), vp(nloc,nd) |
26 |
real wt(nloc,nd), water(nloc,nd), evap(nloc,nd) |
27 |
real ment(nloc,nd,nd), qent(nloc,nd,nd), elij(nloc,nd,nd) |
28 |
real uent(nloc,nd,nd), vent(nloc,nd,nd) |
29 |
real tv(nloc,nd), tvp(nloc,nd) |
30 |
|
31 |
! outputs |
32 |
integer iflag(nloc) ! also an input |
33 |
real cbmf(nloc) ! also an input |
34 |
real wd(nloc), tprime(nloc), qprime(nloc) |
35 |
real precip(nloc) |
36 |
real ft(nloc,nd), fq(nloc,nd), fu(nloc,nd), fv(nloc,nd) |
37 |
real Ma(nloc,nd) |
38 |
real qcondc(nloc,nd) |
39 |
|
40 |
! local variables |
41 |
integer i,j,ij,k,num1 |
42 |
real dpinv,cpinv,awat,fqold,ftold,fuold,fvold,delti |
43 |
real work(nloc), am(nloc),amp1(nloc),ad(nloc) |
44 |
real ents(nloc), uav(nloc),vav(nloc),lvcp(nloc,nd) |
45 |
real qcond(nloc,nd), nqcond(nloc,nd), wa(nloc,nd) ! cld |
46 |
real siga(nloc,nd), ax(nloc,nd), mac(nloc,nd) ! cld |
47 |
|
48 |
|
49 |
! -- initializations: |
50 |
|
51 |
delti = 1.0/delt |
52 |
|
53 |
do 160 i=1,ncum |
54 |
precip(i)=0.0 |
55 |
wd(i)=0.0 |
56 |
tprime(i)=0.0 |
57 |
qprime(i)=0.0 |
58 |
do 170 k=1,nl+1 |
59 |
ft(i,k)=0.0 |
60 |
fu(i,k)=0.0 |
61 |
fv(i,k)=0.0 |
62 |
fq(i,k)=0.0 |
63 |
lvcp(i,k)=lv(i,k)/cpn(i,k) |
64 |
qcondc(i,k)=0.0 ! cld |
65 |
qcond(i,k)=0.0 ! cld |
66 |
nqcond(i,k)=0.0 ! cld |
67 |
170 continue |
68 |
160 continue |
69 |
|
70 |
! |
71 |
! *** Calculate surface precipitation in mm/day *** |
72 |
! |
73 |
do 1190 i=1,ncum |
74 |
if(iflag(i).le.1)then |
75 |
precip(i) = wt(i,1)*sigd*water(i,1)*86400/g |
76 |
endif |
77 |
1190 continue |
78 |
! |
79 |
! |
80 |
! *** Calculate downdraft velocity scale and surface temperature and *** |
81 |
! *** water vapor fluctuations *** |
82 |
! |
83 |
do i=1,ncum |
84 |
wd(i)=betad*abs(mp(i,icb(i)))*0.01*rrd*t(i,icb(i)) & |
85 |
/(sigd*p(i,icb(i))) |
86 |
qprime(i)=0.5*(qp(i,1)-q(i,1)) |
87 |
tprime(i)=lv0*qprime(i)/cpd |
88 |
enddo |
89 |
! |
90 |
! *** Calculate tendencies of lowest level potential temperature *** |
91 |
! *** and mixing ratio *** |
92 |
! |
93 |
do 1200 i=1,ncum |
94 |
work(i)=0.01/(ph(i,1)-ph(i,2)) |
95 |
am(i)=0.0 |
96 |
1200 continue |
97 |
do 1220 k=2,nl |
98 |
do 1210 i=1,ncum |
99 |
if((nk(i).eq.1).and.(k.le.inb(i)).and.(nk(i).eq.1))then |
100 |
am(i)=am(i)+m(i,k) |
101 |
endif |
102 |
1210 continue |
103 |
1220 continue |
104 |
do 1240 i=1,ncum |
105 |
if((g*work(i)*am(i)).ge.delti)iflag(i)=1 |
106 |
ft(i,1)=ft(i,1)+g*work(i)*am(i)*(t(i,2)-t(i,1) & |
107 |
+(gz(i,2)-gz(i,1))/cpn(i,1)) |
108 |
ft(i,1)=ft(i,1)-lvcp(i,1)*sigd*evap(i,1) |
109 |
ft(i,1)=ft(i,1)+sigd*wt(i,2)*(cl-cpd)*water(i,2)*(t(i,2) & |
110 |
-t(i,1))*work(i)/cpn(i,1) |
111 |
fq(i,1)=fq(i,1)+g*mp(i,2)*(qp(i,2)-q(i,1))* & |
112 |
work(i)+sigd*evap(i,1) |
113 |
fq(i,1)=fq(i,1)+g*am(i)*(q(i,2)-q(i,1))*work(i) |
114 |
fu(i,1)=fu(i,1)+g*work(i)*(mp(i,2)*(up(i,2)-u(i,1)) & |
115 |
+am(i)*(u(i,2)-u(i,1))) |
116 |
fv(i,1)=fv(i,1)+g*work(i)*(mp(i,2)*(vp(i,2)-v(i,1)) & |
117 |
+am(i)*(v(i,2)-v(i,1))) |
118 |
1240 continue |
119 |
do 1260 j=2,nl |
120 |
do 1250 i=1,ncum |
121 |
if(j.le.inb(i))then |
122 |
fq(i,1)=fq(i,1) & |
123 |
+g*work(i)*ment(i,j,1)*(qent(i,j,1)-q(i,1)) |
124 |
fu(i,1)=fu(i,1) & |
125 |
+g*work(i)*ment(i,j,1)*(uent(i,j,1)-u(i,1)) |
126 |
fv(i,1)=fv(i,1) & |
127 |
+g*work(i)*ment(i,j,1)*(vent(i,j,1)-v(i,1)) |
128 |
endif |
129 |
1250 continue |
130 |
1260 continue |
131 |
! |
132 |
! *** Calculate tendencies of potential temperature and mixing ratio *** |
133 |
! *** at levels above the lowest level *** |
134 |
! |
135 |
! *** First find the net saturated updraft and downdraft mass fluxes *** |
136 |
! *** through each level *** |
137 |
! |
138 |
do 1500 i=2,nl+1 |
139 |
! |
140 |
num1=0 |
141 |
do 1265 ij=1,ncum |
142 |
if(i.le.inb(ij))num1=num1+1 |
143 |
1265 continue |
144 |
if(num1.le.0)go to 1500 |
145 |
! |
146 |
call zilch(amp1,ncum) |
147 |
call zilch(ad,ncum) |
148 |
! |
149 |
do 1280 k=i+1,nl+1 |
150 |
do 1270 ij=1,ncum |
151 |
if((i.ge.nk(ij)).and.(i.le.inb(ij)) & |
152 |
.and.(k.le.(inb(ij)+1)))then |
153 |
amp1(ij)=amp1(ij)+m(ij,k) |
154 |
endif |
155 |
1270 continue |
156 |
1280 continue |
157 |
! |
158 |
do 1310 k=1,i |
159 |
do 1300 j=i+1,nl+1 |
160 |
do 1290 ij=1,ncum |
161 |
if((j.le.(inb(ij)+1)).and.(i.le.inb(ij)))then |
162 |
amp1(ij)=amp1(ij)+ment(ij,k,j) |
163 |
endif |
164 |
1290 continue |
165 |
1300 continue |
166 |
1310 continue |
167 |
do 1340 k=1,i-1 |
168 |
do 1330 j=i,nl+1 |
169 |
do 1320 ij=1,ncum |
170 |
if((i.le.inb(ij)).and.(j.le.inb(ij)))then |
171 |
ad(ij)=ad(ij)+ment(ij,j,k) |
172 |
endif |
173 |
1320 continue |
174 |
1330 continue |
175 |
1340 continue |
176 |
! |
177 |
do 1350 ij=1,ncum |
178 |
if(i.le.inb(ij))then |
179 |
dpinv=0.01/(ph(ij,i)-ph(ij,i+1)) |
180 |
cpinv=1.0/cpn(ij,i) |
181 |
! |
182 |
ft(ij,i)=ft(ij,i) & |
183 |
+g*dpinv*(amp1(ij)*(t(ij,i+1)-t(ij,i) & |
184 |
+(gz(ij,i+1)-gz(ij,i))*cpinv) & |
185 |
-ad(ij)*(t(ij,i)-t(ij,i-1)+(gz(ij,i)-gz(ij,i-1))*cpinv)) & |
186 |
-sigd*lvcp(ij,i)*evap(ij,i) |
187 |
ft(ij,i)=ft(ij,i)+g*dpinv*ment(ij,i,i)*(hp(ij,i)-h(ij,i)+ & |
188 |
t(ij,i)*(cpv-cpd)*(q(ij,i)-qent(ij,i,i)))*cpinv |
189 |
ft(ij,i)=ft(ij,i)+sigd*wt(ij,i+1)*(cl-cpd)*water(ij,i+1)* & |
190 |
(t(ij,i+1)-t(ij,i))*dpinv*cpinv |
191 |
fq(ij,i)=fq(ij,i)+g*dpinv*(amp1(ij)*(q(ij,i+1)-q(ij,i))- & |
192 |
ad(ij)*(q(ij,i)-q(ij,i-1))) |
193 |
fu(ij,i)=fu(ij,i)+g*dpinv*(amp1(ij)*(u(ij,i+1)-u(ij,i))- & |
194 |
ad(ij)*(u(ij,i)-u(ij,i-1))) |
195 |
fv(ij,i)=fv(ij,i)+g*dpinv*(amp1(ij)*(v(ij,i+1)-v(ij,i))- & |
196 |
ad(ij)*(v(ij,i)-v(ij,i-1))) |
197 |
endif |
198 |
1350 continue |
199 |
do 1370 k=1,i-1 |
200 |
do 1360 ij=1,ncum |
201 |
if(i.le.inb(ij))then |
202 |
awat=elij(ij,k,i)-(1.-ep(ij,i))*clw(ij,i) |
203 |
awat=max(awat,0.0) |
204 |
fq(ij,i)=fq(ij,i) & |
205 |
+g*dpinv*ment(ij,k,i)*(qent(ij,k,i)-awat-q(ij,i)) |
206 |
fu(ij,i)=fu(ij,i) & |
207 |
+g*dpinv*ment(ij,k,i)*(uent(ij,k,i)-u(ij,i)) |
208 |
fv(ij,i)=fv(ij,i) & |
209 |
+g*dpinv*ment(ij,k,i)*(vent(ij,k,i)-v(ij,i)) |
210 |
! (saturated updrafts resulting from mixing) ! cld |
211 |
qcond(ij,i)=qcond(ij,i)+(elij(ij,k,i)-awat) ! cld |
212 |
nqcond(ij,i)=nqcond(ij,i)+1. ! cld |
213 |
endif |
214 |
1360 continue |
215 |
1370 continue |
216 |
do 1390 k=i,nl+1 |
217 |
do 1380 ij=1,ncum |
218 |
if((i.le.inb(ij)).and.(k.le.inb(ij)))then |
219 |
fq(ij,i)=fq(ij,i) & |
220 |
+g*dpinv*ment(ij,k,i)*(qent(ij,k,i)-q(ij,i)) |
221 |
fu(ij,i)=fu(ij,i) & |
222 |
+g*dpinv*ment(ij,k,i)*(uent(ij,k,i)-u(ij,i)) |
223 |
fv(ij,i)=fv(ij,i) & |
224 |
+g*dpinv*ment(ij,k,i)*(vent(ij,k,i)-v(ij,i)) |
225 |
endif |
226 |
1380 continue |
227 |
1390 continue |
228 |
do 1400 ij=1,ncum |
229 |
if(i.le.inb(ij))then |
230 |
fq(ij,i)=fq(ij,i) & |
231 |
+sigd*evap(ij,i)+g*(mp(ij,i+1)* & |
232 |
(qp(ij,i+1)-q(ij,i)) & |
233 |
-mp(ij,i)*(qp(ij,i)-q(ij,i-1)))*dpinv |
234 |
fu(ij,i)=fu(ij,i) & |
235 |
+g*(mp(ij,i+1)*(up(ij,i+1)-u(ij,i))-mp(ij,i)* & |
236 |
(up(ij,i)-u(ij,i-1)))*dpinv |
237 |
fv(ij,i)=fv(ij,i) & |
238 |
+g*(mp(ij,i+1)*(vp(ij,i+1)-v(ij,i))-mp(ij,i)* & |
239 |
(vp(ij,i)-v(ij,i-1)))*dpinv |
240 |
! (saturated downdrafts resulting from mixing) ! cld |
241 |
do k=i+1,inb(ij) ! cld |
242 |
qcond(ij,i)=qcond(ij,i)+elij(ij,k,i) ! cld |
243 |
nqcond(ij,i)=nqcond(ij,i)+1. ! cld |
244 |
enddo ! cld |
245 |
! (particular case: no detraining level is found) ! cld |
246 |
if (nent(ij,i).eq.0) then ! cld |
247 |
qcond(ij,i)=qcond(ij,i)+(1.-ep(ij,i))*clw(ij,i) ! cld |
248 |
nqcond(ij,i)=nqcond(ij,i)+1. ! cld |
249 |
endif ! cld |
250 |
if (nqcond(ij,i).ne.0.) then ! cld |
251 |
qcond(ij,i)=qcond(ij,i)/nqcond(ij,i) ! cld |
252 |
endif ! cld |
253 |
endif |
254 |
1400 continue |
255 |
1500 continue |
256 |
! |
257 |
! *** Adjust tendencies at top of convection layer to reflect *** |
258 |
! *** actual position of the level zero cape *** |
259 |
! |
260 |
do 503 ij=1,ncum |
261 |
fqold=fq(ij,inb(ij)) |
262 |
fq(ij,inb(ij))=fq(ij,inb(ij))*(1.-frac(ij)) |
263 |
fq(ij,inb(ij)-1)=fq(ij,inb(ij)-1) & |
264 |
+frac(ij)*fqold*((ph(ij,inb(ij))-ph(ij,inb(ij)+1))/ & |
265 |
(ph(ij,inb(ij)-1)-ph(ij,inb(ij))))*lv(ij,inb(ij)) & |
266 |
/lv(ij,inb(ij)-1) |
267 |
ftold=ft(ij,inb(ij)) |
268 |
ft(ij,inb(ij))=ft(ij,inb(ij))*(1.-frac(ij)) |
269 |
ft(ij,inb(ij)-1)=ft(ij,inb(ij)-1) & |
270 |
+frac(ij)*ftold*((ph(ij,inb(ij))-ph(ij,inb(ij)+1))/ & |
271 |
(ph(ij,inb(ij)-1)-ph(ij,inb(ij))))*cpn(ij,inb(ij)) & |
272 |
/cpn(ij,inb(ij)-1) |
273 |
fuold=fu(ij,inb(ij)) |
274 |
fu(ij,inb(ij))=fu(ij,inb(ij))*(1.-frac(ij)) |
275 |
fu(ij,inb(ij)-1)=fu(ij,inb(ij)-1) & |
276 |
+frac(ij)*fuold*((ph(ij,inb(ij))-ph(ij,inb(ij)+1))/ & |
277 |
(ph(ij,inb(ij)-1)-ph(ij,inb(ij)))) |
278 |
fvold=fv(ij,inb(ij)) |
279 |
fv(ij,inb(ij))=fv(ij,inb(ij))*(1.-frac(ij)) |
280 |
fv(ij,inb(ij)-1)=fv(ij,inb(ij)-1) & |
281 |
+frac(ij)*fvold*((ph(ij,inb(ij))-ph(ij,inb(ij)+1))/ & |
282 |
(ph(ij,inb(ij)-1)-ph(ij,inb(ij)))) |
283 |
503 continue |
284 |
! |
285 |
! *** Very slightly adjust tendencies to force exact *** |
286 |
! *** enthalpy, momentum and tracer conservation *** |
287 |
! |
288 |
do 682 ij=1,ncum |
289 |
ents(ij)=0.0 |
290 |
uav(ij)=0.0 |
291 |
vav(ij)=0.0 |
292 |
do 681 i=1,inb(ij) |
293 |
ents(ij)=ents(ij) & |
294 |
+(cpn(ij,i)*ft(ij,i)+lv(ij,i)*fq(ij,i))*(ph(ij,i)-ph(ij,i+1)) |
295 |
uav(ij)=uav(ij)+fu(ij,i)*(ph(ij,i)-ph(ij,i+1)) |
296 |
vav(ij)=vav(ij)+fv(ij,i)*(ph(ij,i)-ph(ij,i+1)) |
297 |
681 continue |
298 |
682 continue |
299 |
do 683 ij=1,ncum |
300 |
ents(ij)=ents(ij)/(ph(ij,1)-ph(ij,inb(ij)+1)) |
301 |
uav(ij)=uav(ij)/(ph(ij,1)-ph(ij,inb(ij)+1)) |
302 |
vav(ij)=vav(ij)/(ph(ij,1)-ph(ij,inb(ij)+1)) |
303 |
683 continue |
304 |
do 642 ij=1,ncum |
305 |
do 641 i=1,inb(ij) |
306 |
ft(ij,i)=ft(ij,i)-ents(ij)/cpn(ij,i) |
307 |
fu(ij,i)=(1.-cu)*(fu(ij,i)-uav(ij)) |
308 |
fv(ij,i)=(1.-cu)*(fv(ij,i)-vav(ij)) |
309 |
641 continue |
310 |
642 continue |
311 |
! |
312 |
do 1810 k=1,nl+1 |
313 |
do 1800 i=1,ncum |
314 |
if((q(i,k)+delt*fq(i,k)).lt.0.0)iflag(i)=10 |
315 |
1800 continue |
316 |
1810 continue |
317 |
! |
318 |
! |
319 |
do 1900 i=1,ncum |
320 |
if(iflag(i).gt.2)then |
321 |
precip(i)=0.0 |
322 |
cbmf(i)=0.0 |
323 |
endif |
324 |
1900 continue |
325 |
do 1920 k=1,nl |
326 |
do 1910 i=1,ncum |
327 |
if(iflag(i).gt.2)then |
328 |
ft(i,k)=0.0 |
329 |
fq(i,k)=0.0 |
330 |
fu(i,k)=0.0 |
331 |
fv(i,k)=0.0 |
332 |
qcondc(i,k)=0.0 ! cld |
333 |
endif |
334 |
1910 continue |
335 |
1920 continue |
336 |
|
337 |
do k=1,nl+1 |
338 |
do i=1,ncum |
339 |
Ma(i,k) = 0. |
340 |
enddo |
341 |
enddo |
342 |
do k=nl,1,-1 |
343 |
do i=1,ncum |
344 |
Ma(i,k) = Ma(i,k+1)+m(i,k) |
345 |
enddo |
346 |
enddo |
347 |
|
348 |
! |
349 |
! *** diagnose the in-cloud mixing ratio *** |
350 |
! *** of condensed water *** |
351 |
! |
352 |
DO ij=1,ncum |
353 |
do i=1,nd |
354 |
mac(ij,i)=0.0 |
355 |
wa(ij,i)=0.0 |
356 |
siga(ij,i)=0.0 |
357 |
enddo |
358 |
do i=nk(ij),inb(ij) |
359 |
do k=i+1,inb(ij)+1 |
360 |
mac(ij,i)=mac(ij,i)+m(ij,k) |
361 |
enddo |
362 |
enddo |
363 |
do i=icb(ij),inb(ij)-1 |
364 |
ax(ij,i)=0. |
365 |
do j=icb(ij),i |
366 |
ax(ij,i)=ax(ij,i)+rrd*(tvp(ij,j)-tv(ij,j)) & |
367 |
*(ph(ij,j)-ph(ij,j+1))/p(ij,j) |
368 |
enddo |
369 |
if (ax(ij,i).gt.0.0) then |
370 |
wa(ij,i)=sqrt(2.*ax(ij,i)) |
371 |
endif |
372 |
enddo |
373 |
do i=1,nl |
374 |
if (wa(ij,i).gt.0.0) & |
375 |
siga(ij,i)=mac(ij,i)/wa(ij,i) & |
376 |
*rrd*tvp(ij,i)/p(ij,i)/100./delta |
377 |
siga(ij,i) = min(siga(ij,i),1.0) |
378 |
qcondc(ij,i)=siga(ij,i)*clw(ij,i)*(1.-ep(ij,i)) & |
379 |
+ (1.-siga(ij,i))*qcond(ij,i) |
380 |
enddo |
381 |
ENDDO |
382 |
|
383 |
return |
384 |
end |