source: LMDZ5/branches/testing/libf/phylmd/cv3_routines.F90 @ 1999

Last change on this file since 1999 was 1999, checked in by Laurent Fairhead, 10 years ago

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