/[lmdze]/trunk/libf/phylmd/CV3_routines/cv3_yield.f90
ViewVC logotype

Contents of /trunk/libf/phylmd/CV3_routines/cv3_yield.f90

Parent Directory Parent Directory | Revision Log Revision Log

Revision 47 - (show annotations)
Fri Jul 1 15:00:48 2011 UTC (12 years, 10 months ago) by guez
File size: 27971 byte(s) 
Split "thermcell.f" and "cv3_routines.f".
Removed copies of files that are now in "L_util".
Moved "mva9" and "diagetpq" to their own files.
Unified variable names across procedures.
1
2 SUBROUTINE cv3_yield(nloc,ncum,nd,na,ntra &
3 ,icb,inb,delt &
4 ,t,rr,u,v,tra,gz,p,ph,h,hp,lv,cpn,th &
5 ,ep,clw,m,tp,mp,rp,up,vp,trap &
6 ,wt,water,evap,b &
7 ,ment,qent,uent,vent,nent,elij,traent,sig &
8 ,tv,tvp &
9 ,iflag,precip,VPrecip,ft,fr,fu,fv,ftra &
10 ,upwd,dnwd,dnwd0,ma,mike,tls,tps,qcondc,wd)
11 use conema3_m
12 use cvparam3
13 use cvthermo
14 use cvflag
15 implicit none
16
17
18 ! inputs:
19 integer ncum,nd,na,ntra,nloc
20 integer icb(nloc), inb(nloc)
21 real, intent(in):: delt
22 real t(nloc,nd), rr(nloc,nd), u(nloc,nd), v(nloc,nd)
23 real tra(nloc,nd,ntra), sig(nloc,nd)
24 real gz(nloc,na), ph(nloc,nd+1), h(nloc,na), hp(nloc,na)
25 real th(nloc,na), p(nloc,nd), tp(nloc,na)
26 real lv(nloc,na), cpn(nloc,na), ep(nloc,na), clw(nloc,na)
27 real m(nloc,na), mp(nloc,na), rp(nloc,na), up(nloc,na)
28 real vp(nloc,na), wt(nloc,nd), trap(nloc,nd,ntra)
29 real water(nloc,na), evap(nloc,na), b(nloc,na)
30 real ment(nloc,na,na), qent(nloc,na,na), uent(nloc,na,na)
31 !ym real vent(nloc,na,na), nent(nloc,na), elij(nloc,na,na)
32 real vent(nloc,na,na), elij(nloc,na,na)
33 integer nent(nloc,na)
34 real traent(nloc,na,na,ntra)
35 real tv(nloc,nd), tvp(nloc,nd)
36
37 ! input/output:
38 integer iflag(nloc)
39
40 ! outputs:
41 real precip(nloc)
42 real VPrecip(nloc,nd+1)
43 real ft(nloc,nd), fr(nloc,nd), fu(nloc,nd), fv(nloc,nd)
44 real ftra(nloc,nd,ntra)
45 real upwd(nloc,nd), dnwd(nloc,nd), ma(nloc,nd)
46 real dnwd0(nloc,nd), mike(nloc,nd)
47 real tls(nloc,nd), tps(nloc,nd)
48 real qcondc(nloc,nd) ! cld
49 real wd(nloc) ! gust
50
51 ! local variables:
52 integer i,k,il,n,j,num1
53 real rat, awat, delti
54 real ax, bx, cx, dx, ex
55 real cpinv, rdcp, dpinv
56 real lvcp(nloc,na), mke(nloc,na)
57 real am(nloc), work(nloc), ad(nloc), amp1(nloc)
58 !!! real up1(nloc), dn1(nloc)
59 real up1(nloc,nd,nd), dn1(nloc,nd,nd)
60 real asum(nloc), bsum(nloc), csum(nloc), dsum(nloc)
61 real qcond(nloc,nd), nqcond(nloc,nd), wa(nloc,nd) ! cld
62 real siga(nloc,nd), sax(nloc,nd), mac(nloc,nd) ! cld
63
64
65 !-------------------------------------------------------------
66
67 ! initialization:
68
69 delti = 1.0/delt
70
71 do il=1,ncum
72 precip(il)=0.0
73 wd(il)=0.0 ! gust
74 VPrecip(il,nd+1)=0.
75 enddo
76
77 do i=1,nd
78 do il=1,ncum
79 VPrecip(il,i)=0.0
80 ft(il,i)=0.0
81 fr(il,i)=0.0
82 fu(il,i)=0.0
83 fv(il,i)=0.0
84 qcondc(il,i)=0.0 ! cld
85 qcond(il,i)=0.0 ! cld
86 nqcond(il,i)=0.0 ! cld
87 enddo
88 enddo
89
90 ! do j=1,ntra
91 ! do i=1,nd
92 ! do il=1,ncum
93 ! ftra(il,i,j)=0.0
94 ! enddo
95 ! enddo
96 ! enddo
97
98 do i=1,nl
99 do il=1,ncum
100 lvcp(il,i)=lv(il,i)/cpn(il,i)
101 enddo
102 enddo
103
104
105 !
106 ! *** calculate surface precipitation in mm/day ***
107 !
108 do il=1,ncum
109 if(ep(il,inb(il)).ge.0.0001)then
110 if (cvflag_grav) then
111 precip(il)=wt(il,1)*sigd*water(il,1)*86400.*1000./(rowl*grav)
112 else
113 precip(il)=wt(il,1)*sigd*water(il,1)*8640.
114 endif
115 endif
116 enddo
117
118 ! *** CALCULATE VERTICAL PROFILE OF PRECIPITATIONs IN kg/m2/s ===
119 !
120 ! MAF rajout pour lessivage
121 do k=1,nl
122 do il=1,ncum
123 if (k.le.inb(il)) then
124 if (cvflag_grav) then
125 VPrecip(il,k) = wt(il,k)*sigd*water(il,k)/grav
126 else
127 VPrecip(il,k) = wt(il,k)*sigd*water(il,k)/10.
128 endif
129 endif
130 end do
131 end do
132 !
133 !
134 ! *** Calculate downdraft velocity scale ***
135 ! *** NE PAS UTILISER POUR L'INSTANT ***
136 !
137 !! do il=1,ncum
138 !! wd(il)=betad*abs(mp(il,icb(il)))*0.01*rrd*t(il,icb(il))
139 !! : /(sigd*p(il,icb(il)))
140 !! enddo
141
142 !
143 ! *** calculate tendencies of lowest level potential temperature ***
144 ! *** and mixing ratio ***
145 !
146 do il=1,ncum
147 work(il)=1.0/(ph(il,1)-ph(il,2))
148 am(il)=0.0
149 enddo
150
151 do k=2,nl
152 do il=1,ncum
153 if (k.le.inb(il)) then
154 am(il)=am(il)+m(il,k)
155 endif
156 enddo
157 enddo
158
159 do il=1,ncum
160
161 ! convect3 if((0.1*dpinv*am).ge.delti)iflag(il)=4
162 if (cvflag_grav) then
163 if((0.01*grav*work(il)*am(il)).ge.delti)iflag(il)=1!consist vect
164 ft(il,1)=0.01*grav*work(il)*am(il)*(t(il,2)-t(il,1) &
165 +(gz(il,2)-gz(il,1))/cpn(il,1))
166 else
167 if((0.1*work(il)*am(il)).ge.delti)iflag(il)=1 !consistency vect
168 ft(il,1)=0.1*work(il)*am(il)*(t(il,2)-t(il,1) &
169 +(gz(il,2)-gz(il,1))/cpn(il,1))
170 endif
171
172 ft(il,1)=ft(il,1)-0.5*lvcp(il,1)*sigd*(evap(il,1)+evap(il,2))
173
174 if (cvflag_grav) then
175 ft(il,1)=ft(il,1)-0.009*grav*sigd*mp(il,2) &
176 *t(il,1)*b(il,1)*work(il)
177 else
178 ft(il,1)=ft(il,1)-0.09*sigd*mp(il,2)*t(il,1)*b(il,1)*work(il)
179 endif
180
181 ft(il,1)=ft(il,1)+0.01*sigd*wt(il,1)*(cl-cpd)*water(il,2)*(t(il,2) &
182 -t(il,1))*work(il)/cpn(il,1)
183
184 if (cvflag_grav) then
185 !jyg1 Correction pour mieux conserver l'eau (conformite avec CONVECT4.3)
186 ! (sb: pour l'instant, on ne fait que le chgt concernant grav, pas evap)
187 fr(il,1)=0.01*grav*mp(il,2)*(rp(il,2)-rr(il,1))*work(il) &
188 +sigd*0.5*(evap(il,1)+evap(il,2))
189 !+tard : +sigd*evap(il,1)
190
191 fr(il,1)=fr(il,1)+0.01*grav*am(il)*(rr(il,2)-rr(il,1))*work(il)
192
193 fu(il,1)=fu(il,1)+0.01*grav*work(il)*(mp(il,2)*(up(il,2)-u(il,1)) &
194 +am(il)*(u(il,2)-u(il,1)))
195 fv(il,1)=fv(il,1)+0.01*grav*work(il)*(mp(il,2)*(vp(il,2)-v(il,1)) &
196 +am(il)*(v(il,2)-v(il,1)))
197 else ! cvflag_grav
198 fr(il,1)=0.1*mp(il,2)*(rp(il,2)-rr(il,1))*work(il) &
199 +sigd*0.5*(evap(il,1)+evap(il,2))
200 fr(il,1)=fr(il,1)+0.1*am(il)*(rr(il,2)-rr(il,1))*work(il)
201 fu(il,1)=fu(il,1)+0.1*work(il)*(mp(il,2)*(up(il,2)-u(il,1)) &
202 +am(il)*(u(il,2)-u(il,1)))
203 fv(il,1)=fv(il,1)+0.1*work(il)*(mp(il,2)*(vp(il,2)-v(il,1)) &
204 +am(il)*(v(il,2)-v(il,1)))
205 endif ! cvflag_grav
206
207 enddo ! il
208
209 ! do j=1,ntra
210 ! do il=1,ncum
211 ! if (cvflag_grav) then
212 ! ftra(il,1,j)=ftra(il,1,j)+0.01*grav*work(il)
213 ! : *(mp(il,2)*(trap(il,2,j)-tra(il,1,j))
214 ! : +am(il)*(tra(il,2,j)-tra(il,1,j)))
215 ! else
216 ! ftra(il,1,j)=ftra(il,1,j)+0.1*work(il)
217 ! : *(mp(il,2)*(trap(il,2,j)-tra(il,1,j))
218 ! : +am(il)*(tra(il,2,j)-tra(il,1,j)))
219 ! endif
220 ! enddo
221 ! enddo
222
223 do j=2,nl
224 do il=1,ncum
225 if (j.le.inb(il)) then
226 if (cvflag_grav) then
227 fr(il,1)=fr(il,1) &
228 +0.01*grav*work(il)*ment(il,j,1)*(qent(il,j,1)-rr(il,1))
229 fu(il,1)=fu(il,1) &
230 +0.01*grav*work(il)*ment(il,j,1)*(uent(il,j,1)-u(il,1))
231 fv(il,1)=fv(il,1) &
232 +0.01*grav*work(il)*ment(il,j,1)*(vent(il,j,1)-v(il,1))
233 else ! cvflag_grav
234 fr(il,1)=fr(il,1) &
235 +0.1*work(il)*ment(il,j,1)*(qent(il,j,1)-rr(il,1))
236 fu(il,1)=fu(il,1) &
237 +0.1*work(il)*ment(il,j,1)*(uent(il,j,1)-u(il,1))
238 fv(il,1)=fv(il,1) &
239 +0.1*work(il)*ment(il,j,1)*(vent(il,j,1)-v(il,1))
240 endif ! cvflag_grav
241 endif ! j
242 enddo
243 enddo
244
245 ! do k=1,ntra
246 ! do j=2,nl
247 ! do il=1,ncum
248 ! if (j.le.inb(il)) then
249
250 ! if (cvflag_grav) then
251 ! ftra(il,1,k)=ftra(il,1,k)+0.01*grav*work(il)*ment(il,j,1)
252 ! : *(traent(il,j,1,k)-tra(il,1,k))
253 ! else
254 ! ftra(il,1,k)=ftra(il,1,k)+0.1*work(il)*ment(il,j,1)
255 ! : *(traent(il,j,1,k)-tra(il,1,k))
256 ! endif
257
258 ! endif
259 ! enddo
260 ! enddo
261 ! enddo
262
263 !
264 ! *** calculate tendencies of potential temperature and mixing ratio ***
265 ! *** at levels above the lowest level ***
266 !
267 ! *** first find the net saturated updraft and downdraft mass fluxes ***
268 ! *** through each level ***
269 !
270
271 do 500 i=2,nl+1 ! newvecto: mettre nl au lieu nl+1?
272
273 num1=0
274 do il=1,ncum
275 if(i.le.inb(il))num1=num1+1
276 enddo
277 if(num1.le.0)go to 500
278
279 call zilch(amp1,ncum)
280 call zilch(ad,ncum)
281
282 do 440 k=i+1,nl+1
283 do 441 il=1,ncum
284 if (i.le.inb(il) .and. k.le.(inb(il)+1)) then
285 amp1(il)=amp1(il)+m(il,k)
286 endif
287 441 continue
288 440 continue
289
290 do 450 k=1,i
291 do 451 j=i+1,nl+1
292 do 452 il=1,ncum
293 if (i.le.inb(il) .and. j.le.(inb(il)+1)) then
294 amp1(il)=amp1(il)+ment(il,k,j)
295 endif
296 452 continue
297 451 continue
298 450 continue
299
300 do 470 k=1,i-1
301 do 471 j=i,nl+1 ! newvecto: nl au lieu nl+1?
302 do 472 il=1,ncum
303 if (i.le.inb(il) .and. j.le.inb(il)) then
304 ad(il)=ad(il)+ment(il,j,k)
305 endif
306 472 continue
307 471 continue
308 470 continue
309
310 do 1350 il=1,ncum
311 if (i.le.inb(il)) then
312 dpinv=1.0/(ph(il,i)-ph(il,i+1))
313 cpinv=1.0/cpn(il,i)
314
315 ! convect3 if((0.1*dpinv*amp1).ge.delti)iflag(il)=4
316 if (cvflag_grav) then
317 if((0.01*grav*dpinv*amp1(il)).ge.delti)iflag(il)=1 ! vecto
318 else
319 if((0.1*dpinv*amp1(il)).ge.delti)iflag(il)=1 ! vecto
320 endif
321
322 if (cvflag_grav) then
323 ft(il,i)=0.01*grav*dpinv*(amp1(il)*(t(il,i+1)-t(il,i) &
324 +(gz(il,i+1)-gz(il,i))*cpinv) &
325 -ad(il)*(t(il,i)-t(il,i-1)+(gz(il,i)-gz(il,i-1))*cpinv)) &
326 -0.5*sigd*lvcp(il,i)*(evap(il,i)+evap(il,i+1))
327 rat=cpn(il,i-1)*cpinv
328 ft(il,i)=ft(il,i)-0.009*grav*sigd*(mp(il,i+1)*t(il,i)*b(il,i) &
329 -mp(il,i)*t(il,i-1)*rat*b(il,i-1))*dpinv
330 ft(il,i)=ft(il,i)+0.01*grav*dpinv*ment(il,i,i)*(hp(il,i)-h(il,i) &
331 +t(il,i)*(cpv-cpd)*(rr(il,i)-qent(il,i,i)))*cpinv
332 else ! cvflag_grav
333 ft(il,i)=0.1*dpinv*(amp1(il)*(t(il,i+1)-t(il,i) &
334 +(gz(il,i+1)-gz(il,i))*cpinv) &
335 -ad(il)*(t(il,i)-t(il,i-1)+(gz(il,i)-gz(il,i-1))*cpinv)) &
336 -0.5*sigd*lvcp(il,i)*(evap(il,i)+evap(il,i+1))
337 rat=cpn(il,i-1)*cpinv
338 ft(il,i)=ft(il,i)-0.09*sigd*(mp(il,i+1)*t(il,i)*b(il,i) &
339 -mp(il,i)*t(il,i-1)*rat*b(il,i-1))*dpinv
340 ft(il,i)=ft(il,i)+0.1*dpinv*ment(il,i,i)*(hp(il,i)-h(il,i) &
341 +t(il,i)*(cpv-cpd)*(rr(il,i)-qent(il,i,i)))*cpinv
342 endif ! cvflag_grav
343
344
345 ft(il,i)=ft(il,i)+0.01*sigd*wt(il,i)*(cl-cpd)*water(il,i+1) &
346 *(t(il,i+1)-t(il,i))*dpinv*cpinv
347
348 if (cvflag_grav) then
349 fr(il,i)=0.01*grav*dpinv*(amp1(il)*(rr(il,i+1)-rr(il,i)) &
350 -ad(il)*(rr(il,i)-rr(il,i-1)))
351 fu(il,i)=fu(il,i)+0.01*grav*dpinv*(amp1(il)*(u(il,i+1)-u(il,i)) &
352 -ad(il)*(u(il,i)-u(il,i-1)))
353 fv(il,i)=fv(il,i)+0.01*grav*dpinv*(amp1(il)*(v(il,i+1)-v(il,i)) &
354 -ad(il)*(v(il,i)-v(il,i-1)))
355 else ! cvflag_grav
356 fr(il,i)=0.1*dpinv*(amp1(il)*(rr(il,i+1)-rr(il,i)) &
357 -ad(il)*(rr(il,i)-rr(il,i-1)))
358 fu(il,i)=fu(il,i)+0.1*dpinv*(amp1(il)*(u(il,i+1)-u(il,i)) &
359 -ad(il)*(u(il,i)-u(il,i-1)))
360 fv(il,i)=fv(il,i)+0.1*dpinv*(amp1(il)*(v(il,i+1)-v(il,i)) &
361 -ad(il)*(v(il,i)-v(il,i-1)))
362 endif ! cvflag_grav
363
364 endif ! i
365 1350 continue
366
367 ! do k=1,ntra
368 ! do il=1,ncum
369 ! if (i.le.inb(il)) then
370 ! dpinv=1.0/(ph(il,i)-ph(il,i+1))
371 ! cpinv=1.0/cpn(il,i)
372 ! if (cvflag_grav) then
373 ! ftra(il,i,k)=ftra(il,i,k)+0.01*grav*dpinv
374 ! : *(amp1(il)*(tra(il,i+1,k)-tra(il,i,k))
375 ! : -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
376 ! else
377 ! ftra(il,i,k)=ftra(il,i,k)+0.1*dpinv
378 ! : *(amp1(il)*(tra(il,i+1,k)-tra(il,i,k))
379 ! : -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
380 ! endif
381 ! endif
382 ! enddo
383 ! enddo
384
385 do 480 k=1,i-1
386 do 1370 il=1,ncum
387 if (i.le.inb(il)) then
388 dpinv=1.0/(ph(il,i)-ph(il,i+1))
389 cpinv=1.0/cpn(il,i)
390
391 awat=elij(il,k,i)-(1.-ep(il,i))*clw(il,i)
392 awat=amax1(awat,0.0)
393
394 if (cvflag_grav) then
395 fr(il,i)=fr(il,i) &
396 +0.01*grav*dpinv*ment(il,k,i)*(qent(il,k,i)-awat-rr(il,i))
397 fu(il,i)=fu(il,i) &
398 +0.01*grav*dpinv*ment(il,k,i)*(uent(il,k,i)-u(il,i))
399 fv(il,i)=fv(il,i) &
400 +0.01*grav*dpinv*ment(il,k,i)*(vent(il,k,i)-v(il,i))
401 else ! cvflag_grav
402 fr(il,i)=fr(il,i) &
403 +0.1*dpinv*ment(il,k,i)*(qent(il,k,i)-awat-rr(il,i))
404 fu(il,i)=fu(il,i) &
405 +0.01*grav*dpinv*ment(il,k,i)*(uent(il,k,i)-u(il,i))
406 fv(il,i)=fv(il,i) &
407 +0.1*dpinv*ment(il,k,i)*(vent(il,k,i)-v(il,i))
408 endif ! cvflag_grav
409
410 ! (saturated updrafts resulting from mixing) ! cld
411 qcond(il,i)=qcond(il,i)+(elij(il,k,i)-awat) ! cld
412 nqcond(il,i)=nqcond(il,i)+1. ! cld
413 endif ! i
414 1370 continue
415 480 continue
416
417 ! do j=1,ntra
418 ! do k=1,i-1
419 ! do il=1,ncum
420 ! if (i.le.inb(il)) then
421 ! dpinv=1.0/(ph(il,i)-ph(il,i+1))
422 ! cpinv=1.0/cpn(il,i)
423 ! if (cvflag_grav) then
424 ! ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv*ment(il,k,i)
425 ! : *(traent(il,k,i,j)-tra(il,i,j))
426 ! else
427 ! ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv*ment(il,k,i)
428 ! : *(traent(il,k,i,j)-tra(il,i,j))
429 ! endif
430 ! endif
431 ! enddo
432 ! enddo
433 ! enddo
434
435 do 490 k=i,nl+1
436 do 1380 il=1,ncum
437 if (i.le.inb(il) .and. k.le.inb(il)) then
438 dpinv=1.0/(ph(il,i)-ph(il,i+1))
439 cpinv=1.0/cpn(il,i)
440
441 if (cvflag_grav) then
442 fr(il,i)=fr(il,i) &
443 +0.01*grav*dpinv*ment(il,k,i)*(qent(il,k,i)-rr(il,i))
444 fu(il,i)=fu(il,i) &
445 +0.01*grav*dpinv*ment(il,k,i)*(uent(il,k,i)-u(il,i))
446 fv(il,i)=fv(il,i) &
447 +0.01*grav*dpinv*ment(il,k,i)*(vent(il,k,i)-v(il,i))
448 else ! cvflag_grav
449 fr(il,i)=fr(il,i) &
450 +0.1*dpinv*ment(il,k,i)*(qent(il,k,i)-rr(il,i))
451 fu(il,i)=fu(il,i) &
452 +0.1*dpinv*ment(il,k,i)*(uent(il,k,i)-u(il,i))
453 fv(il,i)=fv(il,i) &
454 +0.1*dpinv*ment(il,k,i)*(vent(il,k,i)-v(il,i))
455 endif ! cvflag_grav
456 endif ! i and k
457 1380 continue
458 490 continue
459
460 ! do j=1,ntra
461 ! do k=i,nl+1
462 ! do il=1,ncum
463 ! if (i.le.inb(il) .and. k.le.inb(il)) then
464 ! dpinv=1.0/(ph(il,i)-ph(il,i+1))
465 ! cpinv=1.0/cpn(il,i)
466 ! if (cvflag_grav) then
467 ! ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv*ment(il,k,i)
468 ! : *(traent(il,k,i,j)-tra(il,i,j))
469 ! else
470 ! ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv*ment(il,k,i)
471 ! : *(traent(il,k,i,j)-tra(il,i,j))
472 ! endif
473 ! endif ! i and k
474 ! enddo
475 ! enddo
476 ! enddo
477
478 do 1400 il=1,ncum
479 if (i.le.inb(il)) then
480 dpinv=1.0/(ph(il,i)-ph(il,i+1))
481 cpinv=1.0/cpn(il,i)
482
483 if (cvflag_grav) then
484 ! sb: on ne fait pas encore la correction permettant de mieux
485 ! conserver l'eau:
486 fr(il,i)=fr(il,i)+0.5*sigd*(evap(il,i)+evap(il,i+1)) &
487 +0.01*grav*(mp(il,i+1)*(rp(il,i+1)-rr(il,i))-mp(il,i) &
488 *(rp(il,i)-rr(il,i-1)))*dpinv
489
490 fu(il,i)=fu(il,i)+0.01*grav*(mp(il,i+1)*(up(il,i+1)-u(il,i)) &
491 -mp(il,i)*(up(il,i)-u(il,i-1)))*dpinv
492 fv(il,i)=fv(il,i)+0.01*grav*(mp(il,i+1)*(vp(il,i+1)-v(il,i)) &
493 -mp(il,i)*(vp(il,i)-v(il,i-1)))*dpinv
494 else ! cvflag_grav
495 fr(il,i)=fr(il,i)+0.5*sigd*(evap(il,i)+evap(il,i+1)) &
496 +0.1*(mp(il,i+1)*(rp(il,i+1)-rr(il,i))-mp(il,i) &
497 *(rp(il,i)-rr(il,i-1)))*dpinv
498 fu(il,i)=fu(il,i)+0.1*(mp(il,i+1)*(up(il,i+1)-u(il,i)) &
499 -mp(il,i)*(up(il,i)-u(il,i-1)))*dpinv
500 fv(il,i)=fv(il,i)+0.1*(mp(il,i+1)*(vp(il,i+1)-v(il,i)) &
501 -mp(il,i)*(vp(il,i)-v(il,i-1)))*dpinv
502 endif ! cvflag_grav
503
504 endif ! i
505 1400 continue
506
507 ! sb: interface with the cloud parameterization: ! cld
508
509 do k=i+1,nl
510 do il=1,ncum
511 if (k.le.inb(il) .and. i.le.inb(il)) then ! cld
512 ! (saturated downdrafts resulting from mixing) ! cld
513 qcond(il,i)=qcond(il,i)+elij(il,k,i) ! cld
514 nqcond(il,i)=nqcond(il,i)+1. ! cld
515 endif ! cld
516 enddo ! cld
517 enddo ! cld
518
519 ! (particular case: no detraining level is found) ! cld
520 do il=1,ncum ! cld
521 if (i.le.inb(il) .and. nent(il,i).eq.0) then ! cld
522 qcond(il,i)=qcond(il,i)+(1.-ep(il,i))*clw(il,i) ! cld
523 nqcond(il,i)=nqcond(il,i)+1. ! cld
524 endif ! cld
525 enddo ! cld
526
527 do il=1,ncum ! cld
528 if (i.le.inb(il) .and. nqcond(il,i).ne.0.) then ! cld
529 qcond(il,i)=qcond(il,i)/nqcond(il,i) ! cld
530 endif ! cld
531 enddo
532
533 ! do j=1,ntra
534 ! do il=1,ncum
535 ! if (i.le.inb(il)) then
536 ! dpinv=1.0/(ph(il,i)-ph(il,i+1))
537 ! cpinv=1.0/cpn(il,i)
538
539 ! if (cvflag_grav) then
540 ! ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv
541 ! : *(mp(il,i+1)*(trap(il,i+1,j)-tra(il,i,j))
542 ! : -mp(il,i)*(trap(il,i,j)-tra(il,i-1,j)))
543 ! else
544 ! ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv
545 ! : *(mp(il,i+1)*(trap(il,i+1,j)-tra(il,i,j))
546 ! : -mp(il,i)*(trap(il,i,j)-tra(il,i-1,j)))
547 ! endif
548 ! endif ! i
549 ! enddo
550 ! enddo
551
552 500 continue
553
554
555 ! *** move the detrainment at level inb down to level inb-1 ***
556 ! *** in such a way as to preserve the vertically ***
557 ! *** integrated enthalpy and water tendencies ***
558 !
559 do 503 il=1,ncum
560
561 ax=0.1*ment(il,inb(il),inb(il))*(hp(il,inb(il))-h(il,inb(il)) &
562 +t(il,inb(il))*(cpv-cpd) &
563 *(rr(il,inb(il))-qent(il,inb(il),inb(il)))) &
564 /(cpn(il,inb(il))*(ph(il,inb(il))-ph(il,inb(il)+1)))
565 ft(il,inb(il))=ft(il,inb(il))-ax
566 ft(il,inb(il)-1)=ft(il,inb(il)-1)+ax*cpn(il,inb(il)) &
567 *(ph(il,inb(il))-ph(il,inb(il)+1))/(cpn(il,inb(il)-1) &
568 *(ph(il,inb(il)-1)-ph(il,inb(il))))
569
570 bx=0.1*ment(il,inb(il),inb(il))*(qent(il,inb(il),inb(il)) &
571 -rr(il,inb(il)))/(ph(il,inb(il))-ph(il,inb(il)+1))
572 fr(il,inb(il))=fr(il,inb(il))-bx
573 fr(il,inb(il)-1)=fr(il,inb(il)-1) &
574 +bx*(ph(il,inb(il))-ph(il,inb(il)+1)) &
575 /(ph(il,inb(il)-1)-ph(il,inb(il)))
576
577 cx=0.1*ment(il,inb(il),inb(il))*(uent(il,inb(il),inb(il)) &
578 -u(il,inb(il)))/(ph(il,inb(il))-ph(il,inb(il)+1))
579 fu(il,inb(il))=fu(il,inb(il))-cx
580 fu(il,inb(il)-1)=fu(il,inb(il)-1) &
581 +cx*(ph(il,inb(il))-ph(il,inb(il)+1)) &
582 /(ph(il,inb(il)-1)-ph(il,inb(il)))
583
584 dx=0.1*ment(il,inb(il),inb(il))*(vent(il,inb(il),inb(il)) &
585 -v(il,inb(il)))/(ph(il,inb(il))-ph(il,inb(il)+1))
586 fv(il,inb(il))=fv(il,inb(il))-dx
587 fv(il,inb(il)-1)=fv(il,inb(il)-1) &
588 +dx*(ph(il,inb(il))-ph(il,inb(il)+1)) &
589 /(ph(il,inb(il)-1)-ph(il,inb(il)))
590
591 503 continue
592
593 ! do j=1,ntra
594 ! do il=1,ncum
595 ! ex=0.1*ment(il,inb(il),inb(il))
596 ! : *(traent(il,inb(il),inb(il),j)-tra(il,inb(il),j))
597 ! : /(ph(il,inb(il))-ph(il,inb(il)+1))
598 ! ftra(il,inb(il),j)=ftra(il,inb(il),j)-ex
599 ! ftra(il,inb(il)-1,j)=ftra(il,inb(il)-1,j)
600 ! : +ex*(ph(il,inb(il))-ph(il,inb(il)+1))
601 ! : /(ph(il,inb(il)-1)-ph(il,inb(il)))
602 ! enddo
603 ! enddo
604
605 !
606 ! *** homoginize tendencies below cloud base ***
607 !
608 !
609 do il=1,ncum
610 asum(il)=0.0
611 bsum(il)=0.0
612 csum(il)=0.0
613 dsum(il)=0.0
614 enddo
615
616 do i=1,nl
617 do il=1,ncum
618 if (i.le.(icb(il)-1)) then
619 asum(il)=asum(il)+ft(il,i)*(ph(il,i)-ph(il,i+1))
620 bsum(il)=bsum(il)+fr(il,i)*(lv(il,i)+(cl-cpd)*(t(il,i)-t(il,1))) &
621 *(ph(il,i)-ph(il,i+1))
622 csum(il)=csum(il)+(lv(il,i)+(cl-cpd)*(t(il,i)-t(il,1))) &
623 *(ph(il,i)-ph(il,i+1))
624 dsum(il)=dsum(il)+t(il,i)*(ph(il,i)-ph(il,i+1))/th(il,i)
625 endif
626 enddo
627 enddo
628
629 !!!! do 700 i=1,icb(il)-1
630 do i=1,nl
631 do il=1,ncum
632 if (i.le.(icb(il)-1)) then
633 ft(il,i)=asum(il)*t(il,i)/(th(il,i)*dsum(il))
634 fr(il,i)=bsum(il)/csum(il)
635 endif
636 enddo
637 enddo
638
639 !
640 ! *** reset counter and return ***
641 !
642 do il=1,ncum
643 sig(il,nd)=2.0
644 enddo
645
646
647 do i=1,nd
648 do il=1,ncum
649 upwd(il,i)=0.0
650 dnwd(il,i)=0.0
651 enddo
652 enddo
653
654 do i=1,nl
655 do il=1,ncum
656 dnwd0(il,i)=-mp(il,i)
657 enddo
658 enddo
659 do i=nl+1,nd
660 do il=1,ncum
661 dnwd0(il,i)=0.
662 enddo
663 enddo
664
665
666 do i=1,nl
667 do il=1,ncum
668 if (i.ge.icb(il) .and. i.le.inb(il)) then
669 upwd(il,i)=0.0
670 dnwd(il,i)=0.0
671 endif
672 enddo
673 enddo
674
675 do i=1,nl
676 do k=1,nl
677 do il=1,ncum
678 up1(il,k,i)=0.0
679 dn1(il,k,i)=0.0
680 enddo
681 enddo
682 enddo
683
684 do i=1,nl
685 do k=i,nl
686 do n=1,i-1
687 do il=1,ncum
688 if (i.ge.icb(il).and.i.le.inb(il).and.k.le.inb(il)) then
689 up1(il,k,i)=up1(il,k,i)+ment(il,n,k)
690 dn1(il,k,i)=dn1(il,k,i)-ment(il,k,n)
691 endif
692 enddo
693 enddo
694 enddo
695 enddo
696
697 do i=2,nl
698 do k=i,nl
699 do il=1,ncum
700 !test if (i.ge.icb(il).and.i.le.inb(il).and.k.le.inb(il)) then
701 if (i.le.inb(il).and.k.le.inb(il)) then
702 upwd(il,i)=upwd(il,i)+m(il,k)+up1(il,k,i)
703 dnwd(il,i)=dnwd(il,i)+dn1(il,k,i)
704 endif
705 enddo
706 enddo
707 enddo
708
709
710 !!!! DO il=1,ncum
711 !!!! do i=icb(il),inb(il)
712 !!!!
713 !!!! upwd(il,i)=0.0
714 !!!! dnwd(il,i)=0.0
715 !!!! do k=i,inb(il)
716 !!!! up1=0.0
717 !!!! dn1=0.0
718 !!!! do n=1,i-1
719 !!!! up1=up1+ment(il,n,k)
720 !!!! dn1=dn1-ment(il,k,n)
721 !!!! enddo
722 !!!! upwd(il,i)=upwd(il,i)+m(il,k)+up1
723 !!!! dnwd(il,i)=dnwd(il,i)+dn1
724 !!!! enddo
725 !!!! enddo
726 !!!!
727 !!!! ENDDO
728
729 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
730 ! determination de la variation de flux ascendant entre
731 ! deux niveau non dilue mike
732 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
733
734 do i=1,nl
735 do il=1,ncum
736 mike(il,i)=m(il,i)
737 enddo
738 enddo
739
740 do i=nl+1,nd
741 do il=1,ncum
742 mike(il,i)=0.
743 enddo
744 enddo
745
746 do i=1,nd
747 do il=1,ncum
748 ma(il,i)=0
749 enddo
750 enddo
751
752 do i=1,nl
753 do j=i,nl
754 do il=1,ncum
755 ma(il,i)=ma(il,i)+m(il,j)
756 enddo
757 enddo
758 enddo
759
760 do i=nl+1,nd
761 do il=1,ncum
762 ma(il,i)=0.
763 enddo
764 enddo
765
766 do i=1,nl
767 do il=1,ncum
768 if (i.le.(icb(il)-1)) then
769 ma(il,i)=0
770 endif
771 enddo
772 enddo
773
774 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
775 ! icb represente de niveau ou se trouve la
776 ! base du nuage , et inb le top du nuage
777 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
778
779 do i=1,nd
780 do il=1,ncum
781 mke(il,i)=upwd(il,i)+dnwd(il,i)
782 enddo
783 enddo
784
785 do i=1,nd
786 DO 999 il=1,ncum
787 rdcp=(rrd*(1.-rr(il,i))-rr(il,i)*rrv) &
788 /(cpd*(1.-rr(il,i))+rr(il,i)*cpv)
789 tls(il,i)=t(il,i)*(1000.0/p(il,i))**rdcp
790 tps(il,i)=tp(il,i)
791 999 CONTINUE
792 enddo
793
794 !
795 ! *** diagnose the in-cloud mixing ratio *** ! cld
796 ! *** of condensed water *** ! cld
797 ! ! cld
798
799 do i=1,nd ! cld
800 do il=1,ncum ! cld
801 mac(il,i)=0.0 ! cld
802 wa(il,i)=0.0 ! cld
803 siga(il,i)=0.0 ! cld
804 sax(il,i)=0.0 ! cld
805 enddo ! cld
806 enddo ! cld
807
808 do i=minorig, nl ! cld
809 do k=i+1,nl+1 ! cld
810 do il=1,ncum ! cld
811 if (i.le.inb(il) .and. k.le.(inb(il)+1)) then ! cld
812 mac(il,i)=mac(il,i)+m(il,k) ! cld
813 endif ! cld
814 enddo ! cld
815 enddo ! cld
816 enddo ! cld
817
818 do i=1,nl ! cld
819 do j=1,i ! cld
820 do il=1,ncum ! cld
821 if (i.ge.icb(il) .and. i.le.(inb(il)-1) &
822 .and. j.ge.icb(il) ) then ! cld
823 sax(il,i)=sax(il,i)+rrd*(tvp(il,j)-tv(il,j)) &
824 *(ph(il,j)-ph(il,j+1))/p(il,j) ! cld
825 endif ! cld
826 enddo ! cld
827 enddo ! cld
828 enddo ! cld
829
830 do i=1,nl ! cld
831 do il=1,ncum ! cld
832 if (i.ge.icb(il) .and. i.le.(inb(il)-1) &
833 .and. sax(il,i).gt.0.0 ) then ! cld
834 wa(il,i)=sqrt(2.*sax(il,i)) ! cld
835 endif ! cld
836 enddo ! cld
837 enddo ! cld
838
839 do i=1,nl ! cld
840 do il=1,ncum ! cld
841 if (wa(il,i).gt.0.0) &
842 siga(il,i)=mac(il,i)/wa(il,i) &
843 *rrd*tvp(il,i)/p(il,i)/100./delta ! cld
844 siga(il,i) = min(siga(il,i),1.0) ! cld
845 !IM cf. FH
846 if (iflag_clw.eq.0) then
847 qcondc(il,i)=siga(il,i)*clw(il,i)*(1.-ep(il,i)) &
848 + (1.-siga(il,i))*qcond(il,i) ! cld
849 else if (iflag_clw.eq.1) then
850 qcondc(il,i)=qcond(il,i) ! cld
851 endif
852
853 enddo ! cld
854 enddo ! cld
855
856 return
857 end
  ViewVC Help
Powered by ViewVC 1.1.21