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