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

Last change on this file since 5296 was 5285, checked in by abarral, 9 months ago

As discussed internally, remove generic ONLY: ... for new _mod_h modules

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