/[lmdze]/trunk/Sources/phylmd/CV30_routines/cv30_mixing.f
ViewVC logotype

Contents of /trunk/Sources/phylmd/CV30_routines/cv30_mixing.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 134 - (show annotations)
Wed Apr 29 15:47:56 2015 UTC (9 years ago) by guez
Original Path: trunk/Sources/phylmd/CV3_routines/cv3_mixing.f
File size: 11234 byte(s)
Sources inside, compilation outside.
1 module cv3_mixing_m
2
3 implicit none
4
5 contains
6
7 SUBROUTINE cv3_mixing(nloc,ncum,nd,na,icb,nk,inb &
8 ,ph,t,rr,rs,u,v,h,lv,qnk &
9 ,hp,tv,tvp,ep,clw,m,sig &
10 ,ment,qent,uent,vent, nent, sij,elij,ments,qents)
11 use cv3_param_m
12 use cvthermo
13
14 !---------------------------------------------------------------------
15 ! a faire:
16 ! - changer rr(il,1) -> qnk(il)
17 ! - vectorisation de la partie normalisation des flux (do 789...)
18 !---------------------------------------------------------------------
19
20
21 ! inputs:
22 integer, intent(in):: ncum, nd, na, nloc
23 integer icb(nloc), inb(nloc), nk(nloc)
24 real sig(nloc,nd)
25 real qnk(nloc)
26 real ph(nloc,nd+1)
27 real t(nloc,nd), rr(nloc,nd), rs(nloc,nd)
28 real u(nloc,nd), v(nloc,nd)
29 real lv(nloc,na), h(nloc,na), hp(nloc,na)
30 real tv(nloc,na), tvp(nloc,na), ep(nloc,na), clw(nloc,na)
31 real m(nloc,na) ! input of convect3
32
33 ! outputs:
34 real ment(nloc,na,na), qent(nloc,na,na)
35 real uent(nloc,na,na), vent(nloc,na,na)
36 real sij(nloc,na,na), elij(nloc,na,na)
37 real ments(nloc,nd,nd), qents(nloc,nd,nd)
38 real sigij(nloc,nd,nd)
39 integer nent(nloc,nd)
40
41 ! local variables:
42 integer i, j, k, il, im, jm
43 integer num1, num2
44 real rti, bf2, anum, denom, dei, altem, cwat, stemp, qp
45 real alt, smid, sjmin, sjmax, delp, delm
46 real asij(nloc), smax(nloc), scrit(nloc)
47 real asum(nloc,nd),bsum(nloc,nd),csum(nloc,nd)
48 real wgh
49 real zm(nloc,na)
50 logical lwork(nloc)
51
52 !=====================================================================
53 ! --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS
54 !=====================================================================
55
56 do j=1,nl
57 do i=1,ncum
58 nent(i,j)=0
59 ! in convect3, m is computed in cv3_closure
60 ! ori m(i,1)=0.0
61 end do
62 end do
63
64 do j=1,nl
65 do k=1,nl
66 do i=1,ncum
67 qent(i,k,j)=rr(i,j)
68 uent(i,k,j)=u(i,j)
69 vent(i,k,j)=v(i,j)
70 elij(i,k,j)=0.0
71 !ym ment(i,k,j)=0.0
72 !ym sij(i,k,j)=0.0
73 end do
74 end do
75 end do
76
77 !ym
78 ment(1:ncum,1:nd,1:nd)=0.0
79 sij(1:ncum,1:nd,1:nd)=0.0
80
81 zm(:,:)=0.
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 do i=minorig+1, nl
90
91 do j=minorig,nl
92 do il=1,ncum
93 if( (i.ge.icb(il)).and.(i.le.inb(il)).and. &
94 (j.ge.(icb(il)-1)).and.(j.le.inb(il)))then
95
96 rti=rr(il,1)-ep(il,i)*clw(il,i)
97 bf2=1.+lv(il,j)*lv(il,j)*rs(il,j)/(rrv*t(il,j)*t(il,j)*cpd)
98 anum=h(il,j)-hp(il,i)+(cpv-cpd)*t(il,j)*(rti-rr(il,j))
99 denom=h(il,i)-hp(il,i)+(cpd-cpv)*(rr(il,i)-rti)*t(il,j)
100 dei=denom
101 if(abs(dei).lt.0.01)dei=0.01
102 sij(il,i,j)=anum/dei
103 sij(il,i,i)=1.0
104 altem=sij(il,i,j)*rr(il,i)+(1.-sij(il,i,j))*rti-rs(il,j)
105 altem=altem/bf2
106 cwat=clw(il,j)*(1.-ep(il,j))
107 stemp=sij(il,i,j)
108 if((stemp.lt.0.0.or.stemp.gt.1.0.or.altem.gt.cwat) &
109 .and.j.gt.i)then
110 anum=anum-lv(il,j)*(rti-rs(il,j)-cwat*bf2)
111 denom=denom+lv(il,j)*(rr(il,i)-rti)
112 if(abs(denom).lt.0.01)denom=0.01
113 sij(il,i,j)=anum/denom
114 altem=sij(il,i,j)*rr(il,i)+(1.-sij(il,i,j))*rti-rs(il,j)
115 altem=altem-(bf2-1.)*cwat
116 end if
117 if(sij(il,i,j).gt.0.0.and.sij(il,i,j).lt.0.95)then
118 qent(il,i,j)=sij(il,i,j)*rr(il,i)+(1.-sij(il,i,j))*rti
119 uent(il,i,j)=sij(il,i,j)*u(il,i)+(1.-sij(il,i,j))*u(il,nk(il))
120 vent(il,i,j)=sij(il,i,j)*v(il,i)+(1.-sij(il,i,j))*v(il,nk(il))
121 elij(il,i,j)=altem
122 elij(il,i,j)=amax1(0.0,elij(il,i,j))
123 ment(il,i,j)=m(il,i)/(1.-sij(il,i,j))
124 nent(il,i)=nent(il,i)+1
125 end if
126 sij(il,i,j)=amax1(0.0,sij(il,i,j))
127 sij(il,i,j)=amin1(1.0,sij(il,i,j))
128 endif ! new
129 end do
130 end do
131
132 !
133 ! *** if no air can entrain at level i assume that updraft detrains ***
134 ! *** at that level and calculate detrained air flux and properties ***
135 !
136
137 !@ do 170 i=icb(il),inb(il)
138
139 do il=1,ncum
140 if ((i.ge.icb(il)).and.(i.le.inb(il)).and.(nent(il,i).eq.0)) then
141 !@ if(nent(il,i).eq.0)then
142 ment(il,i,i)=m(il,i)
143 qent(il,i,i)=rr(il,nk(il))-ep(il,i)*clw(il,i)
144 uent(il,i,i)=u(il,nk(il))
145 vent(il,i,i)=v(il,nk(il))
146 elij(il,i,i)=clw(il,i)
147 !MAF sij(il,i,i)=1.0
148 sij(il,i,i)=0.0
149 end if
150 end do
151 end do
152
153 do j=minorig,nl
154 do i=minorig,nl
155 do il=1,ncum
156 if ((j.ge.(icb(il)-1)).and.(j.le.inb(il)) &
157 .and.(i.ge.icb(il)).and.(i.le.inb(il)))then
158 sigij(il,i,j)=sij(il,i,j)
159 endif
160 end do
161 end do
162 end do
163 !@ enddo
164
165 !@170 continue
166
167 !=====================================================================
168 ! --- NORMALIZE ENTRAINED AIR MASS FLUXES
169 ! --- TO REPRESENT EQUAL PROBABILITIES OF MIXING
170 !=====================================================================
171
172 call zilch(asum,nloc*nd)
173 call zilch(csum,nloc*nd)
174 call zilch(csum,nloc*nd)
175
176 do il=1,ncum
177 lwork(il) = .FALSE.
178 enddo
179
180 DO i=minorig+1,nl
181
182 num1=0
183 do il=1,ncum
184 if ( i.ge.icb(il) .and. i.le.inb(il) ) num1=num1+1
185 enddo
186 if (num1.le.0) cycle
187
188
189 do il=1,ncum
190 if ( i.ge.icb(il) .and. i.le.inb(il) ) then
191 lwork(il)=(nent(il,i).ne.0)
192 qp=rr(il,1)-ep(il,i)*clw(il,i)
193 anum=h(il,i)-hp(il,i)-lv(il,i)*(qp-rs(il,i)) &
194 +(cpv-cpd)*t(il,i)*(qp-rr(il,i))
195 denom=h(il,i)-hp(il,i)+lv(il,i)*(rr(il,i)-qp) &
196 +(cpd-cpv)*t(il,i)*(rr(il,i)-qp)
197 if(abs(denom).lt.0.01)denom=0.01
198 scrit(il)=anum/denom
199 alt=qp-rs(il,i)+scrit(il)*(rr(il,i)-qp)
200 if(scrit(il).le.0.0.or.alt.le.0.0)scrit(il)=1.0
201 smax(il)=0.0
202 asij(il)=0.0
203 endif
204 end do
205
206 do j=nl,minorig,-1
207
208 num2=0
209 do il=1,ncum
210 if ( i.ge.icb(il) .and. i.le.inb(il) .and. &
211 j.ge.(icb(il)-1) .and. j.le.inb(il) &
212 .and. lwork(il) ) num2=num2+1
213 enddo
214 if (num2.le.0) cycle
215
216 do il=1,ncum
217 if ( i.ge.icb(il) .and. i.le.inb(il) .and. &
218 j.ge.(icb(il)-1) .and. j.le.inb(il) &
219 .and. lwork(il) ) then
220
221 if(sij(il,i,j).gt.1.0e-16.and.sij(il,i,j).lt.0.95)then
222 wgh=1.0
223 if(j.gt.i)then
224 sjmax=amax1(sij(il,i,j+1),smax(il))
225 sjmax=amin1(sjmax,scrit(il))
226 smax(il)=amax1(sij(il,i,j),smax(il))
227 sjmin=amax1(sij(il,i,j-1),smax(il))
228 sjmin=amin1(sjmin,scrit(il))
229 if(sij(il,i,j).lt.(smax(il)-1.0e-16))wgh=0.0
230 smid=amin1(sij(il,i,j),scrit(il))
231 else
232 sjmax=amax1(sij(il,i,j+1),scrit(il))
233 smid=amax1(sij(il,i,j),scrit(il))
234 sjmin=0.0
235 if(j.gt.1)sjmin=sij(il,i,j-1)
236 sjmin=amax1(sjmin,scrit(il))
237 endif
238 delp=abs(sjmax-smid)
239 delm=abs(sjmin-smid)
240 asij(il)=asij(il)+wgh*(delp+delm)
241 ment(il,i,j)=ment(il,i,j)*(delp+delm)*wgh
242 endif
243 endif
244 end do
245
246 end do
247
248 do il=1,ncum
249 if (i.ge.icb(il).and.i.le.inb(il).and.lwork(il)) then
250 asij(il)=amax1(1.0e-16,asij(il))
251 asij(il)=1.0/asij(il)
252 asum(il,i)=0.0
253 bsum(il,i)=0.0
254 csum(il,i)=0.0
255 endif
256 enddo
257
258 do j=minorig,nl
259 do il=1,ncum
260 if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il) &
261 .and. j.ge.(icb(il)-1) .and. j.le.inb(il) ) then
262 ment(il,i,j)=ment(il,i,j)*asij(il)
263 endif
264 enddo
265 end do
266
267 do j=minorig,nl
268 do il=1,ncum
269 if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il) &
270 .and. j.ge.(icb(il)-1) .and. j.le.inb(il) ) then
271 asum(il,i)=asum(il,i)+ment(il,i,j)
272 ment(il,i,j)=ment(il,i,j)*sig(il,j)
273 bsum(il,i)=bsum(il,i)+ment(il,i,j)
274 endif
275 enddo
276 end do
277
278 do il=1,ncum
279 if (i.ge.icb(il).and.i.le.inb(il).and.lwork(il)) then
280 bsum(il,i)=amax1(bsum(il,i),1.0e-16)
281 bsum(il,i)=1.0/bsum(il,i)
282 endif
283 enddo
284
285 do j=minorig,nl
286 do il=1,ncum
287 if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il) &
288 .and. j.ge.(icb(il)-1) .and. j.le.inb(il) ) then
289 ment(il,i,j)=ment(il,i,j)*asum(il,i)*bsum(il,i)
290 endif
291 enddo
292 end do
293
294 do j=minorig,nl
295 do il=1,ncum
296 if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il) &
297 .and. j.ge.(icb(il)-1) .and. j.le.inb(il) ) then
298 csum(il,i)=csum(il,i)+ment(il,i,j)
299 endif
300 enddo
301 end do
302
303 do il=1,ncum
304 if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il) &
305 .and. csum(il,i).lt.m(il,i) ) then
306 nent(il,i)=0
307 ment(il,i,i)=m(il,i)
308 qent(il,i,i)=rr(il,1)-ep(il,i)*clw(il,i)
309 uent(il,i,i)=u(il,nk(il))
310 vent(il,i,i)=v(il,nk(il))
311 elij(il,i,i)=clw(il,i)
312 !MAF sij(il,i,i)=1.0
313 sij(il,i,i)=0.0
314 endif
315 enddo ! il
316
317 end DO
318 !
319 ! MAF: renormalisation de MENT
320 do jm=1,nd
321 do im=1,nd
322 do il=1,ncum
323 zm(il,im)=zm(il,im)+(1.-sij(il,im,jm))*ment(il,im,jm)
324 end do
325 end do
326 end do
327 !
328 do jm=1,nd
329 do im=1,nd
330 do il=1,ncum
331 if(zm(il,im).ne.0.) then
332 ment(il,im,jm)=ment(il,im,jm)*m(il,im)/zm(il,im)
333 endif
334 end do
335 end do
336 end do
337 !
338 do jm=1,nd
339 do im=1,nd
340 do il=1,ncum
341 qents(il,im,jm)=qent(il,im,jm)
342 ments(il,im,jm)=ment(il,im,jm)
343 end do
344 enddo
345 enddo
346
347 end SUBROUTINE cv3_mixing
348
349 end module cv3_mixing_m

  ViewVC Help
Powered by ViewVC 1.1.21