/[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 186 - (show annotations)
Mon Mar 21 15:36:26 2016 UTC (8 years, 1 month ago) by guez
File size: 10471 byte(s)
Removed variables nlm and nlp of module cv30_param_m. We do not
believe much in the benefit of these intermediary variables so we go
for clarity.

Removed variable noff of module cv30_param_m. Never used anywhere
else. Just set the value of nl explicitly in cv30_param.

Removed argument nd of cv30_param. Only called with nd = klev.

Replaced calls to zilch by array assignments. There was a strange
double call to zilch with the same arguments in cv30_mixing.

Removed procedure cv_flag. Just set the value of variable cvflag_grav
of module cvflag at declaration.

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

  ViewVC Help
Powered by ViewVC 1.1.21