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

Last change on this file since 5747 was 5707, checked in by yann meurdesoif, 3 weeks ago

Convection GPU porting : suppress potential dependency between columns, may change results (cv3_closure)

YM

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