source: LMDZ5/branches/testing/libf/phylmd/cv3_routines.F90 @ 5431

Last change on this file since 5431 was 2839, checked in by Laurent Fairhead, 8 years ago

Merged trunk changes r2785:2838 into testing branch

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