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

Last change on this file since 5445 was 5346, checked in by fhourdin, 3 weeks ago

Debut de replaysation de la convection profonde.

Regroupement de cvparam, cv3param et cvthermo (récemment
passés de statut de .h à module, dans un unique module
lmdz_cv_ini.f90

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