My Project
 All Classes Files Functions Variables Macros
cv30_routines.F
Go to the documentation of this file.
1 !
2 ! $Id: cv30_routines.F 1744 2013-04-09 06:31:40Z emillour $
3 !
4 c
5 c
6  SUBROUTINE cv30_param(nd,delt)
7  implicit none
8 
9 c------------------------------------------------------------
10 c Set parameters for convectL for iflag_con = 3
11 c------------------------------------------------------------
12 
13 C
14 C *** PBCRIT IS THE CRITICAL CLOUD DEPTH (MB) BENEATH WHICH THE ***
15 C *** PRECIPITATION EFFICIENCY IS ASSUMED TO BE ZERO ***
16 C *** PTCRIT IS THE CLOUD DEPTH (MB) ABOVE WHICH THE PRECIP. ***
17 C *** EFFICIENCY IS ASSUMED TO BE UNITY ***
18 C *** SIGD IS THE FRACTIONAL AREA COVERED BY UNSATURATED DNDRAFT ***
19 C *** SPFAC IS THE FRACTION OF PRECIPITATION FALLING OUTSIDE ***
20 C *** OF CLOUD ***
21 C
22 C [TAU: CHARACTERISTIC TIMESCALE USED TO COMPUTE ALPHA & BETA]
23 C *** ALPHA AND BETA ARE PARAMETERS THAT CONTROL THE RATE OF ***
24 C *** APPROACH TO QUASI-EQUILIBRIUM ***
25 C *** (THEIR STANDARD VALUES ARE 1.0 AND 0.96, RESPECTIVELY) ***
26 C *** (BETA MUST BE LESS THAN OR EQUAL TO 1) ***
27 C
28 C *** DTCRIT IS THE CRITICAL BUOYANCY (K) USED TO ADJUST THE ***
29 C *** APPROACH TO QUASI-EQUILIBRIUM ***
30 C *** IT MUST BE LESS THAN 0 ***
31 
32 #include "cv30param.h"
33 #include "conema3.h"
34 
35  integer nd
36  real delt ! timestep (seconds)
37 
38 c noff: integer limit for convection (nd-noff)
39 c minorig: First level of convection
40 
41 c -- limit levels for convection:
42 
43  noff = 1
44  minorig = 1
45  nl=nd-noff
46  nlp=nl+1
47  nlm=nl-1
48 
49 c -- "microphysical" parameters:
50 
51  sigd = 0.01
52  spfac = 0.15
53  pbcrit = 150.0
54  ptcrit = 500.0
55 cIM cf. FH epmax = 0.993
56 
57  omtrain = 45.0 ! used also for snow (no disctinction rain/snow)
58 
59 c -- misc:
60 
61  dtovsh = -0.2 ! dT for overshoot
62  dpbase = -40. ! definition cloud base (400m above LCL)
63  dttrig = 5. ! (loose) condition for triggering
64 
65 c -- rate of approach to quasi-equilibrium:
66 
67  dtcrit = -2.0
68  tau = 8000.
69  beta = 1.0 - delt/tau
70  alpha = 1.5e-3 * delt/tau
71 c increase alpha to compensate W decrease:
72  alpha = alpha*1.5
73 
74 c -- interface cloud parameterization:
75 
76  delta=0.01 ! cld
77 
78 c -- interface with boundary-layer (gust factor): (sb)
79 
80  betad=10.0 ! original value (from convect 4.3)
81 
82  return
83  end
84 
85  SUBROUTINE cv30_prelim(len,nd,ndp1,t,q,p,ph
86  : ,lv,cpn,tv,gz,h,hm,th)
87  implicit none
88 
89 !=====================================================================
90 ! --- CALCULATE ARRAYS OF GEOPOTENTIAL, HEAT CAPACITY & STATIC ENERGY
91 ! "ori": from convect4.3 (vectorized)
92 ! "convect3": to be exactly consistent with convect3
93 !=====================================================================
94 
95 c inputs:
96  integer len, nd, ndp1
97  real t(len,nd), q(len,nd), p(len,nd), ph(len,ndp1)
98 
99 c outputs:
100  real lv(len,nd), cpn(len,nd), tv(len,nd)
101  real gz(len,nd), h(len,nd), hm(len,nd)
102  real th(len,nd)
103 
104 c local variables:
105  integer k, i
106  real rdcp
107  real tvx,tvy ! convect3
108  real cpx(len,nd)
109 
110 #include "cvthermo.h"
111 #include "cv30param.h"
112 
113 
114 c ori do 110 k=1,nlp
115  do 110 k=1,nl ! convect3
116  do 100 i=1,len
117 cdebug lv(i,k)= lv0-clmcpv*(t(i,k)-t0)
118  lv(i,k)= lv0-clmcpv*(t(i,k)-273.15)
119  cpn(i,k)=cpd*(1.0-q(i,k))+cpv*q(i,k)
120  cpx(i,k)=cpd*(1.0-q(i,k))+cl*q(i,k)
121 c ori tv(i,k)=t(i,k)*(1.0+q(i,k)*epsim1)
122  tv(i,k)=t(i,k)*(1.0+q(i,k)/eps-q(i,k))
123  rdcp=(rrd*(1.-q(i,k))+q(i,k)*rrv)/cpn(i,k)
124  th(i,k)=t(i,k)*(1000.0/p(i,k))**rdcp
125  100 continue
126  110 continue
127 c
128 c gz = phi at the full levels (same as p).
129 c
130  do 120 i=1,len
131  gz(i,1)=0.0
132  120 continue
133 c ori do 140 k=2,nlp
134  do 140 k=2,nl ! convect3
135  do 130 i=1,len
136  tvx=t(i,k)*(1.+q(i,k)/eps-q(i,k)) !convect3
137  tvy=t(i,k-1)*(1.+q(i,k-1)/eps-q(i,k-1)) !convect3
138  gz(i,k)=gz(i,k-1)+0.5*rrd*(tvx+tvy) !convect3
139  & *(p(i,k-1)-p(i,k))/ph(i,k) !convect3
140 
141 c ori gz(i,k)=gz(i,k-1)+hrd*(tv(i,k-1)+tv(i,k))
142 c ori & *(p(i,k-1)-p(i,k))/ph(i,k)
143  130 continue
144  140 continue
145 c
146 c h = phi + cpT (dry static energy).
147 c hm = phi + cp(T-Tbase)+Lq
148 c
149 c ori do 170 k=1,nlp
150  do 170 k=1,nl ! convect3
151  do 160 i=1,len
152  h(i,k)=gz(i,k)+cpn(i,k)*t(i,k)
153  hm(i,k)=gz(i,k)+cpx(i,k)*(t(i,k)-t(i,1))+lv(i,k)*q(i,k)
154  160 continue
155  170 continue
156 
157  return
158  end
159 
160  SUBROUTINE cv30_feed(len,nd,t,q,qs,p,ph,hm,gz
161  : ,nk,icb,icbmax,iflag,tnk,qnk,gznk,plcl)
162  implicit none
163 
164 C================================================================
165 C Purpose: CONVECTIVE FEED
166 C
167 C Main differences with cv_feed:
168 C - ph added in input
169 C - here, nk(i)=minorig
170 C - icb defined differently (plcl compared with ph instead of p)
171 C
172 C Main differences with convect3:
173 C - we do not compute dplcldt and dplcldr of CLIFT anymore
174 C - values iflag different (but tests identical)
175 C - A,B explicitely defined (!...)
176 C================================================================
177 
178 #include "cv30param.h"
179 
180 c inputs:
181  integer len, nd
182  real t(len,nd), q(len,nd), qs(len,nd), p(len,nd)
183  real hm(len,nd), gz(len,nd)
184  real ph(len,nd+1)
185 
186 c outputs:
187  integer iflag(len), nk(len), icb(len), icbmax
188  real tnk(len), qnk(len), gznk(len), plcl(len)
189 
190 c local variables:
191  integer i, k
192  integer ihmin(len)
193  real work(len)
194  real pnk(len), qsnk(len), rh(len), chi(len)
195  real a, b ! convect3
196 cym
197  plcl=0.0
198 c@ !-------------------------------------------------------------------
199 c@ ! --- Find level of minimum moist static energy
200 c@ ! --- If level of minimum moist static energy coincides with
201 c@ ! --- or is lower than minimum allowable parcel origin level,
202 c@ ! --- set iflag to 6.
203 c@ !-------------------------------------------------------------------
204 c@
205 c@ do 180 i=1,len
206 c@ work(i)=1.0e12
207 c@ ihmin(i)=nl
208 c@ 180 continue
209 c@ do 200 k=2,nlp
210 c@ do 190 i=1,len
211 c@ if((hm(i,k).lt.work(i)).and.
212 c@ & (hm(i,k).lt.hm(i,k-1)))then
213 c@ work(i)=hm(i,k)
214 c@ ihmin(i)=k
215 c@ endif
216 c@ 190 continue
217 c@ 200 continue
218 c@ do 210 i=1,len
219 c@ ihmin(i)=min(ihmin(i),nlm)
220 c@ if(ihmin(i).le.minorig)then
221 c@ iflag(i)=6
222 c@ endif
223 c@ 210 continue
224 c@ c
225 c@ !-------------------------------------------------------------------
226 c@ ! --- Find that model level below the level of minimum moist static
227 c@ ! --- energy that has the maximum value of moist static energy
228 c@ !-------------------------------------------------------------------
229 c@
230 c@ do 220 i=1,len
231 c@ work(i)=hm(i,minorig)
232 c@ nk(i)=minorig
233 c@ 220 continue
234 c@ do 240 k=minorig+1,nl
235 c@ do 230 i=1,len
236 c@ if((hm(i,k).gt.work(i)).and.(k.le.ihmin(i)))then
237 c@ work(i)=hm(i,k)
238 c@ nk(i)=k
239 c@ endif
240 c@ 230 continue
241 c@ 240 continue
242 
243 !-------------------------------------------------------------------
244 ! --- Origin level of ascending parcels for convect3:
245 !-------------------------------------------------------------------
246 
247  do 220 i=1,len
248  nk(i)=minorig
249  220 continue
250 
251 !-------------------------------------------------------------------
252 ! --- Check whether parcel level temperature and specific humidity
253 ! --- are reasonable
254 !-------------------------------------------------------------------
255  do 250 i=1,len
256  if( ( ( t(i,nk(i)).lt.250.0 )
257  & .or.( q(i,nk(i)).le.0.0 ) )
258 c@ & .or.( p(i,ihmin(i)).lt.400.0 ) )
259  & .and.
260  & ( iflag(i).eq.0) ) iflag(i)=7
261  250 continue
262 !-------------------------------------------------------------------
263 ! --- Calculate lifted condensation level of air at parcel origin level
264 ! --- (Within 0.2% of formula of Bolton, MON. WEA. REV.,1980)
265 !-------------------------------------------------------------------
266 
267  a = 1669.0 ! convect3
268  b = 122.0 ! convect3
269 
270  do 260 i=1,len
271 
272  if (iflag(i).ne.7) then ! modif sb Jun7th 2002
273 
274  tnk(i)=t(i,nk(i))
275  qnk(i)=q(i,nk(i))
276  gznk(i)=gz(i,nk(i))
277  pnk(i)=p(i,nk(i))
278  qsnk(i)=qs(i,nk(i))
279 c
280  rh(i)=qnk(i)/qsnk(i)
281 c ori rh(i)=min(1.0,rh(i)) ! removed for convect3
282 c ori chi(i)=tnk(i)/(1669.0-122.0*rh(i)-tnk(i))
283  chi(i)=tnk(i)/(a-b*rh(i)-tnk(i)) ! convect3
284  plcl(i)=pnk(i)*(rh(i)**chi(i))
285  if(((plcl(i).lt.200.0).or.(plcl(i).ge.2000.0))
286  & .and.(iflag(i).eq.0))iflag(i)=8
287 
288  endif ! iflag=7
289 
290  260 continue
291 
292 !-------------------------------------------------------------------
293 ! --- Calculate first level above lcl (=icb)
294 !-------------------------------------------------------------------
295 
296 c@ do 270 i=1,len
297 c@ icb(i)=nlm
298 c@ 270 continue
299 c@c
300 c@ do 290 k=minorig,nl
301 c@ do 280 i=1,len
302 c@ if((k.ge.(nk(i)+1)).and.(p(i,k).lt.plcl(i)))
303 c@ & icb(i)=min(icb(i),k)
304 c@ 280 continue
305 c@ 290 continue
306 c@c
307 c@ do 300 i=1,len
308 c@ if((icb(i).ge.nlm).and.(iflag(i).eq.0))iflag(i)=9
309 c@ 300 continue
310 
311  do 270 i=1,len
312  icb(i)=nlm
313  270 continue
314 c
315 c la modification consiste a comparer plcl a ph et non a p:
316 c icb est defini par : ph(icb)<plcl<ph(icb-1)
317 c@ do 290 k=minorig,nl
318  do 290 k=3,nl-1 ! modif pour que icb soit sup/egal a 2
319  do 280 i=1,len
320  if( ph(i,k).lt.plcl(i) ) icb(i)=min(icb(i),k)
321  280 continue
322  290 continue
323 c
324  do 300 i=1,len
325 c@ if((icb(i).ge.nlm).and.(iflag(i).eq.0))iflag(i)=9
326  if((icb(i).eq.nlm).and.(iflag(i).eq.0))iflag(i)=9
327  300 continue
328 
329  do 400 i=1,len
330  icb(i) = icb(i)-1 ! icb sup ou egal a 2
331  400 continue
332 c
333 c Compute icbmax.
334 c
335  icbmax=2
336  do 310 i=1,len
337 c! icbmax=max(icbmax,icb(i))
338  if (iflag(i).lt.7) icbmax=max(icbmax,icb(i)) ! sb Jun7th02
339  310 continue
340 
341  return
342  end
343 
344  SUBROUTINE cv30_undilute1(len,nd,t,q,qs,gz,plcl,p,nk,icb
345  : ,tp,tvp,clw,icbs)
346  implicit none
347 
348 !----------------------------------------------------------------
349 ! Equivalent de TLIFT entre NK et ICB+1 inclus
350 !
351 ! Differences with convect4:
352 ! - specify plcl in input
353 ! - icbs is the first level above LCL (may differ from icb)
354 ! - in the iterations, used x(icbs) instead x(icb)
355 ! - many minor differences in the iterations
356 ! - tvp is computed in only one time
357 ! - icbs: first level above Plcl (IMIN de TLIFT) in output
358 ! - if icbs=icb, compute also tp(icb+1),tvp(icb+1) & clw(icb+1)
359 !----------------------------------------------------------------
360 
361 #include "cvthermo.h"
362 #include "cv30param.h"
363 
364 c inputs:
365  integer len, nd
366  integer nk(len), icb(len)
367  real t(len,nd), q(len,nd), qs(len,nd), gz(len,nd)
368  real p(len,nd)
369  real plcl(len) ! convect3
370 
371 c outputs:
372  real tp(len,nd), tvp(len,nd), clw(len,nd)
373 
374 c local variables:
375  integer i, k
376  integer icb1(len), icbs(len), icbsmax2 ! convect3
377  real tg, qg, alv, s, ahg, tc, denom, es, rg
378  real ah0(len), cpp(len)
379  real tnk(len), qnk(len), gznk(len), ticb(len), gzicb(len)
380  real qsicb(len) ! convect3
381  real cpinv(len) ! convect3
382 
383 !-------------------------------------------------------------------
384 ! --- Calculates the lifted parcel virtual temperature at nk,
385 ! --- the actual temperature, and the adiabatic
386 ! --- liquid water content. The procedure is to solve the equation.
387 ! cp*tp+L*qp+phi=cp*tnk+L*qnk+gznk.
388 !-------------------------------------------------------------------
389 
390  do 320 i=1,len
391  tnk(i)=t(i,nk(i))
392  qnk(i)=q(i,nk(i))
393  gznk(i)=gz(i,nk(i))
394 c ori ticb(i)=t(i,icb(i))
395 c ori gzicb(i)=gz(i,icb(i))
396  320 continue
397 c
398 c *** Calculate certain parcel quantities, including static energy ***
399 c
400  do 330 i=1,len
401  ah0(i)=(cpd*(1.-qnk(i))+cl*qnk(i))*tnk(i)
402  & +qnk(i)*(lv0-clmcpv*(tnk(i)-273.15))+gznk(i)
403  cpp(i)=cpd*(1.-qnk(i))+qnk(i)*cpv
404  cpinv(i)=1./cpp(i)
405  330 continue
406 c
407 c *** Calculate lifted parcel quantities below cloud base ***
408 c
409  do i=1,len !convect3
410  icb1(i)=max(icb(i),2) !convect3
411  icb1(i)=min(icb(i),nl) !convect3
412 c if icb is below LCL, start loop at ICB+1:
413 c (icbs est le premier niveau au-dessus du LCL)
414  icbs(i)=icb1(i) !convect3
415  if (plcl(i).lt.p(i,icb1(i))) then
416  icbs(i)=min(icbs(i)+1,nl) !convect3
417  endif
418  enddo !convect3
419 
420  do i=1,len !convect3
421  ticb(i)=t(i,icbs(i)) !convect3
422  gzicb(i)=gz(i,icbs(i)) !convect3
423  qsicb(i)=qs(i,icbs(i)) !convect3
424  enddo !convect3
425 
426 c
427 c Re-compute icbsmax (icbsmax2): !convect3
428 c !convect3
429  icbsmax2=2 !convect3
430  do 310 i=1,len !convect3
431  icbsmax2=max(icbsmax2,icbs(i)) !convect3
432  310 continue !convect3
433 
434 c initialization outputs:
435 
436  do k=1,icbsmax2 ! convect3
437  do i=1,len ! convect3
438  tp(i,k) = 0.0 ! convect3
439  tvp(i,k) = 0.0 ! convect3
440  clw(i,k) = 0.0 ! convect3
441  enddo ! convect3
442  enddo ! convect3
443 
444 c tp and tvp below cloud base:
445 
446  do 350 k=minorig,icbsmax2-1
447  do 340 i=1,len
448  tp(i,k)=tnk(i)-(gz(i,k)-gznk(i))*cpinv(i)
449  tvp(i,k)=tp(i,k)*(1.+qnk(i)/eps-qnk(i)) !whole thing (convect3)
450  340 continue
451  350 continue
452 c
453 c *** Find lifted parcel quantities above cloud base ***
454 c
455  do 360 i=1,len
456  tg=ticb(i)
457 c ori qg=qs(i,icb(i))
458  qg=qsicb(i) ! convect3
459 cdebug alv=lv0-clmcpv*(ticb(i)-t0)
460  alv=lv0-clmcpv*(ticb(i)-273.15)
461 c
462 c First iteration.
463 c
464 c ori s=cpd+alv*alv*qg/(rrv*ticb(i)*ticb(i))
465  s=cpd*(1.-qnk(i))+cl*qnk(i) ! convect3
466  : +alv*alv*qg/(rrv*ticb(i)*ticb(i)) ! convect3
467  s=1./s
468 c ori ahg=cpd*tg+(cl-cpd)*qnk(i)*ticb(i)+alv*qg+gzicb(i)
469  ahg=cpd*tg+(cl-cpd)*qnk(i)*tg+alv*qg+gzicb(i) ! convect3
470  tg=tg+s*(ah0(i)-ahg)
471 c ori tg=max(tg,35.0)
472 cdebug tc=tg-t0
473  tc=tg-273.15
474  denom=243.5+tc
475  denom=max(denom,1.0) ! convect3
476 c ori if(tc.ge.0.0)then
477  es=6.112*exp(17.67*tc/denom)
478 c ori else
479 c ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
480 c ori endif
481 c ori qg=eps*es/(p(i,icb(i))-es*(1.-eps))
482  qg=eps*es/(p(i,icbs(i))-es*(1.-eps))
483 c
484 c Second iteration.
485 c
486 
487 c ori s=cpd+alv*alv*qg/(rrv*ticb(i)*ticb(i))
488 c ori s=1./s
489 c ori ahg=cpd*tg+(cl-cpd)*qnk(i)*ticb(i)+alv*qg+gzicb(i)
490  ahg=cpd*tg+(cl-cpd)*qnk(i)*tg+alv*qg+gzicb(i) ! convect3
491  tg=tg+s*(ah0(i)-ahg)
492 c ori tg=max(tg,35.0)
493 cdebug tc=tg-t0
494  tc=tg-273.15
495  denom=243.5+tc
496  denom=max(denom,1.0) ! convect3
497 c ori if(tc.ge.0.0)then
498  es=6.112*exp(17.67*tc/denom)
499 c ori else
500 c ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
501 c ori end if
502 c ori qg=eps*es/(p(i,icb(i))-es*(1.-eps))
503  qg=eps*es/(p(i,icbs(i))-es*(1.-eps))
504 
505  alv=lv0-clmcpv*(ticb(i)-273.15)
506 
507 c ori c approximation here:
508 c ori tp(i,icb(i))=(ah0(i)-(cl-cpd)*qnk(i)*ticb(i)
509 c ori & -gz(i,icb(i))-alv*qg)/cpd
510 
511 c convect3: no approximation:
512  tp(i,icbs(i))=(ah0(i)-gz(i,icbs(i))-alv*qg)
513  : /(cpd+(cl-cpd)*qnk(i))
514 
515 c ori clw(i,icb(i))=qnk(i)-qg
516 c ori clw(i,icb(i))=max(0.0,clw(i,icb(i)))
517  clw(i,icbs(i))=qnk(i)-qg
518  clw(i,icbs(i))=max(0.0,clw(i,icbs(i)))
519 
520  rg=qg/(1.-qnk(i))
521 c ori tvp(i,icb(i))=tp(i,icb(i))*(1.+rg*epsi)
522 c convect3: (qg utilise au lieu du vrai mixing ratio rg)
523  tvp(i,icbs(i))=tp(i,icbs(i))*(1.+qg/eps-qnk(i)) !whole thing
524 
525  360 continue
526 c
527 c ori do 380 k=minorig,icbsmax2
528 c ori do 370 i=1,len
529 c ori tvp(i,k)=tvp(i,k)-tp(i,k)*qnk(i)
530 c ori 370 continue
531 c ori 380 continue
532 c
533 
534 c -- The following is only for convect3:
535 c
536 c * icbs is the first level above the LCL:
537 c if plcl<p(icb), then icbs=icb+1
538 c if plcl>p(icb), then icbs=icb
539 c
540 c * the routine above computes tvp from minorig to icbs (included).
541 c
542 c * to compute buoybase (in cv3_trigger.F), both tvp(icb) and tvp(icb+1)
543 c must be known. This is the case if icbs=icb+1, but not if icbs=icb.
544 c
545 c * therefore, in the case icbs=icb, we compute tvp at level icb+1
546 c (tvp at other levels will be computed in cv3_undilute2.F)
547 c
548 
549  do i=1,len
550  ticb(i)=t(i,icb(i)+1)
551  gzicb(i)=gz(i,icb(i)+1)
552  qsicb(i)=qs(i,icb(i)+1)
553  enddo
554 
555  do 460 i=1,len
556  tg=ticb(i)
557  qg=qsicb(i) ! convect3
558 cdebug alv=lv0-clmcpv*(ticb(i)-t0)
559  alv=lv0-clmcpv*(ticb(i)-273.15)
560 c
561 c First iteration.
562 c
563 c ori s=cpd+alv*alv*qg/(rrv*ticb(i)*ticb(i))
564  s=cpd*(1.-qnk(i))+cl*qnk(i) ! convect3
565  : +alv*alv*qg/(rrv*ticb(i)*ticb(i)) ! convect3
566  s=1./s
567 c ori ahg=cpd*tg+(cl-cpd)*qnk(i)*ticb(i)+alv*qg+gzicb(i)
568  ahg=cpd*tg+(cl-cpd)*qnk(i)*tg+alv*qg+gzicb(i) ! convect3
569  tg=tg+s*(ah0(i)-ahg)
570 c ori tg=max(tg,35.0)
571 cdebug tc=tg-t0
572  tc=tg-273.15
573  denom=243.5+tc
574  denom=max(denom,1.0) ! convect3
575 c ori if(tc.ge.0.0)then
576  es=6.112*exp(17.67*tc/denom)
577 c ori else
578 c ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
579 c ori endif
580 c ori qg=eps*es/(p(i,icb(i))-es*(1.-eps))
581  qg=eps*es/(p(i,icb(i)+1)-es*(1.-eps))
582 c
583 c Second iteration.
584 c
585 
586 c ori s=cpd+alv*alv*qg/(rrv*ticb(i)*ticb(i))
587 c ori s=1./s
588 c ori ahg=cpd*tg+(cl-cpd)*qnk(i)*ticb(i)+alv*qg+gzicb(i)
589  ahg=cpd*tg+(cl-cpd)*qnk(i)*tg+alv*qg+gzicb(i) ! convect3
590  tg=tg+s*(ah0(i)-ahg)
591 c ori tg=max(tg,35.0)
592 cdebug tc=tg-t0
593  tc=tg-273.15
594  denom=243.5+tc
595  denom=max(denom,1.0) ! convect3
596 c ori if(tc.ge.0.0)then
597  es=6.112*exp(17.67*tc/denom)
598 c ori else
599 c ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
600 c ori end if
601 c ori qg=eps*es/(p(i,icb(i))-es*(1.-eps))
602  qg=eps*es/(p(i,icb(i)+1)-es*(1.-eps))
603 
604  alv=lv0-clmcpv*(ticb(i)-273.15)
605 
606 c ori c approximation here:
607 c ori tp(i,icb(i))=(ah0(i)-(cl-cpd)*qnk(i)*ticb(i)
608 c ori & -gz(i,icb(i))-alv*qg)/cpd
609 
610 c convect3: no approximation:
611  tp(i,icb(i)+1)=(ah0(i)-gz(i,icb(i)+1)-alv*qg)
612  : /(cpd+(cl-cpd)*qnk(i))
613 
614 c ori clw(i,icb(i))=qnk(i)-qg
615 c ori clw(i,icb(i))=max(0.0,clw(i,icb(i)))
616  clw(i,icb(i)+1)=qnk(i)-qg
617  clw(i,icb(i)+1)=max(0.0,clw(i,icb(i)+1))
618 
619  rg=qg/(1.-qnk(i))
620 c ori tvp(i,icb(i))=tp(i,icb(i))*(1.+rg*epsi)
621 c convect3: (qg utilise au lieu du vrai mixing ratio rg)
622  tvp(i,icb(i)+1)=tp(i,icb(i)+1)*(1.+qg/eps-qnk(i)) !whole thing
623 
624  460 continue
625 
626  return
627  end
628 
629  SUBROUTINE cv30_trigger(len,nd,icb,plcl,p,th,tv,tvp
630  o ,pbase,buoybase,iflag,sig,w0)
631  implicit none
632 
633 !-------------------------------------------------------------------
634 ! --- TRIGGERING
635 !
636 ! - computes the cloud base
637 ! - triggering (crude in this version)
638 ! - relaxation of sig and w0 when no convection
639 !
640 ! Caution1: if no convection, we set iflag=4
641 ! (it used to be 0 in convect3)
642 !
643 ! Caution2: at this stage, tvp (and thus buoy) are know up
644 ! through icb only!
645 ! -> the buoyancy below cloud base not (yet) set to the cloud base buoyancy
646 !-------------------------------------------------------------------
647 
648 #include "cv30param.h"
649 
650 c input:
651  integer len, nd
652  integer icb(len)
653  real plcl(len), p(len,nd)
654  real th(len,nd), tv(len,nd), tvp(len,nd)
655 
656 c output:
657  real pbase(len), buoybase(len)
658 
659 c input AND output:
660  integer iflag(len)
661  real sig(len,nd), w0(len,nd)
662 
663 c local variables:
664  integer i,k
665  real tvpbase, tvbase, tdif, ath, ath1
666 
667 c
668 c *** set cloud base buoyancy at (plcl+dpbase) level buoyancy
669 c
670  do 100 i=1,len
671  pbase(i) = plcl(i) + dpbase
672  tvpbase = tvp(i,icb(i))*(pbase(i)-p(i,icb(i)+1))
673  : /(p(i,icb(i))-p(i,icb(i)+1))
674  : + tvp(i,icb(i)+1)*(p(i,icb(i))-pbase(i))
675  : /(p(i,icb(i))-p(i,icb(i)+1))
676  tvbase = tv(i,icb(i))*(pbase(i)-p(i,icb(i)+1))
677  : /(p(i,icb(i))-p(i,icb(i)+1))
678  : + tv(i,icb(i)+1)*(p(i,icb(i))-pbase(i))
679  : /(p(i,icb(i))-p(i,icb(i)+1))
680  buoybase(i) = tvpbase - tvbase
681 100 continue
682 
683 c
684 c *** make sure that column is dry adiabatic between the surface ***
685 c *** and cloud base, and that lifted air is positively buoyant ***
686 c *** at cloud base ***
687 c *** if not, return to calling program after resetting ***
688 c *** sig(i) and w0(i) ***
689 c
690 
691 c oct3 do 200 i=1,len
692 c oct3
693 c oct3 tdif = buoybase(i)
694 c oct3 ath1 = th(i,1)
695 c oct3 ath = th(i,icb(i)-1) - dttrig
696 c oct3
697 c oct3 if (tdif.lt.dtcrit .or. ath.gt.ath1) then
698 c oct3 do 60 k=1,nl
699 c oct3 sig(i,k) = beta*sig(i,k) - 2.*alpha*tdif*tdif
700 c oct3 sig(i,k) = AMAX1(sig(i,k),0.0)
701 c oct3 w0(i,k) = beta*w0(i,k)
702 c oct3 60 continue
703 c oct3 iflag(i)=4 ! pour version vectorisee
704 c oct3c convect3 iflag(i)=0
705 c oct3cccc return
706 c oct3 endif
707 c oct3
708 c oct3200 continue
709 
710 c -- oct3: on reecrit la boucle 200 (pour la vectorisation)
711 
712  do 60 k=1,nl
713  do 200 i=1,len
714 
715  tdif = buoybase(i)
716  ath1 = th(i,1)
717  ath = th(i,icb(i)-1) - dttrig
718 
719  if (tdif.lt.dtcrit .or. ath.gt.ath1) then
720  sig(i,k) = beta*sig(i,k) - 2.*alpha*tdif*tdif
721  sig(i,k) = amax1(sig(i,k),0.0)
722  w0(i,k) = beta*w0(i,k)
723  iflag(i)=4 ! pour version vectorisee
724 c convect3 iflag(i)=0
725  endif
726 
727 200 continue
728  60 continue
729 
730 c fin oct3 --
731 
732  return
733  end
734 
735  SUBROUTINE cv30_compress( len,nloc,ncum,nd,ntra
736  : ,iflag1,nk1,icb1,icbs1
737  : ,plcl1,tnk1,qnk1,gznk1,pbase1,buoybase1
738  : ,t1,q1,qs1,u1,v1,gz1,th1
739  : ,tra1
740  : ,h1,lv1,cpn1,p1,ph1,tv1,tp1,tvp1,clw1
741  : ,sig1,w01
742  o ,iflag,nk,icb,icbs
743  o ,plcl,tnk,qnk,gznk,pbase,buoybase
744  o ,t,q,qs,u,v,gz,th
745  o ,tra
746  o ,h,lv,cpn,p,ph,tv,tp,tvp,clw
747  o ,sig,w0 )
748  implicit none
749 
750 #include "cv30param.h"
751  include 'iniprint.h'
752 
753 c inputs:
754  integer len,ncum,nd,ntra,nloc
755  integer iflag1(len),nk1(len),icb1(len),icbs1(len)
756  real plcl1(len),tnk1(len),qnk1(len),gznk1(len)
757  real pbase1(len),buoybase1(len)
758  real t1(len,nd),q1(len,nd),qs1(len,nd),u1(len,nd),v1(len,nd)
759  real gz1(len,nd),h1(len,nd),lv1(len,nd),cpn1(len,nd)
760  real p1(len,nd),ph1(len,nd+1),tv1(len,nd),tp1(len,nd)
761  real tvp1(len,nd),clw1(len,nd)
762  real th1(len,nd)
763  real sig1(len,nd), w01(len,nd)
764  real tra1(len,nd,ntra)
765 
766 c outputs:
767 c en fait, on a nloc=len pour l'instant (cf cv_driver)
768  integer iflag(nloc),nk(nloc),icb(nloc),icbs(nloc)
769  real plcl(nloc),tnk(nloc),qnk(nloc),gznk(nloc)
770  real pbase(nloc),buoybase(nloc)
771  real t(nloc,nd),q(nloc,nd),qs(nloc,nd),u(nloc,nd),v(nloc,nd)
772  real gz(nloc,nd),h(nloc,nd),lv(nloc,nd),cpn(nloc,nd)
773  real p(nloc,nd),ph(nloc,nd+1),tv(nloc,nd),tp(nloc,nd)
774  real tvp(nloc,nd),clw(nloc,nd)
775  real th(nloc,nd)
776  real sig(nloc,nd), w0(nloc,nd)
777  real tra(nloc,nd,ntra)
778 
779 c local variables:
780  integer i,k,nn,j
781 
782  CHARACTER (LEN=20) :: modname='cv30_compress'
783  CHARACTER (LEN=80) :: abort_message
784 
785 
786  do 110 k=1,nl+1
787  nn=0
788  do 100 i=1,len
789  if(iflag1(i).eq.0)then
790  nn=nn+1
791  sig(nn,k)=sig1(i,k)
792  w0(nn,k)=w01(i,k)
793  t(nn,k)=t1(i,k)
794  q(nn,k)=q1(i,k)
795  qs(nn,k)=qs1(i,k)
796  u(nn,k)=u1(i,k)
797  v(nn,k)=v1(i,k)
798  gz(nn,k)=gz1(i,k)
799  h(nn,k)=h1(i,k)
800  lv(nn,k)=lv1(i,k)
801  cpn(nn,k)=cpn1(i,k)
802  p(nn,k)=p1(i,k)
803  ph(nn,k)=ph1(i,k)
804  tv(nn,k)=tv1(i,k)
805  tp(nn,k)=tp1(i,k)
806  tvp(nn,k)=tvp1(i,k)
807  clw(nn,k)=clw1(i,k)
808  th(nn,k)=th1(i,k)
809  endif
810  100 continue
811  110 continue
812 
813 c do 121 j=1,ntra
814 c do 111 k=1,nd
815 c nn=0
816 c do 101 i=1,len
817 c if(iflag1(i).eq.0)then
818 c nn=nn+1
819 c tra(nn,k,j)=tra1(i,k,j)
820 c endif
821 c 101 continue
822 c 111 continue
823 c 121 continue
824 
825  if (nn.ne.ncum) then
826  write(lunout,*)'strange! nn not equal to ncum: ',nn,ncum
827  abort_message = ''
828  CALL abort_gcm(modname,abort_message,1)
829  endif
830 
831  nn=0
832  do 150 i=1,len
833  if(iflag1(i).eq.0)then
834  nn=nn+1
835  pbase(nn)=pbase1(i)
836  buoybase(nn)=buoybase1(i)
837  plcl(nn)=plcl1(i)
838  tnk(nn)=tnk1(i)
839  qnk(nn)=qnk1(i)
840  gznk(nn)=gznk1(i)
841  nk(nn)=nk1(i)
842  icb(nn)=icb1(i)
843  icbs(nn)=icbs1(i)
844  iflag(nn)=iflag1(i)
845  endif
846  150 continue
847 
848  return
849  end
850 
851  SUBROUTINE cv30_undilute2(nloc,ncum,nd,icb,icbs,nk
852  : ,tnk,qnk,gznk,t,q,qs,gz
853  : ,p,h,tv,lv,pbase,buoybase,plcl
854  o ,inb,tp,tvp,clw,hp,ep,sigp,buoy)
855  implicit none
856 
857 C---------------------------------------------------------------------
858 C Purpose:
859 C FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
860 C &
861 C COMPUTE THE PRECIPITATION EFFICIENCIES AND THE
862 C FRACTION OF PRECIPITATION FALLING OUTSIDE OF CLOUD
863 C &
864 C FIND THE LEVEL OF NEUTRAL BUOYANCY
865 C
866 C Main differences convect3/convect4:
867 C - icbs (input) is the first level above LCL (may differ from icb)
868 C - many minor differences in the iterations
869 C - condensed water not removed from tvp in convect3
870 C - vertical profile of buoyancy computed here (use of buoybase)
871 C - the determination of inb is different
872 C - no inb1, only inb in output
873 C---------------------------------------------------------------------
874 
875 #include "cvthermo.h"
876 #include "cv30param.h"
877 #include "conema3.h"
878 
879 c inputs:
880  integer ncum, nd, nloc
881  integer icb(nloc), icbs(nloc), nk(nloc)
882  real t(nloc,nd), q(nloc,nd), qs(nloc,nd), gz(nloc,nd)
883  real p(nloc,nd)
884  real tnk(nloc), qnk(nloc), gznk(nloc)
885  real lv(nloc,nd), tv(nloc,nd), h(nloc,nd)
886  real pbase(nloc), buoybase(nloc), plcl(nloc)
887 
888 c outputs:
889  integer inb(nloc)
890  real tp(nloc,nd), tvp(nloc,nd), clw(nloc,nd)
891  real ep(nloc,nd), sigp(nloc,nd), hp(nloc,nd)
892  real buoy(nloc,nd)
893 
894 c local variables:
895  integer i, k
896  real tg,qg,ahg,alv,s,tc,es,denom,rg,tca,elacrit
897  real by, defrac, pden
898  real ah0(nloc), cape(nloc), capem(nloc), byp(nloc)
899  logical lcape(nloc)
900 
901 !=====================================================================
902 ! --- SOME INITIALIZATIONS
903 !=====================================================================
904 
905  do 170 k=1,nl
906  do 160 i=1,ncum
907  ep(i,k)=0.0
908  sigp(i,k)=spfac
909  160 continue
910  170 continue
911 
912 !=====================================================================
913 ! --- FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
914 !=====================================================================
915 c
916 c --- The procedure is to solve the equation.
917 c cp*tp+L*qp+phi=cp*tnk+L*qnk+gznk.
918 c
919 c *** Calculate certain parcel quantities, including static energy ***
920 c
921 c
922  do 240 i=1,ncum
923  ah0(i)=(cpd*(1.-qnk(i))+cl*qnk(i))*tnk(i)
924 cdebug & +qnk(i)*(lv0-clmcpv*(tnk(i)-t0))+gznk(i)
925  & +qnk(i)*(lv0-clmcpv*(tnk(i)-273.15))+gznk(i)
926  240 continue
927 c
928 c
929 c *** Find lifted parcel quantities above cloud base ***
930 c
931 c
932  do 300 k=minorig+1,nl
933  do 290 i=1,ncum
934 c ori if(k.ge.(icb(i)+1))then
935  if(k.ge.(icbs(i)+1))then ! convect3
936  tg=t(i,k)
937  qg=qs(i,k)
938 cdebug alv=lv0-clmcpv*(t(i,k)-t0)
939  alv=lv0-clmcpv*(t(i,k)-273.15)
940 c
941 c First iteration.
942 c
943 c ori s=cpd+alv*alv*qg/(rrv*t(i,k)*t(i,k))
944  s=cpd*(1.-qnk(i))+cl*qnk(i) ! convect3
945  : +alv*alv*qg/(rrv*t(i,k)*t(i,k)) ! convect3
946  s=1./s
947 c ori ahg=cpd*tg+(cl-cpd)*qnk(i)*t(i,k)+alv*qg+gz(i,k)
948  ahg=cpd*tg+(cl-cpd)*qnk(i)*tg+alv*qg+gz(i,k) ! convect3
949  tg=tg+s*(ah0(i)-ahg)
950 c ori tg=max(tg,35.0)
951 cdebug tc=tg-t0
952  tc=tg-273.15
953  denom=243.5+tc
954  denom=max(denom,1.0) ! convect3
955 c ori if(tc.ge.0.0)then
956  es=6.112*exp(17.67*tc/denom)
957 c ori else
958 c ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
959 c ori endif
960  qg=eps*es/(p(i,k)-es*(1.-eps))
961 c
962 c Second iteration.
963 c
964 c ori s=cpd+alv*alv*qg/(rrv*t(i,k)*t(i,k))
965 c ori s=1./s
966 c ori ahg=cpd*tg+(cl-cpd)*qnk(i)*t(i,k)+alv*qg+gz(i,k)
967  ahg=cpd*tg+(cl-cpd)*qnk(i)*tg+alv*qg+gz(i,k) ! convect3
968  tg=tg+s*(ah0(i)-ahg)
969 c ori tg=max(tg,35.0)
970 cdebug tc=tg-t0
971  tc=tg-273.15
972  denom=243.5+tc
973  denom=max(denom,1.0) ! convect3
974 c ori if(tc.ge.0.0)then
975  es=6.112*exp(17.67*tc/denom)
976 c ori else
977 c ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
978 c ori endif
979  qg=eps*es/(p(i,k)-es*(1.-eps))
980 c
981 cdebug alv=lv0-clmcpv*(t(i,k)-t0)
982  alv=lv0-clmcpv*(t(i,k)-273.15)
983 c print*,'cpd dans convect2 ',cpd
984 c print*,'tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd'
985 c print*,tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd
986 
987 c ori c approximation here:
988 c ori tp(i,k)=(ah0(i)-(cl-cpd)*qnk(i)*t(i,k)-gz(i,k)-alv*qg)/cpd
989 
990 c convect3: no approximation:
991  tp(i,k)=(ah0(i)-gz(i,k)-alv*qg)/(cpd+(cl-cpd)*qnk(i))
992 
993  clw(i,k)=qnk(i)-qg
994  clw(i,k)=max(0.0,clw(i,k))
995  rg=qg/(1.-qnk(i))
996 c ori tvp(i,k)=tp(i,k)*(1.+rg*epsi)
997 c convect3: (qg utilise au lieu du vrai mixing ratio rg):
998  tvp(i,k)=tp(i,k)*(1.+qg/eps-qnk(i)) ! whole thing
999  endif
1000  290 continue
1001  300 continue
1002 c
1003 !=====================================================================
1004 ! --- SET THE PRECIPITATION EFFICIENCIES AND THE FRACTION OF
1005 ! --- PRECIPITATION FALLING OUTSIDE OF CLOUD
1006 ! --- THESE MAY BE FUNCTIONS OF TP(I), P(I) AND CLW(I)
1007 !=====================================================================
1008 c
1009 c ori do 320 k=minorig+1,nl
1010  do 320 k=1,nl ! convect3
1011  do 310 i=1,ncum
1012  pden=ptcrit-pbcrit
1013  ep(i,k)=(plcl(i)-p(i,k)-pbcrit)/pden*epmax
1014  ep(i,k)=amax1(ep(i,k),0.0)
1015  ep(i,k)=amin1(ep(i,k),epmax)
1016  sigp(i,k)=spfac
1017 c ori if(k.ge.(nk(i)+1))then
1018 c ori tca=tp(i,k)-t0
1019 c ori if(tca.ge.0.0)then
1020 c ori elacrit=elcrit
1021 c ori else
1022 c ori elacrit=elcrit*(1.0-tca/tlcrit)
1023 c ori endif
1024 c ori elacrit=max(elacrit,0.0)
1025 c ori ep(i,k)=1.0-elacrit/max(clw(i,k),1.0e-8)
1026 c ori ep(i,k)=max(ep(i,k),0.0 )
1027 c ori ep(i,k)=min(ep(i,k),1.0 )
1028 c ori sigp(i,k)=sigs
1029 c ori endif
1030  310 continue
1031  320 continue
1032 c
1033 !=====================================================================
1034 ! --- CALCULATE VIRTUAL TEMPERATURE AND LIFTED PARCEL
1035 ! --- VIRTUAL TEMPERATURE
1036 !=====================================================================
1037 c
1038 c dans convect3, tvp est calcule en une seule fois, et sans retirer
1039 c l'eau condensee (~> reversible CAPE)
1040 c
1041 c ori do 340 k=minorig+1,nl
1042 c ori do 330 i=1,ncum
1043 c ori if(k.ge.(icb(i)+1))then
1044 c ori tvp(i,k)=tvp(i,k)*(1.0-qnk(i)+ep(i,k)*clw(i,k))
1045 c oric print*,'i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)'
1046 c oric print*, i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)
1047 c ori endif
1048 c ori 330 continue
1049 c ori 340 continue
1050 
1051 c ori do 350 i=1,ncum
1052 c ori tvp(i,nlp)=tvp(i,nl)-(gz(i,nlp)-gz(i,nl))/cpd
1053 c ori 350 continue
1054 
1055  do 350 i=1,ncum ! convect3
1056  tp(i,nlp)=tp(i,nl) ! convect3
1057  350 continue ! convect3
1058 c
1059 c=====================================================================
1060 c --- EFFECTIVE VERTICAL PROFILE OF BUOYANCY (convect3 only):
1061 c=====================================================================
1062 
1063 c-- this is for convect3 only:
1064 
1065 c first estimate of buoyancy:
1066 
1067  do 500 i=1,ncum
1068  do 501 k=1,nl
1069  buoy(i,k)=tvp(i,k)-tv(i,k)
1070  501 continue
1071  500 continue
1072 
1073 c set buoyancy=buoybase for all levels below base
1074 c for safety, set buoy(icb)=buoybase
1075 
1076  do 505 i=1,ncum
1077  do 506 k=1,nl
1078  if((k.ge.icb(i)).and.(k.le.nl).and.(p(i,k).ge.pbase(i)))then
1079  buoy(i,k)=buoybase(i)
1080  endif
1081  506 continue
1082 cIM cf. CRio/JYG 270807 buoy(icb(i),k)=buoybase(i)
1083  buoy(i,icb(i))=buoybase(i)
1084  505 continue
1085 
1086 c-- end convect3
1087 
1088 c=====================================================================
1089 c --- FIND THE FIRST MODEL LEVEL (INB) ABOVE THE PARCEL'S
1090 c --- LEVEL OF NEUTRAL BUOYANCY
1091 c=====================================================================
1092 c
1093 c-- this is for convect3 only:
1094 
1095  do 510 i=1,ncum
1096  inb(i)=nl-1
1097  510 continue
1098 
1099  do 530 i=1,ncum
1100  do 535 k=1,nl-1
1101  if ((k.ge.icb(i)).and.(buoy(i,k).lt.dtovsh)) then
1102  inb(i)=min(inb(i),k)
1103  endif
1104  535 continue
1105  530 continue
1106 
1107 c-- end convect3
1108 
1109 c ori do 510 i=1,ncum
1110 c ori cape(i)=0.0
1111 c ori capem(i)=0.0
1112 c ori inb(i)=icb(i)+1
1113 c ori inb1(i)=inb(i)
1114 c ori 510 continue
1115 c
1116 c Originial Code
1117 c
1118 c do 530 k=minorig+1,nl-1
1119 c do 520 i=1,ncum
1120 c if(k.ge.(icb(i)+1))then
1121 c by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
1122 c byp=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
1123 c cape(i)=cape(i)+by
1124 c if(by.ge.0.0)inb1(i)=k+1
1125 c if(cape(i).gt.0.0)then
1126 c inb(i)=k+1
1127 c capem(i)=cape(i)
1128 c endif
1129 c endif
1130 c520 continue
1131 c530 continue
1132 c do 540 i=1,ncum
1133 c byp=(tvp(i,nl)-tv(i,nl))*dph(i,nl)/p(i,nl)
1134 c cape(i)=capem(i)+byp
1135 c defrac=capem(i)-cape(i)
1136 c defrac=max(defrac,0.001)
1137 c frac(i)=-cape(i)/defrac
1138 c frac(i)=min(frac(i),1.0)
1139 c frac(i)=max(frac(i),0.0)
1140 c540 continue
1141 c
1142 c K Emanuel fix
1143 c
1144 c call zilch(byp,ncum)
1145 c do 530 k=minorig+1,nl-1
1146 c do 520 i=1,ncum
1147 c if(k.ge.(icb(i)+1))then
1148 c by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
1149 c cape(i)=cape(i)+by
1150 c if(by.ge.0.0)inb1(i)=k+1
1151 c if(cape(i).gt.0.0)then
1152 c inb(i)=k+1
1153 c capem(i)=cape(i)
1154 c byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
1155 c endif
1156 c endif
1157 c520 continue
1158 c530 continue
1159 c do 540 i=1,ncum
1160 c inb(i)=max(inb(i),inb1(i))
1161 c cape(i)=capem(i)+byp(i)
1162 c defrac=capem(i)-cape(i)
1163 c defrac=max(defrac,0.001)
1164 c frac(i)=-cape(i)/defrac
1165 c frac(i)=min(frac(i),1.0)
1166 c frac(i)=max(frac(i),0.0)
1167 c540 continue
1168 c
1169 c J Teixeira fix
1170 c
1171 c ori call zilch(byp,ncum)
1172 c ori do 515 i=1,ncum
1173 c ori lcape(i)=.true.
1174 c ori 515 continue
1175 c ori do 530 k=minorig+1,nl-1
1176 c ori do 520 i=1,ncum
1177 c ori if(cape(i).lt.0.0)lcape(i)=.false.
1178 c ori if((k.ge.(icb(i)+1)).and.lcape(i))then
1179 c ori by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
1180 c ori byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
1181 c ori cape(i)=cape(i)+by
1182 c ori if(by.ge.0.0)inb1(i)=k+1
1183 c ori if(cape(i).gt.0.0)then
1184 c ori inb(i)=k+1
1185 c ori capem(i)=cape(i)
1186 c ori endif
1187 c ori endif
1188 c ori 520 continue
1189 c ori 530 continue
1190 c ori do 540 i=1,ncum
1191 c ori cape(i)=capem(i)+byp(i)
1192 c ori defrac=capem(i)-cape(i)
1193 c ori defrac=max(defrac,0.001)
1194 c ori frac(i)=-cape(i)/defrac
1195 c ori frac(i)=min(frac(i),1.0)
1196 c ori frac(i)=max(frac(i),0.0)
1197 c ori 540 continue
1198 c
1199 c=====================================================================
1200 c --- CALCULATE LIQUID WATER STATIC ENERGY OF LIFTED PARCEL
1201 c=====================================================================
1202 c
1203 cym do i=1,ncum*nlp
1204 cym hp(i,1)=h(i,1)
1205 cym enddo
1206 
1207  do k=1,nlp
1208  do i=1,ncum
1209  hp(i,k)=h(i,k)
1210  enddo
1211  enddo
1212 
1213  do 600 k=minorig+1,nl
1214  do 590 i=1,ncum
1215  if((k.ge.icb(i)).and.(k.le.inb(i)))then
1216  hp(i,k)=h(i,nk(i))+(lv(i,k)+(cpd-cpv)*t(i,k))*ep(i,k)*clw(i,k)
1217  endif
1218  590 continue
1219  600 continue
1220 
1221  return
1222  end
1223 
1224  SUBROUTINE cv30_closure(nloc,ncum,nd,icb,inb
1225  : ,pbase,p,ph,tv,buoy
1226  o ,sig,w0,cape,m)
1227  implicit none
1228 
1229 !===================================================================
1230 ! --- CLOSURE OF CONVECT3
1231 !
1232 ! vectorization: S. Bony
1233 !===================================================================
1234 
1235 #include "cvthermo.h"
1236 #include "cv30param.h"
1237 
1238 c input:
1239  integer ncum, nd, nloc
1240  integer icb(nloc), inb(nloc)
1241  real pbase(nloc)
1242  real p(nloc,nd), ph(nloc,nd+1)
1243  real tv(nloc,nd), buoy(nloc,nd)
1244 
1245 c input/output:
1246  real sig(nloc,nd), w0(nloc,nd)
1247 
1248 c output:
1249  real cape(nloc)
1250  real m(nloc,nd)
1251 
1252 c local variables:
1253  integer i, j, k, icbmax
1254  real deltap, fac, w, amu
1255  real dtmin(nloc,nd), sigold(nloc,nd)
1256 
1257 
1258 c -------------------------------------------------------
1259 c -- Initialization
1260 c -------------------------------------------------------
1261 
1262  do k=1,nl
1263  do i=1,ncum
1264  m(i,k)=0.0
1265  enddo
1266  enddo
1267 
1268 c -------------------------------------------------------
1269 c -- Reset sig(i) and w0(i) for i>inb and i<icb
1270 c -------------------------------------------------------
1271 
1272 c update sig and w0 above LNB:
1273 
1274  do 100 k=1,nl-1
1275  do 110 i=1,ncum
1276  if ((inb(i).lt.(nl-1)).and.(k.ge.(inb(i)+1)))then
1277  sig(i,k)=beta*sig(i,k)
1278  : +2.*alpha*buoy(i,inb(i))*abs(buoy(i,inb(i)))
1279  sig(i,k)=amax1(sig(i,k),0.0)
1280  w0(i,k)=beta*w0(i,k)
1281  endif
1282  110 continue
1283  100 continue
1284 
1285 c compute icbmax:
1286 
1287  icbmax=2
1288  do 200 i=1,ncum
1289  icbmax=max(icbmax,icb(i))
1290  200 continue
1291 
1292 c update sig and w0 below cloud base:
1293 
1294  do 300 k=1,icbmax
1295  do 310 i=1,ncum
1296  if (k.le.icb(i))then
1297  sig(i,k)=beta*sig(i,k)-2.*alpha*buoy(i,icb(i))*buoy(i,icb(i))
1298  sig(i,k)=amax1(sig(i,k),0.0)
1299  w0(i,k)=beta*w0(i,k)
1300  endif
1301 310 continue
1302 300 continue
1303 
1304 c! if(inb.lt.(nl-1))then
1305 c! do 85 i=inb+1,nl-1
1306 c! sig(i)=beta*sig(i)+2.*alpha*buoy(inb)*
1307 c! 1 abs(buoy(inb))
1308 c! sig(i)=amax1(sig(i),0.0)
1309 c! w0(i)=beta*w0(i)
1310 c! 85 continue
1311 c! end if
1312 
1313 c! do 87 i=1,icb
1314 c! sig(i)=beta*sig(i)-2.*alpha*buoy(icb)*buoy(icb)
1315 c! sig(i)=amax1(sig(i),0.0)
1316 c! w0(i)=beta*w0(i)
1317 c! 87 continue
1318 
1319 c -------------------------------------------------------------
1320 c -- Reset fractional areas of updrafts and w0 at initial time
1321 c -- and after 10 time steps of no convection
1322 c -------------------------------------------------------------
1323 
1324  do 400 k=1,nl-1
1325  do 410 i=1,ncum
1326  if (sig(i,nd).lt.1.5.or.sig(i,nd).gt.12.0)then
1327  sig(i,k)=0.0
1328  w0(i,k)=0.0
1329  endif
1330  410 continue
1331  400 continue
1332 
1333 c -------------------------------------------------------------
1334 c -- Calculate convective available potential energy (cape),
1335 c -- vertical velocity (w), fractional area covered by
1336 c -- undilute updraft (sig), and updraft mass flux (m)
1337 c -------------------------------------------------------------
1338 
1339  do 500 i=1,ncum
1340  cape(i)=0.0
1341  500 continue
1342 
1343 c compute dtmin (minimum buoyancy between ICB and given level k):
1344 
1345  do i=1,ncum
1346  do k=1,nl
1347  dtmin(i,k)=100.0
1348  enddo
1349  enddo
1350 
1351  do 550 i=1,ncum
1352  do 560 k=1,nl
1353  do 570 j=minorig,nl
1354  if ( (k.ge.(icb(i)+1)).and.(k.le.inb(i)).and.
1355  : (j.ge.icb(i)).and.(j.le.(k-1)) )then
1356  dtmin(i,k)=amin1(dtmin(i,k),buoy(i,j))
1357  endif
1358  570 continue
1359  560 continue
1360  550 continue
1361 
1362 c the interval on which cape is computed starts at pbase :
1363 
1364  do 600 k=1,nl
1365  do 610 i=1,ncum
1366 
1367  if ((k.ge.(icb(i)+1)).and.(k.le.inb(i))) then
1368 
1369  deltap = min(pbase(i),ph(i,k-1))-min(pbase(i),ph(i,k))
1370  cape(i)=cape(i)+rrd*buoy(i,k-1)*deltap/p(i,k-1)
1371  cape(i)=amax1(0.0,cape(i))
1372  sigold(i,k)=sig(i,k)
1373 
1374 c dtmin(i,k)=100.0
1375 c do 97 j=icb(i),k-1 ! mauvaise vectorisation
1376 c dtmin(i,k)=AMIN1(dtmin(i,k),buoy(i,j))
1377 c 97 continue
1378 
1379  sig(i,k)=beta*sig(i,k)+alpha*dtmin(i,k)*abs(dtmin(i,k))
1380  sig(i,k)=amax1(sig(i,k),0.0)
1381  sig(i,k)=amin1(sig(i,k),0.01)
1382  fac=amin1(((dtcrit-dtmin(i,k))/dtcrit),1.0)
1383  w=(1.-beta)*fac*sqrt(cape(i))+beta*w0(i,k)
1384  amu=0.5*(sig(i,k)+sigold(i,k))*w
1385  m(i,k)=amu*0.007*p(i,k)*(ph(i,k)-ph(i,k+1))/tv(i,k)
1386  w0(i,k)=w
1387  endif
1388 
1389  610 continue
1390  600 continue
1391 
1392  do 700 i=1,ncum
1393  w0(i,icb(i))=0.5*w0(i,icb(i)+1)
1394  m(i,icb(i))=0.5*m(i,icb(i)+1)
1395  : *(ph(i,icb(i))-ph(i,icb(i)+1))
1396  : /(ph(i,icb(i)+1)-ph(i,icb(i)+2))
1397  sig(i,icb(i))=sig(i,icb(i)+1)
1398  sig(i,icb(i)-1)=sig(i,icb(i))
1399  700 continue
1400 
1401 
1402 c! cape=0.0
1403 c! do 98 i=icb+1,inb
1404 c! deltap = min(pbase,ph(i-1))-min(pbase,ph(i))
1405 c! cape=cape+rrd*buoy(i-1)*deltap/p(i-1)
1406 c! dcape=rrd*buoy(i-1)*deltap/p(i-1)
1407 c! dlnp=deltap/p(i-1)
1408 c! cape=amax1(0.0,cape)
1409 c! sigold=sig(i)
1410 
1411 c! dtmin=100.0
1412 c! do 97 j=icb,i-1
1413 c! dtmin=amin1(dtmin,buoy(j))
1414 c! 97 continue
1415 
1416 c! sig(i)=beta*sig(i)+alpha*dtmin*abs(dtmin)
1417 c! sig(i)=amax1(sig(i),0.0)
1418 c! sig(i)=amin1(sig(i),0.01)
1419 c! fac=amin1(((dtcrit-dtmin)/dtcrit),1.0)
1420 c! w=(1.-beta)*fac*sqrt(cape)+beta*w0(i)
1421 c! amu=0.5*(sig(i)+sigold)*w
1422 c! m(i)=amu*0.007*p(i)*(ph(i)-ph(i+1))/tv(i)
1423 c! w0(i)=w
1424 c! 98 continue
1425 c! w0(icb)=0.5*w0(icb+1)
1426 c! m(icb)=0.5*m(icb+1)*(ph(icb)-ph(icb+1))/(ph(icb+1)-ph(icb+2))
1427 c! sig(icb)=sig(icb+1)
1428 c! sig(icb-1)=sig(icb)
1429 
1430  return
1431  end
1432 
1433  SUBROUTINE cv30_mixing(nloc,ncum,nd,na,ntra,icb,nk,inb
1434  : ,ph,t,rr,rs,u,v,tra,h,lv,qnk
1435  : ,hp,tv,tvp,ep,clw,m,sig
1436  : ,ment,qent,uent,vent,sij,elij,ments,qents,traent)
1437  implicit none
1438 
1439 !---------------------------------------------------------------------
1440 ! a faire:
1441 ! - changer rr(il,1) -> qnk(il)
1442 ! - vectorisation de la partie normalisation des flux (do 789...)
1443 !---------------------------------------------------------------------
1444 
1445 #include "cvthermo.h"
1446 #include "cv30param.h"
1447 
1448 c inputs:
1449  integer ncum, nd, na, ntra, nloc
1450  integer icb(nloc), inb(nloc), nk(nloc)
1451  real sig(nloc,nd)
1452  real qnk(nloc)
1453  real ph(nloc,nd+1)
1454  real t(nloc,nd), rr(nloc,nd), rs(nloc,nd)
1455  real u(nloc,nd), v(nloc,nd)
1456  real tra(nloc,nd,ntra) ! input of convect3
1457  real lv(nloc,na), h(nloc,na), hp(nloc,na)
1458  real tv(nloc,na), tvp(nloc,na), ep(nloc,na), clw(nloc,na)
1459  real m(nloc,na) ! input of convect3
1460 
1461 c outputs:
1462  real ment(nloc,na,na), qent(nloc,na,na)
1463  real uent(nloc,na,na), vent(nloc,na,na)
1464  real sij(nloc,na,na), elij(nloc,na,na)
1465  real traent(nloc,nd,nd,ntra)
1466  real ments(nloc,nd,nd), qents(nloc,nd,nd)
1467  real sigij(nloc,nd,nd)
1468 
1469 c local variables:
1470  integer i, j, k, il, im, jm
1471  integer num1, num2
1472  integer nent(nloc,na)
1473  real rti, bf2, anum, denom, dei, altem, cwat, stemp, qp
1474  real alt, smid, sjmin, sjmax, delp, delm
1475  real asij(nloc), smax(nloc), scrit(nloc)
1476  real asum(nloc,nd),bsum(nloc,nd),csum(nloc,nd)
1477  real wgh
1478  real zm(nloc,na)
1479  logical lwork(nloc)
1480 
1481 c=====================================================================
1482 c --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS
1483 c=====================================================================
1484 
1485 c ori do 360 i=1,ncum*nlp
1486  do 361 j=1,nl
1487  do 360 i=1,ncum
1488  nent(i,j)=0
1489 c in convect3, m is computed in cv3_closure
1490 c ori m(i,1)=0.0
1491  360 continue
1492  361 continue
1493 
1494 c ori do 400 k=1,nlp
1495 c ori do 390 j=1,nlp
1496  do 400 j=1,nl
1497  do 390 k=1,nl
1498  do 385 i=1,ncum
1499  qent(i,k,j)=rr(i,j)
1500  uent(i,k,j)=u(i,j)
1501  vent(i,k,j)=v(i,j)
1502  elij(i,k,j)=0.0
1503 cym ment(i,k,j)=0.0
1504 cym sij(i,k,j)=0.0
1505  385 continue
1506  390 continue
1507  400 continue
1508 
1509 cym
1510  ment(1:ncum,1:nd,1:nd)=0.0
1511  sij(1:ncum,1:nd,1:nd)=0.0
1512 
1513 c do k=1,ntra
1514 c do j=1,nd ! instead nlp
1515 c do i=1,nd ! instead nlp
1516 c do il=1,ncum
1517 c traent(il,i,j,k)=tra(il,j,k)
1518 c enddo
1519 c enddo
1520 c enddo
1521 c enddo
1522  zm(:,:)=0.
1523 
1524 c=====================================================================
1525 c --- CALCULATE ENTRAINED AIR MASS FLUX (ment), TOTAL WATER MIXING
1526 c --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
1527 c --- FRACTION (sij)
1528 c=====================================================================
1529 
1530  do 750 i=minorig+1, nl
1531 
1532  do 710 j=minorig,nl
1533  do 700 il=1,ncum
1534  if( (i.ge.icb(il)).and.(i.le.inb(il)).and.
1535  : (j.ge.(icb(il)-1)).and.(j.le.inb(il)))then
1536 
1537  rti=rr(il,1)-ep(il,i)*clw(il,i)
1538  bf2=1.+lv(il,j)*lv(il,j)*rs(il,j)/(rrv*t(il,j)*t(il,j)*cpd)
1539  anum=h(il,j)-hp(il,i)+(cpv-cpd)*t(il,j)*(rti-rr(il,j))
1540  denom=h(il,i)-hp(il,i)+(cpd-cpv)*(rr(il,i)-rti)*t(il,j)
1541  dei=denom
1542  if(abs(dei).lt.0.01)dei=0.01
1543  sij(il,i,j)=anum/dei
1544  sij(il,i,i)=1.0
1545  altem=sij(il,i,j)*rr(il,i)+(1.-sij(il,i,j))*rti-rs(il,j)
1546  altem=altem/bf2
1547  cwat=clw(il,j)*(1.-ep(il,j))
1548  stemp=sij(il,i,j)
1549  if((stemp.lt.0.0.or.stemp.gt.1.0.or.altem.gt.cwat)
1550  : .and.j.gt.i)then
1551  anum=anum-lv(il,j)*(rti-rs(il,j)-cwat*bf2)
1552  denom=denom+lv(il,j)*(rr(il,i)-rti)
1553  if(abs(denom).lt.0.01)denom=0.01
1554  sij(il,i,j)=anum/denom
1555  altem=sij(il,i,j)*rr(il,i)+(1.-sij(il,i,j))*rti-rs(il,j)
1556  altem=altem-(bf2-1.)*cwat
1557  end if
1558  if(sij(il,i,j).gt.0.0.and.sij(il,i,j).lt.0.95)then
1559  qent(il,i,j)=sij(il,i,j)*rr(il,i)+(1.-sij(il,i,j))*rti
1560  uent(il,i,j)=sij(il,i,j)*u(il,i)+(1.-sij(il,i,j))*u(il,nk(il))
1561  vent(il,i,j)=sij(il,i,j)*v(il,i)+(1.-sij(il,i,j))*v(il,nk(il))
1562 c!!! do k=1,ntra
1563 c!!! traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k)
1564 c!!! : +(1.-sij(il,i,j))*tra(il,nk(il),k)
1565 c!!! end do
1566  elij(il,i,j)=altem
1567  elij(il,i,j)=amax1(0.0,elij(il,i,j))
1568  ment(il,i,j)=m(il,i)/(1.-sij(il,i,j))
1569  nent(il,i)=nent(il,i)+1
1570  end if
1571  sij(il,i,j)=amax1(0.0,sij(il,i,j))
1572  sij(il,i,j)=amin1(1.0,sij(il,i,j))
1573  endif ! new
1574  700 continue
1575  710 continue
1576 
1577 c do k=1,ntra
1578 c do j=minorig,nl
1579 c do il=1,ncum
1580 c if( (i.ge.icb(il)).and.(i.le.inb(il)).and.
1581 c : (j.ge.(icb(il)-1)).and.(j.le.inb(il)))then
1582 c traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k)
1583 c : +(1.-sij(il,i,j))*tra(il,nk(il),k)
1584 c endif
1585 c enddo
1586 c enddo
1587 c enddo
1588 
1589 c
1590 c *** if no air can entrain at level i assume that updraft detrains ***
1591 c *** at that level and calculate detrained air flux and properties ***
1592 c
1593 
1594 c@ do 170 i=icb(il),inb(il)
1595 
1596  do 740 il=1,ncum
1597  if ((i.ge.icb(il)).and.(i.le.inb(il)).and.(nent(il,i).eq.0)) then
1598 c@ if(nent(il,i).eq.0)then
1599  ment(il,i,i)=m(il,i)
1600  qent(il,i,i)=rr(il,nk(il))-ep(il,i)*clw(il,i)
1601  uent(il,i,i)=u(il,nk(il))
1602  vent(il,i,i)=v(il,nk(il))
1603  elij(il,i,i)=clw(il,i)
1604 cMAF sij(il,i,i)=1.0
1605  sij(il,i,i)=0.0
1606  end if
1607  740 continue
1608  750 continue
1609 
1610 c do j=1,ntra
1611 c do i=minorig+1,nl
1612 c do il=1,ncum
1613 c if (i.ge.icb(il) .and. i.le.inb(il) .and. nent(il,i).eq.0) then
1614 c traent(il,i,i,j)=tra(il,nk(il),j)
1615 c endif
1616 c enddo
1617 c enddo
1618 c enddo
1619 
1620  do 100 j=minorig,nl
1621  do 101 i=minorig,nl
1622  do 102 il=1,ncum
1623  if ((j.ge.(icb(il)-1)).and.(j.le.inb(il))
1624  : .and.(i.ge.icb(il)).and.(i.le.inb(il)))then
1625  sigij(il,i,j)=sij(il,i,j)
1626  endif
1627  102 continue
1628  101 continue
1629  100 continue
1630 c@ enddo
1631 
1632 c@170 continue
1633 
1634 c=====================================================================
1635 c --- NORMALIZE ENTRAINED AIR MASS FLUXES
1636 c --- TO REPRESENT EQUAL PROBABILITIES OF MIXING
1637 c=====================================================================
1638 
1639 cym call zilch(asum,ncum*nd)
1640 cym call zilch(bsum,ncum*nd)
1641 cym call zilch(csum,ncum*nd)
1642  call zilch(asum,nloc*nd)
1643  call zilch(csum,nloc*nd)
1644  call zilch(csum,nloc*nd)
1645 
1646  do il=1,ncum
1647  lwork(il) = .false.
1648  enddo
1649 
1650  DO 789 i=minorig+1,nl
1651 
1652  num1=0
1653  do il=1,ncum
1654  if ( i.ge.icb(il) .and. i.le.inb(il) ) num1=num1+1
1655  enddo
1656  if (num1.le.0) goto 789
1657 
1658 
1659  do 781 il=1,ncum
1660  if ( i.ge.icb(il) .and. i.le.inb(il) ) then
1661  lwork(il)=(nent(il,i).ne.0)
1662  qp=rr(il,1)-ep(il,i)*clw(il,i)
1663  anum=h(il,i)-hp(il,i)-lv(il,i)*(qp-rs(il,i))
1664  : +(cpv-cpd)*t(il,i)*(qp-rr(il,i))
1665  denom=h(il,i)-hp(il,i)+lv(il,i)*(rr(il,i)-qp)
1666  : +(cpd-cpv)*t(il,i)*(rr(il,i)-qp)
1667  if(abs(denom).lt.0.01)denom=0.01
1668  scrit(il)=anum/denom
1669  alt=qp-rs(il,i)+scrit(il)*(rr(il,i)-qp)
1670  if(scrit(il).le.0.0.or.alt.le.0.0)scrit(il)=1.0
1671  smax(il)=0.0
1672  asij(il)=0.0
1673  endif
1674 781 continue
1675 
1676  do 175 j=nl,minorig,-1
1677 
1678  num2=0
1679  do il=1,ncum
1680  if ( i.ge.icb(il) .and. i.le.inb(il) .and.
1681  : j.ge.(icb(il)-1) .and. j.le.inb(il)
1682  : .and. lwork(il) ) num2=num2+1
1683  enddo
1684  if (num2.le.0) goto 175
1685 
1686  do 782 il=1,ncum
1687  if ( i.ge.icb(il) .and. i.le.inb(il) .and.
1688  : j.ge.(icb(il)-1) .and. j.le.inb(il)
1689  : .and. lwork(il) ) then
1690 
1691  if(sij(il,i,j).gt.1.0e-16.and.sij(il,i,j).lt.0.95)then
1692  wgh=1.0
1693  if(j.gt.i)then
1694  sjmax=amax1(sij(il,i,j+1),smax(il))
1695  sjmax=amin1(sjmax,scrit(il))
1696  smax(il)=amax1(sij(il,i,j),smax(il))
1697  sjmin=amax1(sij(il,i,j-1),smax(il))
1698  sjmin=amin1(sjmin,scrit(il))
1699  if(sij(il,i,j).lt.(smax(il)-1.0e-16))wgh=0.0
1700  smid=amin1(sij(il,i,j),scrit(il))
1701  else
1702  sjmax=amax1(sij(il,i,j+1),scrit(il))
1703  smid=amax1(sij(il,i,j),scrit(il))
1704  sjmin=0.0
1705  if(j.gt.1)sjmin=sij(il,i,j-1)
1706  sjmin=amax1(sjmin,scrit(il))
1707  endif
1708  delp=abs(sjmax-smid)
1709  delm=abs(sjmin-smid)
1710  asij(il)=asij(il)+wgh*(delp+delm)
1711  ment(il,i,j)=ment(il,i,j)*(delp+delm)*wgh
1712  endif
1713  endif
1714 782 continue
1715 
1716 175 continue
1717 
1718  do il=1,ncum
1719  if (i.ge.icb(il).and.i.le.inb(il).and.lwork(il)) then
1720  asij(il)=amax1(1.0e-16,asij(il))
1721  asij(il)=1.0/asij(il)
1722  asum(il,i)=0.0
1723  bsum(il,i)=0.0
1724  csum(il,i)=0.0
1725  endif
1726  enddo
1727 
1728  do 180 j=minorig,nl
1729  do il=1,ncum
1730  if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il)
1731  : .and. j.ge.(icb(il)-1) .and. j.le.inb(il) ) then
1732  ment(il,i,j)=ment(il,i,j)*asij(il)
1733  endif
1734  enddo
1735 180 continue
1736 
1737  do 190 j=minorig,nl
1738  do il=1,ncum
1739  if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il)
1740  : .and. j.ge.(icb(il)-1) .and. j.le.inb(il) ) then
1741  asum(il,i)=asum(il,i)+ment(il,i,j)
1742  ment(il,i,j)=ment(il,i,j)*sig(il,j)
1743  bsum(il,i)=bsum(il,i)+ment(il,i,j)
1744  endif
1745  enddo
1746 190 continue
1747 
1748  do il=1,ncum
1749  if (i.ge.icb(il).and.i.le.inb(il).and.lwork(il)) then
1750  bsum(il,i)=amax1(bsum(il,i),1.0e-16)
1751  bsum(il,i)=1.0/bsum(il,i)
1752  endif
1753  enddo
1754 
1755  do 195 j=minorig,nl
1756  do il=1,ncum
1757  if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il)
1758  : .and. j.ge.(icb(il)-1) .and. j.le.inb(il) ) then
1759  ment(il,i,j)=ment(il,i,j)*asum(il,i)*bsum(il,i)
1760  endif
1761  enddo
1762 195 continue
1763 
1764  do 197 j=minorig,nl
1765  do il=1,ncum
1766  if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il)
1767  : .and. j.ge.(icb(il)-1) .and. j.le.inb(il) ) then
1768  csum(il,i)=csum(il,i)+ment(il,i,j)
1769  endif
1770  enddo
1771 197 continue
1772 
1773  do il=1,ncum
1774  if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il)
1775  : .and. csum(il,i).lt.m(il,i) ) then
1776  nent(il,i)=0
1777  ment(il,i,i)=m(il,i)
1778  qent(il,i,i)=rr(il,1)-ep(il,i)*clw(il,i)
1779  uent(il,i,i)=u(il,nk(il))
1780  vent(il,i,i)=v(il,nk(il))
1781  elij(il,i,i)=clw(il,i)
1782 cMAF sij(il,i,i)=1.0
1783  sij(il,i,i)=0.0
1784  endif
1785  enddo ! il
1786 
1787 c do j=1,ntra
1788 c do il=1,ncum
1789 c if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il)
1790 c : .and. csum(il,i).lt.m(il,i) ) then
1791 c traent(il,i,i,j)=tra(il,nk(il),j)
1792 c endif
1793 c enddo
1794 c enddo
1795 789 continue
1796 c
1797 c MAF: renormalisation de MENT
1798  do jm=1,nd
1799  do im=1,nd
1800  do il=1,ncum
1801  zm(il,im)=zm(il,im)+(1.-sij(il,im,jm))*ment(il,im,jm)
1802  end do
1803  end do
1804  end do
1805 c
1806  do jm=1,nd
1807  do im=1,nd
1808  do il=1,ncum
1809  if(zm(il,im).ne.0.) then
1810  ment(il,im,jm)=ment(il,im,jm)*m(il,im)/zm(il,im)
1811  endif
1812  end do
1813  end do
1814  end do
1815 c
1816  do jm=1,nd
1817  do im=1,nd
1818  do 999 il=1,ncum
1819  qents(il,im,jm)=qent(il,im,jm)
1820  ments(il,im,jm)=ment(il,im,jm)
1821 999 continue
1822  enddo
1823  enddo
1824 
1825  return
1826  end
1827 
1828 
1829  SUBROUTINE cv30_unsat(nloc,ncum,nd,na,ntra,icb,inb
1830  : ,t,rr,rs,gz,u,v,tra,p,ph
1831  : ,th,tv,lv,cpn,ep,sigp,clw
1832  : ,m,ment,elij,delt,plcl
1833  : ,mp,rp,up,vp,trap,wt,water,evap,b ! RomP-jyg
1834  : ,wdtraina,wdtrainm) ! 26/08/10 RomP-jyg
1835  implicit none
1836 
1837 
1838 #include "cvthermo.h"
1839 #include "cv30param.h"
1840 #include "cvflag.h"
1841 
1842 c inputs:
1843  integer ncum, nd, na, ntra, nloc
1844  integer icb(nloc), inb(nloc)
1845  real delt, plcl(nloc)
1846  real t(nloc,nd), rr(nloc,nd), rs(nloc,nd)
1847  real u(nloc,nd), v(nloc,nd)
1848  real tra(nloc,nd,ntra)
1849  real p(nloc,nd), ph(nloc,nd+1)
1850  real th(nloc,na), gz(nloc,na)
1851  real lv(nloc,na), ep(nloc,na), sigp(nloc,na), clw(nloc,na)
1852  real cpn(nloc,na), tv(nloc,na)
1853  real m(nloc,na), ment(nloc,na,na), elij(nloc,na,na)
1854 
1855 c outputs:
1856  real mp(nloc,na), rp(nloc,na), up(nloc,na), vp(nloc,na)
1857  real water(nloc,na), evap(nloc,na), wt(nloc,na)
1858  real trap(nloc,na,ntra)
1859  real b(nloc,na)
1860 ! 25/08/10 - RomP---- ajout des masses precipitantes ejectees
1861 ! lascendance adiabatique et des flux melanges Pa et Pm.
1862 ! Distinction des wdtrain
1863 ! Pa = wdtrainA Pm = wdtrainM
1864  real wdtraina(nloc,na), wdtrainm(nloc,na)
1865 
1866 c local variables
1867  integer i,j,k,il,num1
1868  real tinv, delti
1869  real awat, afac, afac1, afac2, bfac
1870  real pr1, pr2, sigt, b6, c6, revap, tevap, delth
1871  real amfac, amp2, xf, tf, fac2, ur, sru, fac, d, af, bf
1872  real ampmax
1873  real lvcp(nloc,na)
1874  real wdtrain(nloc)
1875  logical lwork(nloc)
1876 
1877 
1878 c------------------------------------------------------
1879 
1880  delti = 1./delt
1881  tinv=1./3.
1882 
1883  mp(:,:)=0.
1884 
1885  do i=1,nl
1886  do il=1,ncum
1887  mp(il,i)=0.0
1888  rp(il,i)=rr(il,i)
1889  up(il,i)=u(il,i)
1890  vp(il,i)=v(il,i)
1891  wt(il,i)=0.001
1892  water(il,i)=0.0
1893  evap(il,i)=0.0
1894  b(il,i)=0.0
1895  lvcp(il,i)=lv(il,i)/cpn(il,i)
1896  enddo
1897  enddo
1898 
1899 c do k=1,ntra
1900 c do i=1,nd
1901 c do il=1,ncum
1902 c trap(il,i,k)=tra(il,i,k)
1903 c enddo
1904 c enddo
1905 c enddo
1906 !! RomP >>>
1907  do i=1,nd
1908  do il=1,ncum
1909  wdtraina(il,i)=0.0
1910  wdtrainm(il,i)=0.0
1911  enddo
1912  enddo
1913 !! RomP <<<
1914 c
1915 c *** check whether ep(inb)=0, if so, skip precipitating ***
1916 c *** downdraft calculation ***
1917 c
1918 
1919  do il=1,ncum
1920  lwork(il)=.true.
1921  if(ep(il,inb(il)).lt.0.0001)lwork(il)=.false.
1922  enddo
1923 
1924  call zilch(wdtrain,ncum)
1925 
1926  DO 400 i=nl+1,1,-1
1927 
1928  num1=0
1929  do il=1,ncum
1930  if ( i.le.inb(il) .and. lwork(il) ) num1=num1+1
1931  enddo
1932  if (num1.le.0) goto 400
1933 
1934 c
1935 c *** integrate liquid water equation to find condensed water ***
1936 c *** and condensed water flux ***
1937 c
1938 
1939 c
1940 c *** begin downdraft loop ***
1941 c
1942 
1943 c
1944 c *** calculate detrained precipitation ***
1945 c
1946  do il=1,ncum
1947  if (i.le.inb(il) .and. lwork(il)) then
1948  if (cvflag_grav) then
1949  wdtrain(il)=grav*ep(il,i)*m(il,i)*clw(il,i)
1950  wdtraina(il,i) = wdtrain(il)/grav ! Pa 26/08/10 RomP
1951  else
1952  wdtrain(il)=10.0*ep(il,i)*m(il,i)*clw(il,i)
1953  wdtraina(il,i) = wdtrain(il)/10. ! Pa 26/08/10 RomP
1954  endif
1955  endif
1956  enddo
1957 
1958  if(i.gt.1)then
1959 
1960  do 320 j=1,i-1
1961  do il=1,ncum
1962  if (i.le.inb(il) .and. lwork(il)) then
1963  awat=elij(il,j,i)-(1.-ep(il,i))*clw(il,i)
1964  awat=amax1(awat,0.0)
1965  if (cvflag_grav) then
1966  wdtrain(il)=wdtrain(il)+grav*awat*ment(il,j,i)
1967  else
1968  wdtrain(il)=wdtrain(il)+10.0*awat*ment(il,j,i)
1969  endif
1970  endif
1971  enddo
1972 320 continue
1973  do il=1,ncum
1974  if (cvflag_grav) then
1975  wdtrainm(il,i) = wdtrain(il)/grav-wdtraina(il,i) ! Pm 26/08/10 RomP
1976  else
1977  wdtrainm(il,i) = wdtrain(il)/10.-wdtraina(il,i) ! Pm 26/08/10 RomP
1978  endif
1979  enddo
1980 
1981  endif
1982 
1983 c
1984 c *** find rain water and evaporation using provisional ***
1985 c *** estimates of rp(i)and rp(i-1) ***
1986 c
1987 
1988  do 999 il=1,ncum
1989 
1990  if (i.le.inb(il) .and. lwork(il)) then
1991 
1992  wt(il,i)=45.0
1993 
1994  if(i.lt.inb(il))then
1995  rp(il,i)=rp(il,i+1)
1996  : +(cpd*(t(il,i+1)-t(il,i))+gz(il,i+1)-gz(il,i))/lv(il,i)
1997  rp(il,i)=0.5*(rp(il,i)+rr(il,i))
1998  endif
1999  rp(il,i)=amax1(rp(il,i),0.0)
2000  rp(il,i)=amin1(rp(il,i),rs(il,i))
2001  rp(il,inb(il))=rr(il,inb(il))
2002 
2003  if(i.eq.1)then
2004  afac=p(il,1)*(rs(il,1)-rp(il,1))/(1.0e4+2000.0*p(il,1)*rs(il,1))
2005  else
2006  rp(il,i-1)=rp(il,i)
2007  : +(cpd*(t(il,i)-t(il,i-1))+gz(il,i)-gz(il,i-1))/lv(il,i)
2008  rp(il,i-1)=0.5*(rp(il,i-1)+rr(il,i-1))
2009  rp(il,i-1)=amin1(rp(il,i-1),rs(il,i-1))
2010  rp(il,i-1)=amax1(rp(il,i-1),0.0)
2011  afac1=p(il,i)*(rs(il,i)-rp(il,i))/(1.0e4+2000.0*p(il,i)*rs(il,i))
2012  afac2=p(il,i-1)*(rs(il,i-1)-rp(il,i-1))
2013  : /(1.0e4+2000.0*p(il,i-1)*rs(il,i-1))
2014  afac=0.5*(afac1+afac2)
2015  endif
2016  if(i.eq.inb(il))afac=0.0
2017  afac=amax1(afac,0.0)
2018  bfac=1./(sigd*wt(il,i))
2019 c
2020 cjyg1
2021 ccc sigt=1.0
2022 ccc if(i.ge.icb)sigt=sigp(i)
2023 c prise en compte de la variation progressive de sigt dans
2024 c les couches icb et icb-1:
2025 c pour plcl<ph(i+1), pr1=0 & pr2=1
2026 c pour plcl>ph(i), pr1=1 & pr2=0
2027 c pour ph(i+1)<plcl<ph(i), pr1 est la proportion a cheval
2028 c sur le nuage, et pr2 est la proportion sous la base du
2029 c nuage.
2030  pr1=(plcl(il)-ph(il,i+1))/(ph(il,i)-ph(il,i+1))
2031  pr1=max(0.,min(1.,pr1))
2032  pr2=(ph(il,i)-plcl(il))/(ph(il,i)-ph(il,i+1))
2033  pr2=max(0.,min(1.,pr2))
2034  sigt=sigp(il,i)*pr1+pr2
2035 cjyg2
2036 c
2037  b6=bfac*50.*sigd*(ph(il,i)-ph(il,i+1))*sigt*afac
2038  c6=water(il,i+1)+bfac*wdtrain(il)
2039  : -50.*sigd*bfac*(ph(il,i)-ph(il,i+1))*evap(il,i+1)
2040  if(c6.gt.0.0)then
2041  revap=0.5*(-b6+sqrt(b6*b6+4.*c6))
2042  evap(il,i)=sigt*afac*revap
2043  water(il,i)=revap*revap
2044  else
2045  evap(il,i)=-evap(il,i+1)
2046  : +0.02*(wdtrain(il)+sigd*wt(il,i)*water(il,i+1))
2047  : /(sigd*(ph(il,i)-ph(il,i+1)))
2048  end if
2049 c
2050 c *** calculate precipitating downdraft mass flux under ***
2051 c *** hydrostatic approximation ***
2052 c
2053  if (i.ne.1) then
2054 
2055  tevap=amax1(0.0,evap(il,i))
2056  delth=amax1(0.001,(th(il,i)-th(il,i-1)))
2057  if (cvflag_grav) then
2058  mp(il,i)=100.*ginv*lvcp(il,i)*sigd*tevap
2059  : *(p(il,i-1)-p(il,i))/delth
2060  else
2061  mp(il,i)=10.*lvcp(il,i)*sigd*tevap*(p(il,i-1)-p(il,i))/delth
2062  endif
2063 c
2064 c *** if hydrostatic assumption fails, ***
2065 c *** solve cubic difference equation for downdraft theta ***
2066 c *** and mass flux from two simultaneous differential eqns ***
2067 c
2068  amfac=sigd*sigd*70.0*ph(il,i)*(p(il,i-1)-p(il,i))
2069  : *(th(il,i)-th(il,i-1))/(tv(il,i)*th(il,i))
2070  amp2=abs(mp(il,i+1)*mp(il,i+1)-mp(il,i)*mp(il,i))
2071  if(amp2.gt.(0.1*amfac))then
2072  xf=100.0*sigd*sigd*sigd*(ph(il,i)-ph(il,i+1))
2073  tf=b(il,i)-5.0*(th(il,i)-th(il,i-1))*t(il,i)
2074  : /(lvcp(il,i)*sigd*th(il,i))
2075  af=xf*tf+mp(il,i+1)*mp(il,i+1)*tinv
2076  bf=2.*(tinv*mp(il,i+1))**3+tinv*mp(il,i+1)*xf*tf
2077  : +50.*(p(il,i-1)-p(il,i))*xf*tevap
2078  fac2=1.0
2079  if(bf.lt.0.0)fac2=-1.0
2080  bf=abs(bf)
2081  ur=0.25*bf*bf-af*af*af*tinv*tinv*tinv
2082  if(ur.ge.0.0)then
2083  sru=sqrt(ur)
2084  fac=1.0
2085  if((0.5*bf-sru).lt.0.0)fac=-1.0
2086  mp(il,i)=mp(il,i+1)*tinv+(0.5*bf+sru)**tinv
2087  : +fac*(abs(0.5*bf-sru))**tinv
2088  else
2089  d=atan(2.*sqrt(-ur)/(bf+1.0e-28))
2090  if(fac2.lt.0.0)d=3.14159-d
2091  mp(il,i)=mp(il,i+1)*tinv+2.*sqrt(af*tinv)*cos(d*tinv)
2092  endif
2093  mp(il,i)=amax1(0.0,mp(il,i))
2094 
2095  if (cvflag_grav) then
2096 Cjyg : il y a vraisemblablement une erreur dans la ligne 2 suivante:
2097 C il faut diviser par (mp(il,i)*sigd*grav) et non par (mp(il,i)+sigd*0.1).
2098 C Et il faut bien revoir les facteurs 100.
2099  b(il,i-1)=b(il,i)+100.0*(p(il,i-1)-p(il,i))*tevap
2100  2 /(mp(il,i)+sigd*0.1)
2101  3 -10.0*(th(il,i)-th(il,i-1))*t(il,i)/(lvcp(il,i)*sigd*th(il,i))
2102  else
2103  b(il,i-1)=b(il,i)+100.0*(p(il,i-1)-p(il,i))*tevap
2104  2 /(mp(il,i)+sigd*0.1)
2105  3 -10.0*(th(il,i)-th(il,i-1))*t(il,i)/(lvcp(il,i)*sigd*th(il,i))
2106  endif
2107  b(il,i-1)=amax1(b(il,i-1),0.0)
2108  endif
2109 c
2110 c *** limit magnitude of mp(i) to meet cfl condition ***
2111 c
2112  ampmax=2.0*(ph(il,i)-ph(il,i+1))*delti
2113  amp2=2.0*(ph(il,i-1)-ph(il,i))*delti
2114  ampmax=amin1(ampmax,amp2)
2115  mp(il,i)=amin1(mp(il,i),ampmax)
2116 c
2117 c *** force mp to decrease linearly to zero ***
2118 c *** between cloud base and the surface ***
2119 c
2120  if(p(il,i).gt.p(il,icb(il)))then
2121  mp(il,i)=mp(il,icb(il))*(p(il,1)-p(il,i))/(p(il,1)-p(il,icb(il)))
2122  endif
2123 
2124 360 continue
2125  endif ! i.eq.1
2126 c
2127 c *** find mixing ratio of precipitating downdraft ***
2128 c
2129 
2130  if (i.ne.inb(il)) then
2131 
2132  rp(il,i)=rr(il,i)
2133 
2134  if(mp(il,i).gt.mp(il,i+1))then
2135 
2136  if (cvflag_grav) then
2137  rp(il,i)=rp(il,i+1)*mp(il,i+1)+rr(il,i)*(mp(il,i)-mp(il,i+1))
2138  : +100.*ginv*0.5*sigd*(ph(il,i)-ph(il,i+1))
2139  : *(evap(il,i+1)+evap(il,i))
2140  else
2141  rp(il,i)=rp(il,i+1)*mp(il,i+1)+rr(il,i)*(mp(il,i)-mp(il,i+1))
2142  : +5.*sigd*(ph(il,i)-ph(il,i+1))
2143  : *(evap(il,i+1)+evap(il,i))
2144  endif
2145  rp(il,i)=rp(il,i)/mp(il,i)
2146  up(il,i)=up(il,i+1)*mp(il,i+1)+u(il,i)*(mp(il,i)-mp(il,i+1))
2147  up(il,i)=up(il,i)/mp(il,i)
2148  vp(il,i)=vp(il,i+1)*mp(il,i+1)+v(il,i)*(mp(il,i)-mp(il,i+1))
2149  vp(il,i)=vp(il,i)/mp(il,i)
2150 
2151 c do j=1,ntra
2152 c trap(il,i,j)=trap(il,i+1,j)*mp(il,i+1)
2153 ctestmaf : +trap(il,i,j)*(mp(il,i)-mp(il,i+1))
2154 c : +tra(il,i,j)*(mp(il,i)-mp(il,i+1))
2155 c trap(il,i,j)=trap(il,i,j)/mp(il,i)
2156 c end do
2157 
2158  else
2159 
2160  if(mp(il,i+1).gt.1.0e-16)then
2161  if (cvflag_grav) then
2162  rp(il,i)=rp(il,i+1)
2163  : +100.*ginv*0.5*sigd*(ph(il,i)-ph(il,i+1))
2164  : *(evap(il,i+1)+evap(il,i))/mp(il,i+1)
2165  else
2166  rp(il,i)=rp(il,i+1)
2167  : +5.*sigd*(ph(il,i)-ph(il,i+1))
2168  : *(evap(il,i+1)+evap(il,i))/mp(il,i+1)
2169  endif
2170  up(il,i)=up(il,i+1)
2171  vp(il,i)=vp(il,i+1)
2172 
2173 c do j=1,ntra
2174 c trap(il,i,j)=trap(il,i+1,j)
2175 c end do
2176 
2177  endif
2178  endif
2179  rp(il,i)=amin1(rp(il,i),rs(il,i))
2180  rp(il,i)=amax1(rp(il,i),0.0)
2181 
2182  endif
2183  endif
2184 999 continue
2185 
2186 400 continue
2187 
2188  return
2189  end
2190 
2191  SUBROUTINE cv30_yield(nloc,ncum,nd,na,ntra
2192  : ,icb,inb,delt
2193  : ,t,rr,u,v,tra,gz,p,ph,h,hp,lv,cpn,th
2194  : ,ep,clw,m,tp,mp,rp,up,vp,trap
2195  : ,wt,water,evap,b
2196  : ,ment,qent,uent,vent,nent,elij,traent,sig
2197  : ,tv,tvp
2198  : ,iflag,precip,vprecip,ft,fr,fu,fv,ftra
2199  : ,upwd,dnwd,dnwd0,ma,mike,tls,tps,qcondc,wd)
2200  implicit none
2201 
2202 #include "cvthermo.h"
2203 #include "cv30param.h"
2204 #include "cvflag.h"
2205 #include "conema3.h"
2206 
2207 c inputs:
2208  integer ncum,nd,na,ntra,nloc
2209  integer icb(nloc), inb(nloc)
2210  real delt
2211  real t(nloc,nd), rr(nloc,nd), u(nloc,nd), v(nloc,nd)
2212  real tra(nloc,nd,ntra), sig(nloc,nd)
2213  real gz(nloc,na), ph(nloc,nd+1), h(nloc,na), hp(nloc,na)
2214  real th(nloc,na), p(nloc,nd), tp(nloc,na)
2215  real lv(nloc,na), cpn(nloc,na), ep(nloc,na), clw(nloc,na)
2216  real m(nloc,na), mp(nloc,na), rp(nloc,na), up(nloc,na)
2217  real vp(nloc,na), wt(nloc,nd), trap(nloc,nd,ntra)
2218  real water(nloc,na), evap(nloc,na), b(nloc,na)
2219  real ment(nloc,na,na), qent(nloc,na,na), uent(nloc,na,na)
2220 cym real vent(nloc,na,na), nent(nloc,na), elij(nloc,na,na)
2221  real vent(nloc,na,na), elij(nloc,na,na)
2222  integer nent(nloc,na)
2223  real traent(nloc,na,na,ntra)
2224  real tv(nloc,nd), tvp(nloc,nd)
2225 
2226 c input/output:
2227  integer iflag(nloc)
2228 
2229 c outputs:
2230  real precip(nloc)
2231  real vprecip(nloc,nd+1)
2232  real ft(nloc,nd), fr(nloc,nd), fu(nloc,nd), fv(nloc,nd)
2233  real ftra(nloc,nd,ntra)
2234  real upwd(nloc,nd), dnwd(nloc,nd), ma(nloc,nd)
2235  real dnwd0(nloc,nd), mike(nloc,nd)
2236  real tls(nloc,nd), tps(nloc,nd)
2237  real qcondc(nloc,nd) ! cld
2238  real wd(nloc) ! gust
2239 
2240 c local variables:
2241  integer i,k,il,n,j,num1
2242  real rat, awat, delti
2243  real ax, bx, cx, dx, ex
2244  real cpinv, rdcp, dpinv
2245  real lvcp(nloc,na), mke(nloc,na)
2246  real am(nloc), work(nloc), ad(nloc), amp1(nloc)
2247 c!! real up1(nloc), dn1(nloc)
2248  real up1(nloc,nd,nd), dn1(nloc,nd,nd)
2249  real asum(nloc), bsum(nloc), csum(nloc), dsum(nloc)
2250  real qcond(nloc,nd), nqcond(nloc,nd), wa(nloc,nd) ! cld
2251  real siga(nloc,nd), sax(nloc,nd), mac(nloc,nd) ! cld
2252 
2253 
2254 c-------------------------------------------------------------
2255 
2256 c initialization:
2257 
2258  delti = 1.0/delt
2259 
2260  do il=1,ncum
2261  precip(il)=0.0
2262  wd(il)=0.0 ! gust
2263  vprecip(il,nd+1)=0.
2264  enddo
2265 
2266  do i=1,nd
2267  do il=1,ncum
2268  vprecip(il,i)=0.0
2269  ft(il,i)=0.0
2270  fr(il,i)=0.0
2271  fu(il,i)=0.0
2272  fv(il,i)=0.0
2273  qcondc(il,i)=0.0 ! cld
2274  qcond(il,i)=0.0 ! cld
2275  nqcond(il,i)=0.0 ! cld
2276  enddo
2277  enddo
2278 
2279 c do j=1,ntra
2280 c do i=1,nd
2281 c do il=1,ncum
2282 c ftra(il,i,j)=0.0
2283 c enddo
2284 c enddo
2285 c enddo
2286 
2287  do i=1,nl
2288  do il=1,ncum
2289  lvcp(il,i)=lv(il,i)/cpn(il,i)
2290  enddo
2291  enddo
2292 
2293 
2294 c
2295 c *** calculate surface precipitation in mm/day ***
2296 c
2297  do il=1,ncum
2298  if(ep(il,inb(il)).ge.0.0001)then
2299  if (cvflag_grav) then
2300  precip(il)=wt(il,1)*sigd*water(il,1)*86400.*1000./(rowl*grav)
2301  else
2302  precip(il)=wt(il,1)*sigd*water(il,1)*8640.
2303  endif
2304  endif
2305  enddo
2306 
2307 C *** CALCULATE VERTICAL PROFILE OF PRECIPITATIONs IN kg/m2/s ===
2308 C
2309 c MAF rajout pour lessivage
2310  do k=1,nl
2311  do il=1,ncum
2312  if (k.le.inb(il)) then
2313  if (cvflag_grav) then
2314  vprecip(il,k) = wt(il,k)*sigd*water(il,k)/grav
2315  else
2316  vprecip(il,k) = wt(il,k)*sigd*water(il,k)/10.
2317  endif
2318  endif
2319  end do
2320  end do
2321 C
2322 c
2323 c *** Calculate downdraft velocity scale ***
2324 c *** NE PAS UTILISER POUR L'INSTANT ***
2325 c
2326 c! do il=1,ncum
2327 c! wd(il)=betad*abs(mp(il,icb(il)))*0.01*rrd*t(il,icb(il))
2328 c! : /(sigd*p(il,icb(il)))
2329 c! enddo
2330 
2331 c
2332 c *** calculate tendencies of lowest level potential temperature ***
2333 c *** and mixing ratio ***
2334 c
2335  do il=1,ncum
2336  work(il)=1.0/(ph(il,1)-ph(il,2))
2337  am(il)=0.0
2338  enddo
2339 
2340  do k=2,nl
2341  do il=1,ncum
2342  if (k.le.inb(il)) then
2343  am(il)=am(il)+m(il,k)
2344  endif
2345  enddo
2346  enddo
2347 
2348  do il=1,ncum
2349 
2350 c convect3 if((0.1*dpinv*am).ge.delti)iflag(il)=4
2351  if (cvflag_grav) then
2352  if((0.01*grav*work(il)*am(il)).ge.delti)iflag(il)=1!consist vect
2353  ft(il,1)=0.01*grav*work(il)*am(il)*(t(il,2)-t(il,1)
2354  : +(gz(il,2)-gz(il,1))/cpn(il,1))
2355  else
2356  if((0.1*work(il)*am(il)).ge.delti)iflag(il)=1 !consistency vect
2357  ft(il,1)=0.1*work(il)*am(il)*(t(il,2)-t(il,1)
2358  : +(gz(il,2)-gz(il,1))/cpn(il,1))
2359  endif
2360 
2361  ft(il,1)=ft(il,1)-0.5*lvcp(il,1)*sigd*(evap(il,1)+evap(il,2))
2362 
2363  if (cvflag_grav) then
2364  ft(il,1)=ft(il,1)-0.009*grav*sigd*mp(il,2)
2365  : *t(il,1)*b(il,1)*work(il)
2366  else
2367  ft(il,1)=ft(il,1)-0.09*sigd*mp(il,2)*t(il,1)*b(il,1)*work(il)
2368  endif
2369 
2370  ft(il,1)=ft(il,1)+0.01*sigd*wt(il,1)*(cl-cpd)*water(il,2)*(t(il,2)
2371  :-t(il,1))*work(il)/cpn(il,1)
2372 
2373  if (cvflag_grav) then
2374 Cjyg1 Correction pour mieux conserver l'eau (conformite avec CONVECT4.3)
2375 c (sb: pour l'instant, on ne fait que le chgt concernant grav, pas evap)
2376  fr(il,1)=0.01*grav*mp(il,2)*(rp(il,2)-rr(il,1))*work(il)
2377  : +sigd*0.5*(evap(il,1)+evap(il,2))
2378 c+tard : +sigd*evap(il,1)
2379 
2380  fr(il,1)=fr(il,1)+0.01*grav*am(il)*(rr(il,2)-rr(il,1))*work(il)
2381 
2382  fu(il,1)=fu(il,1)+0.01*grav*work(il)*(mp(il,2)*(up(il,2)-u(il,1))
2383  : +am(il)*(u(il,2)-u(il,1)))
2384  fv(il,1)=fv(il,1)+0.01*grav*work(il)*(mp(il,2)*(vp(il,2)-v(il,1))
2385  : +am(il)*(v(il,2)-v(il,1)))
2386  else ! cvflag_grav
2387  fr(il,1)=0.1*mp(il,2)*(rp(il,2)-rr(il,1))*work(il)
2388  : +sigd*0.5*(evap(il,1)+evap(il,2))
2389  fr(il,1)=fr(il,1)+0.1*am(il)*(rr(il,2)-rr(il,1))*work(il)
2390  fu(il,1)=fu(il,1)+0.1*work(il)*(mp(il,2)*(up(il,2)-u(il,1))
2391  : +am(il)*(u(il,2)-u(il,1)))
2392  fv(il,1)=fv(il,1)+0.1*work(il)*(mp(il,2)*(vp(il,2)-v(il,1))
2393  : +am(il)*(v(il,2)-v(il,1)))
2394  endif ! cvflag_grav
2395 
2396  enddo ! il
2397 
2398 c do j=1,ntra
2399 c do il=1,ncum
2400 c if (cvflag_grav) then
2401 c ftra(il,1,j)=ftra(il,1,j)+0.01*grav*work(il)
2402 c : *(mp(il,2)*(trap(il,2,j)-tra(il,1,j))
2403 c : +am(il)*(tra(il,2,j)-tra(il,1,j)))
2404 c else
2405 c ftra(il,1,j)=ftra(il,1,j)+0.1*work(il)
2406 c : *(mp(il,2)*(trap(il,2,j)-tra(il,1,j))
2407 c : +am(il)*(tra(il,2,j)-tra(il,1,j)))
2408 c endif
2409 c enddo
2410 c enddo
2411 
2412  do j=2,nl
2413  do il=1,ncum
2414  if (j.le.inb(il)) then
2415  if (cvflag_grav) then
2416  fr(il,1)=fr(il,1)
2417  : +0.01*grav*work(il)*ment(il,j,1)*(qent(il,j,1)-rr(il,1))
2418  fu(il,1)=fu(il,1)
2419  : +0.01*grav*work(il)*ment(il,j,1)*(uent(il,j,1)-u(il,1))
2420  fv(il,1)=fv(il,1)
2421  : +0.01*grav*work(il)*ment(il,j,1)*(vent(il,j,1)-v(il,1))
2422  else ! cvflag_grav
2423  fr(il,1)=fr(il,1)
2424  : +0.1*work(il)*ment(il,j,1)*(qent(il,j,1)-rr(il,1))
2425  fu(il,1)=fu(il,1)
2426  : +0.1*work(il)*ment(il,j,1)*(uent(il,j,1)-u(il,1))
2427  fv(il,1)=fv(il,1)
2428  : +0.1*work(il)*ment(il,j,1)*(vent(il,j,1)-v(il,1))
2429  endif ! cvflag_grav
2430  endif ! j
2431  enddo
2432  enddo
2433 
2434 c do k=1,ntra
2435 c do j=2,nl
2436 c do il=1,ncum
2437 c if (j.le.inb(il)) then
2438 
2439 c if (cvflag_grav) then
2440 c ftra(il,1,k)=ftra(il,1,k)+0.01*grav*work(il)*ment(il,j,1)
2441 c : *(traent(il,j,1,k)-tra(il,1,k))
2442 c else
2443 c ftra(il,1,k)=ftra(il,1,k)+0.1*work(il)*ment(il,j,1)
2444 c : *(traent(il,j,1,k)-tra(il,1,k))
2445 c endif
2446 
2447 c endif
2448 c enddo
2449 c enddo
2450 c enddo
2451 
2452 c
2453 c *** calculate tendencies of potential temperature and mixing ratio ***
2454 c *** at levels above the lowest level ***
2455 c
2456 c *** first find the net saturated updraft and downdraft mass fluxes ***
2457 c *** through each level ***
2458 c
2459 
2460  do 500 i=2,nl+1 ! newvecto: mettre nl au lieu nl+1?
2461 
2462  num1=0
2463  do il=1,ncum
2464  if(i.le.inb(il))num1=num1+1
2465  enddo
2466  if(num1.le.0)go to 500
2467 
2468  call zilch(amp1,ncum)
2469  call zilch(ad,ncum)
2470 
2471  do 440 k=i+1,nl+1
2472  do 441 il=1,ncum
2473  if (i.le.inb(il) .and. k.le.(inb(il)+1)) then
2474  amp1(il)=amp1(il)+m(il,k)
2475  endif
2476  441 continue
2477  440 continue
2478 
2479  do 450 k=1,i
2480  do 451 j=i+1,nl+1
2481  do 452 il=1,ncum
2482  if (i.le.inb(il) .and. j.le.(inb(il)+1)) then
2483  amp1(il)=amp1(il)+ment(il,k,j)
2484  endif
2485 452 continue
2486 451 continue
2487 450 continue
2488 
2489  do 470 k=1,i-1
2490  do 471 j=i,nl+1 ! newvecto: nl au lieu nl+1?
2491  do 472 il=1,ncum
2492  if (i.le.inb(il) .and. j.le.inb(il)) then
2493  ad(il)=ad(il)+ment(il,j,k)
2494  endif
2495 472 continue
2496 471 continue
2497 470 continue
2498 
2499  do 1350 il=1,ncum
2500  if (i.le.inb(il)) then
2501  dpinv=1.0/(ph(il,i)-ph(il,i+1))
2502  cpinv=1.0/cpn(il,i)
2503 
2504 c convect3 if((0.1*dpinv*amp1).ge.delti)iflag(il)=4
2505  if (cvflag_grav) then
2506  if((0.01*grav*dpinv*amp1(il)).ge.delti)iflag(il)=1 ! vecto
2507  else
2508  if((0.1*dpinv*amp1(il)).ge.delti)iflag(il)=1 ! vecto
2509  endif
2510 
2511  if (cvflag_grav) then
2512  ft(il,i)=0.01*grav*dpinv*(amp1(il)*(t(il,i+1)-t(il,i)
2513  : +(gz(il,i+1)-gz(il,i))*cpinv)
2514  : -ad(il)*(t(il,i)-t(il,i-1)+(gz(il,i)-gz(il,i-1))*cpinv))
2515  : -0.5*sigd*lvcp(il,i)*(evap(il,i)+evap(il,i+1))
2516  rat=cpn(il,i-1)*cpinv
2517  ft(il,i)=ft(il,i)-0.009*grav*sigd*(mp(il,i+1)*t(il,i)*b(il,i)
2518  : -mp(il,i)*t(il,i-1)*rat*b(il,i-1))*dpinv
2519  ft(il,i)=ft(il,i)+0.01*grav*dpinv*ment(il,i,i)*(hp(il,i)-h(il,i)
2520  : +t(il,i)*(cpv-cpd)*(rr(il,i)-qent(il,i,i)))*cpinv
2521  else ! cvflag_grav
2522  ft(il,i)=0.1*dpinv*(amp1(il)*(t(il,i+1)-t(il,i)
2523  : +(gz(il,i+1)-gz(il,i))*cpinv)
2524  : -ad(il)*(t(il,i)-t(il,i-1)+(gz(il,i)-gz(il,i-1))*cpinv))
2525  : -0.5*sigd*lvcp(il,i)*(evap(il,i)+evap(il,i+1))
2526  rat=cpn(il,i-1)*cpinv
2527  ft(il,i)=ft(il,i)-0.09*sigd*(mp(il,i+1)*t(il,i)*b(il,i)
2528  : -mp(il,i)*t(il,i-1)*rat*b(il,i-1))*dpinv
2529  ft(il,i)=ft(il,i)+0.1*dpinv*ment(il,i,i)*(hp(il,i)-h(il,i)
2530  : +t(il,i)*(cpv-cpd)*(rr(il,i)-qent(il,i,i)))*cpinv
2531  endif ! cvflag_grav
2532 
2533 
2534  ft(il,i)=ft(il,i)+0.01*sigd*wt(il,i)*(cl-cpd)*water(il,i+1)
2535  : *(t(il,i+1)-t(il,i))*dpinv*cpinv
2536 
2537  if (cvflag_grav) then
2538  fr(il,i)=0.01*grav*dpinv*(amp1(il)*(rr(il,i+1)-rr(il,i))
2539  : -ad(il)*(rr(il,i)-rr(il,i-1)))
2540  fu(il,i)=fu(il,i)+0.01*grav*dpinv*(amp1(il)*(u(il,i+1)-u(il,i))
2541  : -ad(il)*(u(il,i)-u(il,i-1)))
2542  fv(il,i)=fv(il,i)+0.01*grav*dpinv*(amp1(il)*(v(il,i+1)-v(il,i))
2543  : -ad(il)*(v(il,i)-v(il,i-1)))
2544  else ! cvflag_grav
2545  fr(il,i)=0.1*dpinv*(amp1(il)*(rr(il,i+1)-rr(il,i))
2546  : -ad(il)*(rr(il,i)-rr(il,i-1)))
2547  fu(il,i)=fu(il,i)+0.1*dpinv*(amp1(il)*(u(il,i+1)-u(il,i))
2548  : -ad(il)*(u(il,i)-u(il,i-1)))
2549  fv(il,i)=fv(il,i)+0.1*dpinv*(amp1(il)*(v(il,i+1)-v(il,i))
2550  : -ad(il)*(v(il,i)-v(il,i-1)))
2551  endif ! cvflag_grav
2552 
2553  endif ! i
2554 1350 continue
2555 
2556 c do k=1,ntra
2557 c do il=1,ncum
2558 c if (i.le.inb(il)) then
2559 c dpinv=1.0/(ph(il,i)-ph(il,i+1))
2560 c cpinv=1.0/cpn(il,i)
2561 c if (cvflag_grav) then
2562 c ftra(il,i,k)=ftra(il,i,k)+0.01*grav*dpinv
2563 c : *(amp1(il)*(tra(il,i+1,k)-tra(il,i,k))
2564 c : -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
2565 c else
2566 c ftra(il,i,k)=ftra(il,i,k)+0.1*dpinv
2567 c : *(amp1(il)*(tra(il,i+1,k)-tra(il,i,k))
2568 c : -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
2569 c endif
2570 c endif
2571 c enddo
2572 c enddo
2573 
2574  do 480 k=1,i-1
2575  do 1370 il=1,ncum
2576  if (i.le.inb(il)) then
2577  dpinv=1.0/(ph(il,i)-ph(il,i+1))
2578  cpinv=1.0/cpn(il,i)
2579 
2580  awat=elij(il,k,i)-(1.-ep(il,i))*clw(il,i)
2581  awat=amax1(awat,0.0)
2582 
2583  if (cvflag_grav) then
2584  fr(il,i)=fr(il,i)
2585  : +0.01*grav*dpinv*ment(il,k,i)*(qent(il,k,i)-awat-rr(il,i))
2586  fu(il,i)=fu(il,i)
2587  : +0.01*grav*dpinv*ment(il,k,i)*(uent(il,k,i)-u(il,i))
2588  fv(il,i)=fv(il,i)
2589  : +0.01*grav*dpinv*ment(il,k,i)*(vent(il,k,i)-v(il,i))
2590  else ! cvflag_grav
2591  fr(il,i)=fr(il,i)
2592  : +0.1*dpinv*ment(il,k,i)*(qent(il,k,i)-awat-rr(il,i))
2593  fu(il,i)=fu(il,i)
2594  : +0.01*grav*dpinv*ment(il,k,i)*(uent(il,k,i)-u(il,i))
2595  fv(il,i)=fv(il,i)
2596  : +0.1*dpinv*ment(il,k,i)*(vent(il,k,i)-v(il,i))
2597  endif ! cvflag_grav
2598 
2599 c (saturated updrafts resulting from mixing) ! cld
2600  qcond(il,i)=qcond(il,i)+(elij(il,k,i)-awat) ! cld
2601  nqcond(il,i)=nqcond(il,i)+1. ! cld
2602  endif ! i
2603 1370 continue
2604 480 continue
2605 
2606 c do j=1,ntra
2607 c do k=1,i-1
2608 c do il=1,ncum
2609 c if (i.le.inb(il)) then
2610 c dpinv=1.0/(ph(il,i)-ph(il,i+1))
2611 c cpinv=1.0/cpn(il,i)
2612 c if (cvflag_grav) then
2613 c ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv*ment(il,k,i)
2614 c : *(traent(il,k,i,j)-tra(il,i,j))
2615 c else
2616 c ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv*ment(il,k,i)
2617 c : *(traent(il,k,i,j)-tra(il,i,j))
2618 c endif
2619 c endif
2620 c enddo
2621 c enddo
2622 c enddo
2623 
2624  do 490 k=i,nl+1
2625  do 1380 il=1,ncum
2626  if (i.le.inb(il) .and. k.le.inb(il)) then
2627  dpinv=1.0/(ph(il,i)-ph(il,i+1))
2628  cpinv=1.0/cpn(il,i)
2629 
2630  if (cvflag_grav) then
2631  fr(il,i)=fr(il,i)
2632  : +0.01*grav*dpinv*ment(il,k,i)*(qent(il,k,i)-rr(il,i))
2633  fu(il,i)=fu(il,i)
2634  : +0.01*grav*dpinv*ment(il,k,i)*(uent(il,k,i)-u(il,i))
2635  fv(il,i)=fv(il,i)
2636  : +0.01*grav*dpinv*ment(il,k,i)*(vent(il,k,i)-v(il,i))
2637  else ! cvflag_grav
2638  fr(il,i)=fr(il,i)
2639  : +0.1*dpinv*ment(il,k,i)*(qent(il,k,i)-rr(il,i))
2640  fu(il,i)=fu(il,i)
2641  : +0.1*dpinv*ment(il,k,i)*(uent(il,k,i)-u(il,i))
2642  fv(il,i)=fv(il,i)
2643  : +0.1*dpinv*ment(il,k,i)*(vent(il,k,i)-v(il,i))
2644  endif ! cvflag_grav
2645  endif ! i and k
2646 1380 continue
2647 490 continue
2648 
2649 c do j=1,ntra
2650 c do k=i,nl+1
2651 c do il=1,ncum
2652 c if (i.le.inb(il) .and. k.le.inb(il)) then
2653 c dpinv=1.0/(ph(il,i)-ph(il,i+1))
2654 c cpinv=1.0/cpn(il,i)
2655 c if (cvflag_grav) then
2656 c ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv*ment(il,k,i)
2657 c : *(traent(il,k,i,j)-tra(il,i,j))
2658 c else
2659 c ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv*ment(il,k,i)
2660 c : *(traent(il,k,i,j)-tra(il,i,j))
2661 c endif
2662 c endif ! i and k
2663 c enddo
2664 c enddo
2665 c enddo
2666 
2667  do 1400 il=1,ncum
2668  if (i.le.inb(il)) then
2669  dpinv=1.0/(ph(il,i)-ph(il,i+1))
2670  cpinv=1.0/cpn(il,i)
2671 
2672  if (cvflag_grav) then
2673 c sb: on ne fait pas encore la correction permettant de mieux
2674 c conserver l'eau:
2675  fr(il,i)=fr(il,i)+0.5*sigd*(evap(il,i)+evap(il,i+1))
2676  : +0.01*grav*(mp(il,i+1)*(rp(il,i+1)-rr(il,i))-mp(il,i)
2677  : *(rp(il,i)-rr(il,i-1)))*dpinv
2678 
2679  fu(il,i)=fu(il,i)+0.01*grav*(mp(il,i+1)*(up(il,i+1)-u(il,i))
2680  : -mp(il,i)*(up(il,i)-u(il,i-1)))*dpinv
2681  fv(il,i)=fv(il,i)+0.01*grav*(mp(il,i+1)*(vp(il,i+1)-v(il,i))
2682  : -mp(il,i)*(vp(il,i)-v(il,i-1)))*dpinv
2683  else ! cvflag_grav
2684  fr(il,i)=fr(il,i)+0.5*sigd*(evap(il,i)+evap(il,i+1))
2685  : +0.1*(mp(il,i+1)*(rp(il,i+1)-rr(il,i))-mp(il,i)
2686  : *(rp(il,i)-rr(il,i-1)))*dpinv
2687  fu(il,i)=fu(il,i)+0.1*(mp(il,i+1)*(up(il,i+1)-u(il,i))
2688  : -mp(il,i)*(up(il,i)-u(il,i-1)))*dpinv
2689  fv(il,i)=fv(il,i)+0.1*(mp(il,i+1)*(vp(il,i+1)-v(il,i))
2690  : -mp(il,i)*(vp(il,i)-v(il,i-1)))*dpinv
2691  endif ! cvflag_grav
2692 
2693  endif ! i
2694 1400 continue
2695 
2696 c sb: interface with the cloud parameterization: ! cld
2697 
2698  do k=i+1,nl
2699  do il=1,ncum
2700  if (k.le.inb(il) .and. i.le.inb(il)) then ! cld
2701 C (saturated downdrafts resulting from mixing) ! cld
2702  qcond(il,i)=qcond(il,i)+elij(il,k,i) ! cld
2703  nqcond(il,i)=nqcond(il,i)+1. ! cld
2704  endif ! cld
2705  enddo ! cld
2706  enddo ! cld
2707 
2708 C (particular case: no detraining level is found) ! cld
2709  do il=1,ncum ! cld
2710  if (i.le.inb(il) .and. nent(il,i).eq.0) then ! cld
2711  qcond(il,i)=qcond(il,i)+(1.-ep(il,i))*clw(il,i) ! cld
2712  nqcond(il,i)=nqcond(il,i)+1. ! cld
2713  endif ! cld
2714  enddo ! cld
2715 
2716  do il=1,ncum ! cld
2717  if (i.le.inb(il) .and. nqcond(il,i).ne.0.) then ! cld
2718  qcond(il,i)=qcond(il,i)/nqcond(il,i) ! cld
2719  endif ! cld
2720  enddo
2721 
2722 c do j=1,ntra
2723 c do il=1,ncum
2724 c if (i.le.inb(il)) then
2725 c dpinv=1.0/(ph(il,i)-ph(il,i+1))
2726 c cpinv=1.0/cpn(il,i)
2727 
2728 c if (cvflag_grav) then
2729 c ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv
2730 c : *(mp(il,i+1)*(trap(il,i+1,j)-tra(il,i,j))
2731 c : -mp(il,i)*(trap(il,i,j)-tra(il,i-1,j)))
2732 c else
2733 c ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv
2734 c : *(mp(il,i+1)*(trap(il,i+1,j)-tra(il,i,j))
2735 c : -mp(il,i)*(trap(il,i,j)-tra(il,i-1,j)))
2736 c endif
2737 c endif ! i
2738 c enddo
2739 c enddo
2740 
2741 500 continue
2742 
2743 
2744 c *** move the detrainment at level inb down to level inb-1 ***
2745 c *** in such a way as to preserve the vertically ***
2746 c *** integrated enthalpy and water tendencies ***
2747 c
2748  do 503 il=1,ncum
2749 
2750  ax=0.1*ment(il,inb(il),inb(il))*(hp(il,inb(il))-h(il,inb(il))
2751  : +t(il,inb(il))*(cpv-cpd)
2752  : *(rr(il,inb(il))-qent(il,inb(il),inb(il))))
2753  : /(cpn(il,inb(il))*(ph(il,inb(il))-ph(il,inb(il)+1)))
2754  ft(il,inb(il))=ft(il,inb(il))-ax
2755  ft(il,inb(il)-1)=ft(il,inb(il)-1)+ax*cpn(il,inb(il))
2756  : *(ph(il,inb(il))-ph(il,inb(il)+1))/(cpn(il,inb(il)-1)
2757  : *(ph(il,inb(il)-1)-ph(il,inb(il))))
2758 
2759  bx=0.1*ment(il,inb(il),inb(il))*(qent(il,inb(il),inb(il))
2760  : -rr(il,inb(il)))/(ph(il,inb(il))-ph(il,inb(il)+1))
2761  fr(il,inb(il))=fr(il,inb(il))-bx
2762  fr(il,inb(il)-1)=fr(il,inb(il)-1)
2763  : +bx*(ph(il,inb(il))-ph(il,inb(il)+1))
2764  : /(ph(il,inb(il)-1)-ph(il,inb(il)))
2765 
2766  cx=0.1*ment(il,inb(il),inb(il))*(uent(il,inb(il),inb(il))
2767  : -u(il,inb(il)))/(ph(il,inb(il))-ph(il,inb(il)+1))
2768  fu(il,inb(il))=fu(il,inb(il))-cx
2769  fu(il,inb(il)-1)=fu(il,inb(il)-1)
2770  : +cx*(ph(il,inb(il))-ph(il,inb(il)+1))
2771  : /(ph(il,inb(il)-1)-ph(il,inb(il)))
2772 
2773  dx=0.1*ment(il,inb(il),inb(il))*(vent(il,inb(il),inb(il))
2774  : -v(il,inb(il)))/(ph(il,inb(il))-ph(il,inb(il)+1))
2775  fv(il,inb(il))=fv(il,inb(il))-dx
2776  fv(il,inb(il)-1)=fv(il,inb(il)-1)
2777  : +dx*(ph(il,inb(il))-ph(il,inb(il)+1))
2778  : /(ph(il,inb(il)-1)-ph(il,inb(il)))
2779 
2780 503 continue
2781 
2782 c do j=1,ntra
2783 c do il=1,ncum
2784 c ex=0.1*ment(il,inb(il),inb(il))
2785 c : *(traent(il,inb(il),inb(il),j)-tra(il,inb(il),j))
2786 c : /(ph(il,inb(il))-ph(il,inb(il)+1))
2787 c ftra(il,inb(il),j)=ftra(il,inb(il),j)-ex
2788 c ftra(il,inb(il)-1,j)=ftra(il,inb(il)-1,j)
2789 c : +ex*(ph(il,inb(il))-ph(il,inb(il)+1))
2790 c : /(ph(il,inb(il)-1)-ph(il,inb(il)))
2791 c enddo
2792 c enddo
2793 
2794 c
2795 c *** homoginize tendencies below cloud base ***
2796 c
2797 c
2798  do il=1,ncum
2799  asum(il)=0.0
2800  bsum(il)=0.0
2801  csum(il)=0.0
2802  dsum(il)=0.0
2803  enddo
2804 
2805  do i=1,nl
2806  do il=1,ncum
2807  if (i.le.(icb(il)-1)) then
2808  asum(il)=asum(il)+ft(il,i)*(ph(il,i)-ph(il,i+1))
2809  bsum(il)=bsum(il)+fr(il,i)*(lv(il,i)+(cl-cpd)*(t(il,i)-t(il,1)))
2810  : *(ph(il,i)-ph(il,i+1))
2811  csum(il)=csum(il)+(lv(il,i)+(cl-cpd)*(t(il,i)-t(il,1)))
2812  : *(ph(il,i)-ph(il,i+1))
2813  dsum(il)=dsum(il)+t(il,i)*(ph(il,i)-ph(il,i+1))/th(il,i)
2814  endif
2815  enddo
2816  enddo
2817 
2818 c!!! do 700 i=1,icb(il)-1
2819  do i=1,nl
2820  do il=1,ncum
2821  if (i.le.(icb(il)-1)) then
2822  ft(il,i)=asum(il)*t(il,i)/(th(il,i)*dsum(il))
2823  fr(il,i)=bsum(il)/csum(il)
2824  endif
2825  enddo
2826  enddo
2827 
2828 c
2829 c *** reset counter and return ***
2830 c
2831  do il=1,ncum
2832  sig(il,nd)=2.0
2833  enddo
2834 
2835 
2836  do i=1,nd
2837  do il=1,ncum
2838  upwd(il,i)=0.0
2839  dnwd(il,i)=0.0
2840  enddo
2841  enddo
2842 
2843  do i=1,nl
2844  do il=1,ncum
2845  dnwd0(il,i)=-mp(il,i)
2846  enddo
2847  enddo
2848  do i=nl+1,nd
2849  do il=1,ncum
2850  dnwd0(il,i)=0.
2851  enddo
2852  enddo
2853 
2854 
2855  do i=1,nl
2856  do il=1,ncum
2857  if (i.ge.icb(il) .and. i.le.inb(il)) then
2858  upwd(il,i)=0.0
2859  dnwd(il,i)=0.0
2860  endif
2861  enddo
2862  enddo
2863 
2864  do i=1,nl
2865  do k=1,nl
2866  do il=1,ncum
2867  up1(il,k,i)=0.0
2868  dn1(il,k,i)=0.0
2869  enddo
2870  enddo
2871  enddo
2872 
2873  do i=1,nl
2874  do k=i,nl
2875  do n=1,i-1
2876  do il=1,ncum
2877  if (i.ge.icb(il).and.i.le.inb(il).and.k.le.inb(il)) then
2878  up1(il,k,i)=up1(il,k,i)+ment(il,n,k)
2879  dn1(il,k,i)=dn1(il,k,i)-ment(il,k,n)
2880  endif
2881  enddo
2882  enddo
2883  enddo
2884  enddo
2885 
2886  do i=2,nl
2887  do k=i,nl
2888  do il=1,ncum
2889 ctest if (i.ge.icb(il).and.i.le.inb(il).and.k.le.inb(il)) then
2890  if (i.le.inb(il).and.k.le.inb(il)) then
2891  upwd(il,i)=upwd(il,i)+m(il,k)+up1(il,k,i)
2892  dnwd(il,i)=dnwd(il,i)+dn1(il,k,i)
2893  endif
2894  enddo
2895  enddo
2896  enddo
2897 
2898 
2899 c!!! DO il=1,ncum
2900 c!!! do i=icb(il),inb(il)
2901 c!!!
2902 c!!! upwd(il,i)=0.0
2903 c!!! dnwd(il,i)=0.0
2904 c!!! do k=i,inb(il)
2905 c!!! up1=0.0
2906 c!!! dn1=0.0
2907 c!!! do n=1,i-1
2908 c!!! up1=up1+ment(il,n,k)
2909 c!!! dn1=dn1-ment(il,k,n)
2910 c!!! enddo
2911 c!!! upwd(il,i)=upwd(il,i)+m(il,k)+up1
2912 c!!! dnwd(il,i)=dnwd(il,i)+dn1
2913 c!!! enddo
2914 c!!! enddo
2915 c!!!
2916 c!!! ENDDO
2917 
2918 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2919 c determination de la variation de flux ascendant entre
2920 c deux niveau non dilue mike
2921 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2922 
2923  do i=1,nl
2924  do il=1,ncum
2925  mike(il,i)=m(il,i)
2926  enddo
2927  enddo
2928 
2929  do i=nl+1,nd
2930  do il=1,ncum
2931  mike(il,i)=0.
2932  enddo
2933  enddo
2934 
2935  do i=1,nd
2936  do il=1,ncum
2937  ma(il,i)=0
2938  enddo
2939  enddo
2940 
2941  do i=1,nl
2942  do j=i,nl
2943  do il=1,ncum
2944  ma(il,i)=ma(il,i)+m(il,j)
2945  enddo
2946  enddo
2947  enddo
2948 
2949  do i=nl+1,nd
2950  do il=1,ncum
2951  ma(il,i)=0.
2952  enddo
2953  enddo
2954 
2955  do i=1,nl
2956  do il=1,ncum
2957  if (i.le.(icb(il)-1)) then
2958  ma(il,i)=0
2959  endif
2960  enddo
2961  enddo
2962 
2963 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2964 c icb represente de niveau ou se trouve la
2965 c base du nuage , et inb le top du nuage
2966 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2967 
2968  do i=1,nd
2969  do il=1,ncum
2970  mke(il,i)=upwd(il,i)+dnwd(il,i)
2971  enddo
2972  enddo
2973 
2974  do i=1,nd
2975  DO 999 il=1,ncum
2976  rdcp=(rrd*(1.-rr(il,i))-rr(il,i)*rrv)
2977  : /(cpd*(1.-rr(il,i))+rr(il,i)*cpv)
2978  tls(il,i)=t(il,i)*(1000.0/p(il,i))**rdcp
2979  tps(il,i)=tp(il,i)
2980 999 CONTINUE
2981  enddo
2982 
2983 c
2984 c *** diagnose the in-cloud mixing ratio *** ! cld
2985 c *** of condensed water *** ! cld
2986 c ! cld
2987 
2988  do i=1,nd ! cld
2989  do il=1,ncum ! cld
2990  mac(il,i)=0.0 ! cld
2991  wa(il,i)=0.0 ! cld
2992  siga(il,i)=0.0 ! cld
2993  sax(il,i)=0.0 ! cld
2994  enddo ! cld
2995  enddo ! cld
2996 
2997  do i=minorig, nl ! cld
2998  do k=i+1,nl+1 ! cld
2999  do il=1,ncum ! cld
3000  if (i.le.inb(il) .and. k.le.(inb(il)+1)) then ! cld
3001  mac(il,i)=mac(il,i)+m(il,k) ! cld
3002  endif ! cld
3003  enddo ! cld
3004  enddo ! cld
3005  enddo ! cld
3006 
3007  do i=1,nl ! cld
3008  do j=1,i ! cld
3009  do il=1,ncum ! cld
3010  if (i.ge.icb(il) .and. i.le.(inb(il)-1) ! cld
3011  : .and. j.ge.icb(il) ) then ! cld
3012  sax(il,i)=sax(il,i)+rrd*(tvp(il,j)-tv(il,j)) ! cld
3013  : *(ph(il,j)-ph(il,j+1))/p(il,j) ! cld
3014  endif ! cld
3015  enddo ! cld
3016  enddo ! cld
3017  enddo ! cld
3018 
3019  do i=1,nl ! cld
3020  do il=1,ncum ! cld
3021  if (i.ge.icb(il) .and. i.le.(inb(il)-1) ! cld
3022  : .and. sax(il,i).gt.0.0 ) then ! cld
3023  wa(il,i)=sqrt(2.*sax(il,i)) ! cld
3024  endif ! cld
3025  enddo ! cld
3026  enddo ! cld
3027 
3028  do i=1,nl ! cld
3029  do il=1,ncum ! cld
3030  if (wa(il,i).gt.0.0) ! cld
3031  : siga(il,i)=mac(il,i)/wa(il,i) ! cld
3032  : *rrd*tvp(il,i)/p(il,i)/100./delta ! cld
3033  siga(il,i) = min(siga(il,i),1.0) ! cld
3034 cIM cf. FH
3035  if (iflag_clw.eq.0) then
3036  qcondc(il,i)=siga(il,i)*clw(il,i)*(1.-ep(il,i)) ! cld
3037  : + (1.-siga(il,i))*qcond(il,i) ! cld
3038  else if (iflag_clw.eq.1) then
3039  qcondc(il,i)=qcond(il,i) ! cld
3040  endif
3041 
3042  enddo ! cld
3043  enddo ! cld
3044 
3045  return
3046  end
3047 
3048 !!RomP >>>
3049  SUBROUTINE cv30_tracer(nloc,len,ncum,nd,na,
3050  & ment,sij,da,phi,phi2,d1a,dam,
3051  & ep,vprecip,elij,clw,epmlmmm,eplamm,
3052  & icb,inb)
3053  implicit none
3054 
3055 #include "cv30param.h"
3056 
3057 c inputs:
3058  integer ncum, nd, na, nloc,len
3059  real ment(nloc,na,na),sij(nloc,na,na)
3060  real clw(nloc,nd),elij(nloc,na,na)
3061  real ep(nloc,na)
3062  integer icb(nloc),inb(nloc)
3063  real vprecip(nloc,nd+1)
3064 c ouputs:
3065  real da(nloc,na),phi(nloc,na,na)
3066  real phi2(nloc,na,na)
3067  real d1a(nloc,na),dam(nloc,na)
3068  real epmlmmm(nloc,na,na),eplamm(nloc,na)
3069 ! variables pour tracer dans precip de l'AA et des mel
3070 c local variables:
3071  integer i,j,k
3072  real epm(nloc,na,na)
3073 c
3074 ! variables d'Emanuel : du second indice au troisieme
3075 ! ---> tab(i,k,j) -> de l origine k a l arrivee j
3076 ! ment, sij, elij
3077 ! variables personnelles : du troisieme au second indice
3078 ! ---> tab(i,j,k) -> de k a j
3079 ! phi, phi2
3080 !
3081 ! initialisations
3082  do j=1,na
3083  do i=1,ncum
3084  da(i,j)=0.
3085  d1a(i,j)=0.
3086  dam(i,j)=0.
3087  eplamm(i,j)=0.
3088  enddo
3089  enddo
3090  do k=1,na
3091  do j=1,na
3092  do i=1,ncum
3093  epm(i,j,k)=0.
3094  epmlmmm(i,j,k)=0.
3095  phi(i,j,k)=0.
3096  phi2(i,j,k)=0.
3097  enddo
3098  enddo
3099  enddo
3100 c
3101 ! fraction deau condensee dans les melanges convertie en precip : epm
3102 ! et eau condensée précipitée dans masse d'air saturé : l_m*dM_m/dzdz.dzdz
3103  do j=1,na
3104  do k=1,j-1
3105  do i=1,ncum
3106  if(k.ge.icb(i).and.k.le.inb(i).and.
3107  & j.le.inb(i)) then
3108 !!jyg epm(i,j,k)=1.-(1.-ep(i,j))*clw(i,j)/elij(i,k,j)
3109  epm(i,j,k)=1.-(1.-ep(i,j))*clw(i,j)/
3110  & max(elij(i,k,j),1.e-16)
3111 !!
3112  epm(i,j,k)=max(epm(i,j,k),0.0)
3113  endif
3114  end do
3115  end do
3116  end do
3117 !
3118  do j=1,na
3119  do k=1,na
3120  do i=1,ncum
3121  if(k.ge.icb(i).and.k.le.inb(i)) then
3122  eplamm(i,j)=eplamm(i,j) + ep(i,j)*clw(i,j)
3123  & *ment(i,j,k)*(1.-sij(i,j,k))
3124  endif
3125  end do
3126  end do
3127  end do
3128 !
3129  do j=1,na
3130  do k=1,j-1
3131  do i=1,ncum
3132  if(k.ge.icb(i).and.k.le.inb(i).and.
3133  & j.le.inb(i)) then
3134  epmlmmm(i,j,k)=epm(i,j,k)*elij(i,k,j)*ment(i,k,j)
3135  endif
3136  end do
3137  end do
3138  end do
3139 
3140 ! matrices pour calculer la tendance des concentrations dans cvltr.F90
3141  do j=1,na
3142  do k=1,na
3143  do i=1,ncum
3144  da(i,j)=da(i,j)+(1.-sij(i,k,j))*ment(i,k,j)
3145  phi(i,j,k)=sij(i,k,j)*ment(i,k,j)
3146  d1a(i,j)=d1a(i,j)+ment(i,k,j)*ep(i,k)
3147  & *(1.-sij(i,k,j))
3148  end do
3149  end do
3150  end do
3151 
3152  do j=1,na
3153  do k=1,j-1
3154  do i=1,ncum
3155  dam(i,j)=dam(i,j)+ment(i,k,j)
3156  & *epm(i,j,k)*(1.-ep(i,k))*(1.-sij(i,k,j))
3157  phi2(i,j,k)=phi(i,j,k)*epm(i,j,k)
3158  end do
3159  end do
3160  end do
3161 
3162  return
3163  end
3164 !RomP <<<
3165 
3166  SUBROUTINE cv30_uncompress(nloc,len,ncum,nd,ntra,idcum
3167  : ,iflag
3168  : ,precip,vprecip,evap,ep,sig,w0
3169  : ,ft,fq,fu,fv,ftra
3170  : ,inb
3171  : ,ma,upwd,dnwd,dnwd0,qcondc,wd,cape
3172  : ,da,phi,mp,phi2,d1a,dam,sij
3173  : ,elij,clw,epmlmmm,eplamm
3174  : ,wdtraina,wdtrainm
3175  : ,iflag1
3176  : ,precip1,vprecip1,evap1,ep1,sig1,w01
3177  : ,ft1,fq1,fu1,fv1,ftra1
3178  : ,inb1
3179  : ,ma1,upwd1,dnwd1,dnwd01,qcondc1,wd1,cape1
3180  : ,da1,phi1,mp1,phi21,d1a1,dam1,sij1
3181  : ,elij1,clw1,epmlmmm1,eplamm1
3182  : ,wdtraina1,wdtrainm1)
3183  implicit none
3184 
3185 #include "cv30param.h"
3186 
3187 c inputs:
3188  integer len, ncum, nd, ntra, nloc
3189  integer idcum(nloc)
3190  integer iflag(nloc)
3191  integer inb(nloc)
3192  real precip(nloc)
3193  real vprecip(nloc,nd+1),evap(nloc,nd)
3194  real ep(nloc,nd)
3195  real sig(nloc,nd), w0(nloc,nd)
3196  real ft(nloc,nd), fq(nloc,nd), fu(nloc,nd), fv(nloc,nd)
3197  real ftra(nloc,nd,ntra)
3198  real ma(nloc,nd)
3199  real upwd(nloc,nd),dnwd(nloc,nd),dnwd0(nloc,nd)
3200  real qcondc(nloc,nd)
3201  real wd(nloc),cape(nloc)
3202  real da(nloc,nd),phi(nloc,nd,nd),mp(nloc,nd)
3203 !RomP >>>
3204  real phi2(nloc,nd,nd)
3205  real d1a(nloc,nd),dam(nloc,nd)
3206  real wdtraina(nloc,nd), wdtrainm(nloc,nd)
3207  real sij(nloc,nd,nd)
3208  real elij(nloc,nd,nd),clw(nloc,nd)
3209  real epmlmmm(nloc,nd,nd),eplamm(nloc,nd)
3210 !RomP <<<
3211 
3212 c outputs:
3213  integer iflag1(len)
3214  integer inb1(len)
3215  real precip1(len)
3216  real vprecip1(len,nd+1),evap1(len,nd)
3217  real ep1(len,nd)
3218  real sig1(len,nd), w01(len,nd)
3219  real ft1(len,nd), fq1(len,nd), fu1(len,nd), fv1(len,nd)
3220  real ftra1(len,nd,ntra)
3221  real ma1(len,nd)
3222  real upwd1(len,nd),dnwd1(len,nd),dnwd01(len,nd)
3223  real qcondc1(nloc,nd)
3224  real wd1(nloc),cape1(nloc)
3225  real da1(nloc,nd),phi1(nloc,nd,nd),mp1(nloc,nd)
3226 !RomP >>>
3227  real phi21(len,nd,nd)
3228  real d1a1(len,nd),dam1(len,nd)
3229  real wdtraina1(len,nd), wdtrainm1(len,nd)
3230  real sij1(len,nd,nd)
3231  real elij1(len,nd,nd),clw1(len,nd)
3232  real epmlmmm1(len,nd,nd),eplamm1(len,nd)
3233 !RomP <<<
3234 
3235 c local variables:
3236  integer i,k,j
3237 
3238  do 2000 i=1,ncum
3239  precip1(idcum(i))=precip(i)
3240  iflag1(idcum(i))=iflag(i)
3241  wd1(idcum(i))=wd(i)
3242  inb1(idcum(i))=inb(i)
3243  cape1(idcum(i))=cape(i)
3244  2000 continue
3245 
3246  do 2020 k=1,nl
3247  do 2010 i=1,ncum
3248  vprecip1(idcum(i),k)=vprecip(i,k)
3249  evap1(idcum(i),k)=evap(i,k)
3250  sig1(idcum(i),k)=sig(i,k)
3251  w01(idcum(i),k)=w0(i,k)
3252  ft1(idcum(i),k)=ft(i,k)
3253  fq1(idcum(i),k)=fq(i,k)
3254  fu1(idcum(i),k)=fu(i,k)
3255  fv1(idcum(i),k)=fv(i,k)
3256  ma1(idcum(i),k)=ma(i,k)
3257  upwd1(idcum(i),k)=upwd(i,k)
3258  dnwd1(idcum(i),k)=dnwd(i,k)
3259  dnwd01(idcum(i),k)=dnwd0(i,k)
3260  qcondc1(idcum(i),k)=qcondc(i,k)
3261  da1(idcum(i),k)=da(i,k)
3262  mp1(idcum(i),k)=mp(i,k)
3263 !RomP >>>
3264  ep1(idcum(i),k)=ep(i,k)
3265  d1a1(idcum(i),k)=d1a(i,k)
3266  dam1(idcum(i),k)=dam(i,k)
3267  clw1(idcum(i),k)=clw(i,k)
3268  eplamm1(idcum(i),k)=eplamm(i,k)
3269  wdtraina1(idcum(i),k)=wdtraina(i,k)
3270  wdtrainm1(idcum(i),k)=wdtrainm(i,k)
3271 !RomP <<<
3272  2010 continue
3273  2020 continue
3274 
3275  do 2200 i=1,ncum
3276  sig1(idcum(i),nd)=sig(i,nd)
3277 2200 continue
3278 
3279 
3280 c do 2100 j=1,ntra
3281 c do 2110 k=1,nd ! oct3
3282 c do 2120 i=1,ncum
3283 c ftra1(idcum(i),k,j)=ftra(i,k,j)
3284 c 2120 continue
3285 c 2110 continue
3286 c 2100 continue
3287  do j=1,nd
3288  do k=1,nd
3289  do i=1,ncum
3290  sij1(idcum(i),k,j)=sij(i,k,j)
3291  phi1(idcum(i),k,j)=phi(i,k,j)
3292  phi21(idcum(i),k,j)=phi2(i,k,j)
3293  elij1(idcum(i),k,j)=elij(i,k,j)
3294  epmlmmm1(idcum(i),k,j)=epmlmmm(i,k,j)
3295  end do
3296  end do
3297  end do
3298 
3299  return
3300  end
3301