source: LMDZ6/branches/Ocean_skin/libf/phylmd/cv3_routines.F90 @ 3605

Last change on this file since 3605 was 3605, checked in by lguez, 5 years ago

Merge revisions 3427:3600 of trunk into branch Ocean_skin

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