source: LMDZ5/trunk/libf/phylmd/cv3_routines.F90 @ 2306

Last change on this file since 2306 was 2306, checked in by jyg, 9 years ago

Improved diagnostics: pmflxr and pmflxs are now
the true vertical profiles of liquid and solid
convective precipitation (previously they where
merely diagnosed through a test on temperature).
The convective scheme internal variable for ice
precipitation is Vprecipi.

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