source: LMDZ6/trunk/libf/phylmd/cv3_routines.f90 @ 5503

Last change on this file since 5503 was 5502, checked in by fhourdin, 11 days ago

Flag moved to lmdz_cv_ini/cv3_param.

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