1 |
|
2 |
SUBROUTINE cv_mixing(nloc,ncum,nd,icb,nk,inb,inb1 & |
3 |
,ph,t,q,qs,u,v,h,lv,qnk & |
4 |
,hp,tv,tvp,ep,clw,cbmf & |
5 |
,m,ment,qent,uent,vent,nent,sij,elij) |
6 |
use cvthermo |
7 |
use cvparam |
8 |
implicit none |
9 |
|
10 |
|
11 |
! inputs: |
12 |
integer ncum, nd, nloc |
13 |
integer icb(nloc), inb(nloc), inb1(nloc), nk(nloc) |
14 |
real cbmf(nloc), qnk(nloc) |
15 |
real ph(nloc,nd+1) |
16 |
real t(nloc,nd), q(nloc,nd), qs(nloc,nd), lv(nloc,nd) |
17 |
real u(nloc,nd), v(nloc,nd), h(nloc,nd), hp(nloc,nd) |
18 |
real tv(nloc,nd), tvp(nloc,nd), ep(nloc,nd), clw(nloc,nd) |
19 |
|
20 |
! outputs: |
21 |
integer nent(nloc,nd) |
22 |
real m(nloc,nd), ment(nloc,nd,nd), qent(nloc,nd,nd) |
23 |
real uent(nloc,nd,nd), vent(nloc,nd,nd) |
24 |
real sij(nloc,nd,nd), elij(nloc,nd,nd) |
25 |
|
26 |
! local variables: |
27 |
integer i, j, k, ij |
28 |
integer num1, num2 |
29 |
real dbo, qti, bf2, anum, denom, dei, altem, cwat, stemp |
30 |
real alt, qp1, smid, sjmin, sjmax, delp, delm |
31 |
real work(nloc), asij(nloc), smin(nloc), scrit(nloc) |
32 |
real bsum(nloc,nd) |
33 |
logical lwork(nloc) |
34 |
|
35 |
!===================================================================== |
36 |
! --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS |
37 |
!===================================================================== |
38 |
! |
39 |
do 360 i=1,ncum*nlp |
40 |
nent(i,1)=0 |
41 |
m(i,1)=0.0 |
42 |
360 continue |
43 |
! |
44 |
do 400 k=1,nlp |
45 |
do 390 j=1,nlp |
46 |
do 385 i=1,ncum |
47 |
qent(i,k,j)=q(i,j) |
48 |
uent(i,k,j)=u(i,j) |
49 |
vent(i,k,j)=v(i,j) |
50 |
elij(i,k,j)=0.0 |
51 |
ment(i,k,j)=0.0 |
52 |
sij(i,k,j)=0.0 |
53 |
385 continue |
54 |
390 continue |
55 |
400 continue |
56 |
! |
57 |
!------------------------------------------------------------------- |
58 |
! --- Calculate rates of mixing, m(i) |
59 |
!------------------------------------------------------------------- |
60 |
! |
61 |
call zilch(work,ncum) |
62 |
! |
63 |
do 670 j=minorig+1,nl |
64 |
do 660 i=1,ncum |
65 |
if((j.ge.(icb(i)+1)).and.(j.le.inb(i)))then |
66 |
k=min(j,inb1(i)) |
67 |
dbo=abs(tv(i,k+1)-tvp(i,k+1)-tv(i,k-1)+tvp(i,k-1)) & |
68 |
+entp*0.04*(ph(i,k)-ph(i,k+1)) |
69 |
work(i)=work(i)+dbo |
70 |
m(i,j)=cbmf(i)*dbo |
71 |
endif |
72 |
660 continue |
73 |
670 continue |
74 |
do 690 k=minorig+1,nl |
75 |
do 680 i=1,ncum |
76 |
if((k.ge.(icb(i)+1)).and.(k.le.inb(i)))then |
77 |
m(i,k)=m(i,k)/work(i) |
78 |
endif |
79 |
680 continue |
80 |
690 continue |
81 |
! |
82 |
! |
83 |
!===================================================================== |
84 |
! --- CALCULATE ENTRAINED AIR MASS FLUX (ment), TOTAL WATER MIXING |
85 |
! --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING |
86 |
! --- FRACTION (sij) |
87 |
!===================================================================== |
88 |
! |
89 |
! |
90 |
do 750 i=minorig+1,nl |
91 |
do 710 j=minorig+1,nl |
92 |
do 700 ij=1,ncum |
93 |
if((i.ge.(icb(ij)+1)).and.(j.ge.icb(ij)) & |
94 |
.and.(i.le.inb(ij)).and.(j.le.inb(ij)))then |
95 |
qti=qnk(ij)-ep(ij,i)*clw(ij,i) |
96 |
bf2=1.+lv(ij,j)*lv(ij,j)*qs(ij,j) & |
97 |
/(rrv*t(ij,j)*t(ij,j)*cpd) |
98 |
anum=h(ij,j)-hp(ij,i)+(cpv-cpd)*t(ij,j)*(qti-q(ij,j)) |
99 |
denom=h(ij,i)-hp(ij,i)+(cpd-cpv)*(q(ij,i)-qti)*t(ij,j) |
100 |
dei=denom |
101 |
if(abs(dei).lt.0.01)dei=0.01 |
102 |
sij(ij,i,j)=anum/dei |
103 |
sij(ij,i,i)=1.0 |
104 |
altem=sij(ij,i,j)*q(ij,i)+(1.-sij(ij,i,j))*qti-qs(ij,j) |
105 |
altem=altem/bf2 |
106 |
cwat=clw(ij,j)*(1.-ep(ij,j)) |
107 |
stemp=sij(ij,i,j) |
108 |
if((stemp.lt.0.0.or.stemp.gt.1.0.or. & |
109 |
altem.gt.cwat).and.j.gt.i)then |
110 |
anum=anum-lv(ij,j)*(qti-qs(ij,j)-cwat*bf2) |
111 |
denom=denom+lv(ij,j)*(q(ij,i)-qti) |
112 |
if(abs(denom).lt.0.01)denom=0.01 |
113 |
sij(ij,i,j)=anum/denom |
114 |
altem=sij(ij,i,j)*q(ij,i)+(1.-sij(ij,i,j))*qti-qs(ij,j) |
115 |
altem=altem-(bf2-1.)*cwat |
116 |
endif |
117 |
if(sij(ij,i,j).gt.0.0.and.sij(ij,i,j).lt.0.9)then |
118 |
qent(ij,i,j)=sij(ij,i,j)*q(ij,i) & |
119 |
+(1.-sij(ij,i,j))*qti |
120 |
uent(ij,i,j)=sij(ij,i,j)*u(ij,i) & |
121 |
+(1.-sij(ij,i,j))*u(ij,nk(ij)) |
122 |
vent(ij,i,j)=sij(ij,i,j)*v(ij,i) & |
123 |
+(1.-sij(ij,i,j))*v(ij,nk(ij)) |
124 |
elij(ij,i,j)=altem |
125 |
elij(ij,i,j)=max(0.0,elij(ij,i,j)) |
126 |
ment(ij,i,j)=m(ij,i)/(1.-sij(ij,i,j)) |
127 |
nent(ij,i)=nent(ij,i)+1 |
128 |
endif |
129 |
sij(ij,i,j)=max(0.0,sij(ij,i,j)) |
130 |
sij(ij,i,j)=min(1.0,sij(ij,i,j)) |
131 |
endif |
132 |
700 continue |
133 |
710 continue |
134 |
! |
135 |
! *** If no air can entrain at level i assume that updraft detrains *** |
136 |
! *** at that level and calculate detrained air flux and properties *** |
137 |
! |
138 |
do 740 ij=1,ncum |
139 |
if((i.ge.(icb(ij)+1)).and.(i.le.inb(ij)) & |
140 |
.and.(nent(ij,i).eq.0))then |
141 |
ment(ij,i,i)=m(ij,i) |
142 |
qent(ij,i,i)=q(ij,nk(ij))-ep(ij,i)*clw(ij,i) |
143 |
uent(ij,i,i)=u(ij,nk(ij)) |
144 |
vent(ij,i,i)=v(ij,nk(ij)) |
145 |
elij(ij,i,i)=clw(ij,i) |
146 |
sij(ij,i,i)=1.0 |
147 |
endif |
148 |
740 continue |
149 |
750 continue |
150 |
! |
151 |
do 770 i=1,ncum |
152 |
sij(i,inb(i),inb(i))=1.0 |
153 |
770 continue |
154 |
! |
155 |
!===================================================================== |
156 |
! --- NORMALIZE ENTRAINED AIR MASS FLUXES |
157 |
! --- TO REPRESENT EQUAL PROBABILITIES OF MIXING |
158 |
!===================================================================== |
159 |
! |
160 |
call zilch(bsum,ncum*nlp) |
161 |
do 780 ij=1,ncum |
162 |
lwork(ij)=.false. |
163 |
780 continue |
164 |
do 789 i=minorig+1,nl |
165 |
! |
166 |
num1=0 |
167 |
do 779 ij=1,ncum |
168 |
if((i.ge.icb(ij)+1).and.(i.le.inb(ij)))num1=num1+1 |
169 |
779 continue |
170 |
if(num1.le.0)go to 789 |
171 |
! |
172 |
do 781 ij=1,ncum |
173 |
if((i.ge.icb(ij)+1).and.(i.le.inb(ij)))then |
174 |
lwork(ij)=(nent(ij,i).ne.0) |
175 |
qp1=q(ij,nk(ij))-ep(ij,i)*clw(ij,i) |
176 |
anum=h(ij,i)-hp(ij,i)-lv(ij,i)*(qp1-qs(ij,i)) |
177 |
denom=h(ij,i)-hp(ij,i)+lv(ij,i)*(q(ij,i)-qp1) |
178 |
if(abs(denom).lt.0.01)denom=0.01 |
179 |
scrit(ij)=anum/denom |
180 |
alt=qp1-qs(ij,i)+scrit(ij)*(q(ij,i)-qp1) |
181 |
if(scrit(ij).lt.0.0.or.alt.lt.0.0)scrit(ij)=1.0 |
182 |
asij(ij)=0.0 |
183 |
smin(ij)=1.0 |
184 |
endif |
185 |
781 continue |
186 |
do 783 j=minorig,nl |
187 |
! |
188 |
num2=0 |
189 |
do 778 ij=1,ncum |
190 |
if((i.ge.icb(ij)+1).and.(i.le.inb(ij)) & |
191 |
.and.(j.ge.icb(ij)).and.(j.le.inb(ij)) & |
192 |
.and.lwork(ij))num2=num2+1 |
193 |
778 continue |
194 |
if(num2.le.0)go to 783 |
195 |
! |
196 |
do 782 ij=1,ncum |
197 |
if((i.ge.icb(ij)+1).and.(i.le.inb(ij)) & |
198 |
.and.(j.ge.icb(ij)).and.(j.le.inb(ij)).and.lwork(ij))then |
199 |
if(sij(ij,i,j).gt.0.0.and.sij(ij,i,j).lt.0.9)then |
200 |
if(j.gt.i)then |
201 |
smid=min(sij(ij,i,j),scrit(ij)) |
202 |
sjmax=smid |
203 |
sjmin=smid |
204 |
if(smid.lt.smin(ij) & |
205 |
.and.sij(ij,i,j+1).lt.smid)then |
206 |
smin(ij)=smid |
207 |
sjmax=min(sij(ij,i,j+1),sij(ij,i,j),scrit(ij)) |
208 |
sjmin=max(sij(ij,i,j-1),sij(ij,i,j)) |
209 |
sjmin=min(sjmin,scrit(ij)) |
210 |
endif |
211 |
else |
212 |
sjmax=max(sij(ij,i,j+1),scrit(ij)) |
213 |
smid=max(sij(ij,i,j),scrit(ij)) |
214 |
sjmin=0.0 |
215 |
if(j.gt.1)sjmin=sij(ij,i,j-1) |
216 |
sjmin=max(sjmin,scrit(ij)) |
217 |
endif |
218 |
delp=abs(sjmax-smid) |
219 |
delm=abs(sjmin-smid) |
220 |
asij(ij)=asij(ij)+(delp+delm) & |
221 |
*(ph(ij,j)-ph(ij,j+1)) |
222 |
ment(ij,i,j)=ment(ij,i,j)*(delp+delm) & |
223 |
*(ph(ij,j)-ph(ij,j+1)) |
224 |
endif |
225 |
endif |
226 |
782 continue |
227 |
783 continue |
228 |
do 784 ij=1,ncum |
229 |
if((i.ge.icb(ij)+1).and.(i.le.inb(ij)).and.lwork(ij))then |
230 |
asij(ij)=max(1.0e-21,asij(ij)) |
231 |
asij(ij)=1.0/asij(ij) |
232 |
bsum(ij,i)=0.0 |
233 |
endif |
234 |
784 continue |
235 |
do 786 j=minorig,nl+1 |
236 |
do 785 ij=1,ncum |
237 |
if((i.ge.icb(ij)+1).and.(i.le.inb(ij)) & |
238 |
.and.(j.ge.icb(ij)).and.(j.le.inb(ij)) & |
239 |
.and.lwork(ij))then |
240 |
ment(ij,i,j)=ment(ij,i,j)*asij(ij) |
241 |
bsum(ij,i)=bsum(ij,i)+ment(ij,i,j) |
242 |
endif |
243 |
785 continue |
244 |
786 continue |
245 |
do 787 ij=1,ncum |
246 |
if((i.ge.icb(ij)+1).and.(i.le.inb(ij)) & |
247 |
.and.(bsum(ij,i).lt.1.0e-18).and.lwork(ij))then |
248 |
nent(ij,i)=0 |
249 |
ment(ij,i,i)=m(ij,i) |
250 |
qent(ij,i,i)=q(ij,nk(ij))-ep(ij,i)*clw(ij,i) |
251 |
uent(ij,i,i)=u(ij,nk(ij)) |
252 |
vent(ij,i,i)=v(ij,nk(ij)) |
253 |
elij(ij,i,i)=clw(ij,i) |
254 |
sij(ij,i,i)=1.0 |
255 |
endif |
256 |
787 continue |
257 |
789 continue |
258 |
! |
259 |
return |
260 |
end |