/[lmdze]/trunk/libf/phylmd/CV_routines/cv_mixing.f90
ViewVC logotype

Contents of /trunk/libf/phylmd/CV_routines/cv_mixing.f90

Parent Directory Parent Directory | Revision Log Revision Log


Revision 52 - (show annotations)
Fri Sep 23 12:28:01 2011 UTC (12 years, 7 months ago) by guez
File size: 9393 byte(s)
Split "conflx.f" into single-procedure files in directory "Conflx".

Split "cv_routines.f" into single-procedure files in directory
"CV_routines". Made module "cvparam" from included file
"cvparam.h". No included file other than "netcdf.inc" left in LMDZE.

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

  ViewVC Help
Powered by ViewVC 1.1.21