/[lmdze]/trunk/phylmd/CV_routines/cv_mixing.f
ViewVC logotype

Contents of /trunk/phylmd/CV_routines/cv_mixing.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 97 - (show annotations)
Fri Apr 25 14:58:31 2014 UTC (10 years ago) by guez
File size: 9407 byte(s)
Module pressure_var is now only used in gcm. Created local variables
pls and p3d in etat0, added argument p3d to regr_pr_o3.

In leapfrog, moved computation of p3d and exner function immediately
after integrd, for clarity (does not change the execution).

Removed unused arguments: ntra, tra1 and tra of cv3_compress; ntra,
tra and traent of cv3_mixing; ntra, ftra, ftra1 of cv3_uncompress;
ntra, tra, trap of cv3_unsat; ntra, tra, trap, traent, ftra of
cv3_yield; tra, tvp, pbase, bbase, dtvpdt1, dtvpdq1, dplcldt,
dplcldr, ntra of concvl; ndp1, ntra, tra1 of cv_driver

Removed argument d_tra and computation of d_tra in concvl. Removed
argument ftra1 and computation of ftra1 in cv_driver. ftra1 was just
set to 0 in cv_driver, associated to d_tra in concvl, and set again to
zero in concvl.

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, intent(in):: 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

  ViewVC Help
Powered by ViewVC 1.1.21