source: LMDZ6/trunk/libf/phylmd/cv3_routines.F90 @ 3741

Last change on this file since 3741 was 3670, checked in by jyg, 5 years ago

In order to makEmanuel scheme : set iflag=14 in
cv3_trigger when
no convection and reset Cin to Cin_noconv at the
end of cva_driver when iflag=14.

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