source: LMDZ6/branches/DYNAMICO-conv-GC/libf/phylmd/cv3_routines.F90 @ 5214

Last change on this file since 5214 was 3406, checked in by jghattas, 6 years ago

Added all modifications in the model code that were used for the simulations with DYANMICO during the Grand Challeng 2018. Modifications done by Y. Meurdesoif, L. Fairhead and A.K. Traore

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