/[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 145 - (show annotations)
Tue Jun 16 15:23:29 2015 UTC (8 years, 11 months ago) by guez
Original Path: trunk/Sources/phylmd/CV3_routines/cv3_mixing.f
File size: 10885 byte(s)
Renamed bibio to misc.

In procedure fxhyp, use the fact that xf is an odd function of xtild.

In procedure invert_zoom_x, replace linear search in xf by
bisection. Also, use result from previous loop iteration as initial
guess. Variable "it" cannot be equal to 2 * nmax after search.

Unused arguments: hm of cv3_feed; ph, qnk, tv,tvp of cv3_mixing; ppsol
of lw; rconst, temp of vdif_kcay; rconst, plev, temp, ustar, l_mix of
yamada.

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

  ViewVC Help
Powered by ViewVC 1.1.21