Ignore:
Timestamp:
Apr 6, 2014, 2:20:38 PM (10 years ago)
Author:
fhourdin
Message:

Modification pour garantir la conservation des espèces traces.

Modification for tracer conservation

Jean-Yves Grandpeix

File:
1 edited

Legend:

Unmodified
Added
Removed

  • LMDZ5/trunk/libf/phylmd/cv3_routines.F90

    r1992 r2007  
    77  IMPLICIT NONE
    88
    9   ! ------------------------------------------------------------
    10   ! Set parameters for convectL for iflag_con = 3
    11   ! ------------------------------------------------------------
    12 
    13 
    14   ! ***  PBCRIT IS THE CRITICAL CLOUD DEPTH (MB) BENEATH WHICH THE ***
    15   ! ***      PRECIPITATION EFFICIENCY IS ASSUMED TO BE ZERO     ***
    16   ! ***  PTCRIT IS THE CLOUD DEPTH (MB) ABOVE WHICH THE PRECIP. ***
    17   ! ***            EFFICIENCY IS ASSUMED TO BE UNITY            ***
    18   ! ***  SIGD IS THE FRACTIONAL AREA COVERED BY UNSATURATED DNDRAFT  ***
    19   ! ***  SPFAC IS THE FRACTION OF PRECIPITATION FALLING OUTSIDE ***
    20   ! ***                        OF CLOUD                         ***
    21 
    22   ! [TAU: CHARACTERISTIC TIMESCALE USED TO COMPUTE ALPHA & BETA]
    23   ! ***    ALPHA AND BETA ARE PARAMETERS THAT CONTROL THE RATE OF ***
    24   ! ***                 APPROACH TO QUASI-EQUILIBRIUM           ***
    25   ! ***    (THEIR STANDARD VALUES ARE 1.0 AND 0.96, RESPECTIVELY) ***
    26   ! ***           (BETA MUST BE LESS THAN OR EQUAL TO 1)        ***
    27 
    28   ! ***    DTCRIT IS THE CRITICAL BUOYANCY (K) USED TO ADJUST THE ***
    29   ! ***                 APPROACH TO QUASI-EQUILIBRIUM           ***
    30   ! ***                     IT MUST BE LESS THAN 0              ***
     9!------------------------------------------------------------
     10!Set parameters for convectL for iflag_con = 3
     11!------------------------------------------------------------
     12
     13
     14!***  PBCRIT IS THE CRITICAL CLOUD DEPTH (MB) BENEATH WHICH THE ***
     15!***      PRECIPITATION EFFICIENCY IS ASSUMED TO BE ZERO     ***
     16!***  PTCRIT IS THE CLOUD DEPTH (MB) ABOVE WHICH THE PRECIP. ***
     17!***            EFFICIENCY IS ASSUMED TO BE UNITY            ***
     18!***  SIGD IS THE FRACTIONAL AREA COVERED BY UNSATURATED DNDRAFT  ***
     19!***  SPFAC IS THE FRACTION OF PRECIPITATION FALLING OUTSIDE ***
     20!***                        OF CLOUD                         ***
     21
     22![TAU: CHARACTERISTIC TIMESCALE USED TO COMPUTE ALPHA & BETA]
     23!***    ALPHA AND BETA ARE PARAMETERS THAT CONTROL THE RATE OF ***
     24!***                 APPROACH TO QUASI-EQUILIBRIUM           ***
     25!***    (THEIR STANDARD VALUES ARE 1.0 AND 0.96, RESPECTIVELY) ***
     26!***           (BETA MUST BE LESS THAN OR EQUAL TO 1)        ***
     27
     28!***    DTCRIT IS THE CRITICAL BUOYANCY (K) USED TO ADJUST THE ***
     29!***                 APPROACH TO QUASI-EQUILIBRIUM           ***
     30!***                     IT MUST BE LESS THAN 0              ***
    3131
    3232  include "cv3param.h"
    … …  
    4141
    4242  LOGICAL, SAVE :: first = .TRUE.
    43   !$OMP THREADPRIVATE(first)
    44 
    45   ! noff: integer limit for convection (nd-noff)
    46   ! minorig: First level of convection
    47 
    48   ! -- limit levels for convection:
     43!$OMP THREADPRIVATE(first)
     44
     45!glb noff: integer limit for convection (nd-noff)
     46! minorig: First level of convection
     47
     48! -- limit levels for convection:
    4949
    5050  noff = 1
    … …  
    5656  IF (first) THEN
    5757
    58     ! -- "microphysical" parameters:
     58!$OMP MASTER
     59! -- "microphysical" parameters:
    5960    sigdz = 0.01
    6061    spfac = 0.15
    6162    pbcrit = 150.0
    6263    ptcrit = 500.0
    63     ! IM beg: ajout fis. reglage ep
     64! IM beg: ajout fis. reglage ep
    6465    flag_epkeorig = 1
    6566    elcrit = 0.0003
    6667    tlcrit = -55.0
    67     ! IM lu dans physiq.def via conf_phys.F90     epmax  = 0.993
     68! IM lu dans physiq.def via conf_phys.F90     epmax  = 0.993
    6869
    6970    omtrain = 45.0 ! used also for snow (no disctinction rain/snow)
    7071
    71     ! -- misc:
     72! -- misc:
    7273
    7374    dtovsh = -0.2 ! dT for overshoot
    7475    dpbase = -40. ! definition cloud base (400m above LCL)
    75     ! cc      dttrig = 5.   ! (loose) condition for triggering
     76! cc      dttrig = 5.   ! (loose) condition for triggering
    7677    dttrig = 10. ! (loose) condition for triggering
    7778    flag_wb = 1
    7879    wbmax = 6. ! (m/s) adiab updraught speed at LFC (used in cv3p1_closure)
    7980
    80     ! -- rate of approach to quasi-equilibrium:
     81! -- rate of approach to quasi-equilibrium:
    8182
    8283    dtcrit = -2.0
    8384    tau = 8000.
    8485
    85     ! -- interface cloud parameterization:
     86! -- interface cloud parameterization:
    8687
    8788    delta = 0.01 ! cld
    8889
    89     ! -- interface with boundary-layer (gust factor): (sb)
     90! -- interface with boundary-layer (gust factor): (sb)
    9091
    9192    betad = 10.0 ! original value (from convect 4.3)
    9293
    93     OPEN (99, FILE='conv_param.data', STATUS='old', FORM='formatted', &
    94       ERR=9999)
     94    OPEN (99, FILE='conv_param.data', STATUS='old', FORM='formatted', ERR=9999)
    9595    READ (99, *, END=9998) dpbase
    9696    READ (99, *, END=9998) pbcrit
    … …  
    113113    WRITE (*, *) 'wbmax =', wbmax
    114114
    115     ! IM Lecture du fichier ep_param.data
     115! IM Lecture du fichier ep_param.data
    116116    OPEN (79, FILE='ep_param.data', STATUS='old', FORM='formatted', ERR=7999)
    117117    READ (79, *, END=7998) flag_epkeorig
    … …  
    124124    WRITE (*, *) 'elcrit=', elcrit
    125125    WRITE (*, *) 'tlcrit=', tlcrit
    126     ! IM end: ajout fis. reglage ep
     126! IM end: ajout fis. reglage ep
    127127
    128128    first = .FALSE.
    129   END IF
     129!$OMP END MASTER
     130
     131  END IF ! (first)
    130132
    131133  beta = 1.0 - delt/tau
    132134  alpha1 = 1.5E-3
    133   ! jyg    Correction bug alpha
     135!JYG    Correction bug alpha
    134136  alpha1 = alpha1*1.5
    135137  alpha = alpha1*delt/tau
    136   ! jyg    Bug
    137   ! cc increase alpha to compensate W decrease:
    138   ! c      alpha  = alpha*1.5
     138!JYG    Bug
     139! cc increase alpha to compensate W decrease:
     140! c      alpha  = alpha*1.5
    139141
    140142  RETURN
    141143END SUBROUTINE cv3_param
    142144
    143 SUBROUTINE cv3_prelim(len, nd, ndp1, t, q, p, ph, lv, lf, cpn, tv, gz, h, hm, &
    144     th)
     145SUBROUTINE cv3_prelim(len, nd, ndp1, t, q, p, ph, &
     146                      lv, lf, cpn, tv, gz, h, hm, th)
    145147  IMPLICIT NONE
    146148
    147   ! =====================================================================
    148   ! --- CALCULATE ARRAYS OF GEOPOTENTIAL, HEAT CAPACITY & STATIC ENERGY
    149   ! "ori": from convect4.3 (vectorized)
    150   ! "convect3": to be exactly consistent with convect3
    151   ! =====================================================================
    152 
    153   ! inputs:
     149! =====================================================================
     150! --- CALCULATE ARRAYS OF GEOPOTENTIAL, HEAT CAPACITY & STATIC ENERGY
     151! "ori": from convect4.3 (vectorized)
     152! "convect3": to be exactly consistent with convect3
     153! =====================================================================
     154
     155! inputs:
    154156  INTEGER len, nd, ndp1
    155157  REAL t(len, nd), q(len, nd), p(len, nd), ph(len, ndp1)
    156158
    157   ! outputs:
     159! outputs:
    158160  REAL lv(len, nd), lf(len, nd), cpn(len, nd), tv(len, nd)
    159161  REAL gz(len, nd), h(len, nd), hm(len, nd)
    160162  REAL th(len, nd)
    161163
    162   ! local variables:
     164! local variables:
    163165  INTEGER k, i
    164166  REAL rdcp
    … …  
    170172
    171173
    172   ! ori      do 110 k=1,nlp
    173   ! abderr     do 110 k=1,nl ! convect3
     174! ori      do 110 k=1,nlp
     175! abderr     do 110 k=1,nl ! convect3
    174176  DO k = 1, nlp
    175177
    176178    DO i = 1, len
    177       ! debug          lv(i,k)= lv0-clmcpv*(t(i,k)-t0)
     179! debug          lv(i,k)= lv0-clmcpv*(t(i,k)-t0)
    178180      lv(i, k) = lv0 - clmcpv*(t(i,k)-273.15)
    179181      lf(i, k) = lf0 - clmci*(t(i,k)-273.15)
    180182      cpn(i, k) = cpd*(1.0-q(i,k)) + cpv*q(i, k)
    181183      cpx(i, k) = cpd*(1.0-q(i,k)) + cl*q(i, k)
    182       ! ori          tv(i,k)=t(i,k)*(1.0+q(i,k)*epsim1)
     184! ori          tv(i,k)=t(i,k)*(1.0+q(i,k)*epsim1)
    183185      tv(i, k) = t(i, k)*(1.0+q(i,k)/eps-q(i,k))
    184186      rdcp = (rrd*(1.-q(i,k))+q(i,k)*rrv)/cpn(i, k)
    … …  
    187189  END DO
    188190
    189   ! gz = phi at the full levels (same as p).
     191! gz = phi at the full levels (same as p).
    190192
    191193  DO i = 1, len
    192194    gz(i, 1) = 0.0
    193195  END DO
    194   ! ori      do 140 k=2,nlp
     196! ori      do 140 k=2,nlp
    195197  DO k = 2, nl ! convect3
    196198    DO i = 1, len
    197       tvx = t(i, k)*(1.+q(i,k)/eps-q(i,k)) !convect3
    198       tvy = t(i, k-1)*(1.+q(i,k-1)/eps-q(i,k-1)) !convect3
    199       gz(i, k) = gz(i, k-1) + 0.5*rrd*(tvx+tvy) & !convect3
    200         *(p(i,k-1)-p(i,k))/ph(i, k) !convect3
    201 
    202       ! c        print *,' gz(',k,')',gz(i,k),' tvx',tvx,' tvy ',tvy
    203 
    204       ! ori         gz(i,k)=gz(i,k-1)+hrd*(tv(i,k-1)+tv(i,k))
    205       ! ori    &         *(p(i,k-1)-p(i,k))/ph(i,k)
    206     END DO
    207   END DO
    208 
    209   ! h  = phi + cpT (dry static energy).
    210   ! hm = phi + cp(T-Tbase)+Lq
    211 
    212   ! ori      do 170 k=1,nlp
     199      tvx = t(i, k)*(1.+q(i,k)/eps-q(i,k))         !convect3
     200      tvy = t(i, k-1)*(1.+q(i,k-1)/eps-q(i,k-1))   !convect3
     201      gz(i, k) = gz(i, k-1) + 0.5*rrd*(tvx+tvy)* & !convect3
     202                 (p(i,k-1)-p(i,k))/ph(i, k)        !convect3
     203
     204! c        print *,' gz(',k,')',gz(i,k),' tvx',tvx,' tvy ',tvy
     205
     206! ori         gz(i,k)=gz(i,k-1)+hrd*(tv(i,k-1)+tv(i,k))
     207! ori    &         *(p(i,k-1)-p(i,k))/ph(i,k)
     208    END DO
     209  END DO
     210
     211! h  = phi + cpT (dry static energy).
     212! hm = phi + cp(T-Tbase)+Lq
     213
     214! ori      do 170 k=1,nlp
    213215  DO k = 1, nl ! convect3
    214216    DO i = 1, len
    … …  
    221223END SUBROUTINE cv3_prelim
    222224
    223 SUBROUTINE cv3_feed(len, nd, t, q, u, v, p, ph, hm, gz, p1feed, p2feed, wght, &
    224     wghti, tnk, thnk, qnk, qsnk, unk, vnk, cpnk, hnk, nk, icb, icbmax, iflag, &
    225     gznk, plcl)
     225SUBROUTINE cv3_feed(len, nd, ok_conserv_q, &
     226                    t, q, u, v, p, ph, hm, gz, &
     227                    p1feed, p2feed, wght, &
     228                    wghti, tnk, thnk, qnk, qsnk, unk, vnk, &
     229                    cpnk, hnk, nk, icb, icbmax, iflag, gznk, plcl)
    226230  IMPLICIT NONE
    227231
    228   ! ================================================================
    229   ! Purpose: CONVECTIVE FEED
    230 
    231   ! Main differences with cv_feed:
    232   ! - ph added in input
    233   ! - here, nk(i)=minorig
    234   ! - icb defined differently (plcl compared with ph instead of p)
    235 
    236   ! Main differences with convect3:
    237   ! - we do not compute dplcldt and dplcldr of CLIFT anymore
    238   ! - values iflag different (but tests identical)
    239   ! - A,B explicitely defined (!...)
    240   ! ================================================================
     232! ================================================================
     233! Purpose: CONVECTIVE FEED
     234
     235! Main differences with cv_feed:
     236! - ph added in input
     237! - here, nk(i)=minorig
     238! - icb defined differently (plcl compared with ph instead of p)
     239
     240! Main differences with convect3:
     241! - we do not compute dplcldt and dplcldr of CLIFT anymore
     242! - values iflag different (but tests identical)
     243! - A,B explicitely defined (!...)
     244! ================================================================
    241245
    242246  include "cv3param.h"
    243247  include "cvthermo.h"
    244248
    245   ! inputs:
     249!inputs:
    246250  INTEGER len, nd
     251  LOGICAL ok_conserv_q
    247252  REAL t(len, nd), q(len, nd), p(len, nd)
    248253  REAL u(len, nd), v(len, nd)
    … …  
    250255  REAL ph(len, nd+1)
    251256  REAL p1feed(len)
    252   ! ,  wght(len)
     257! ,  wght(len)
    253258  REAL wght(nd)
    254   ! input-output
     259!input-output
    255260  REAL p2feed(len)
    256   ! outputs:
     261!outputs:
    257262  INTEGER iflag(len), nk(len), icb(len), icbmax
    258   ! real   wghti(len)
     263!  real   wghti(len)
    259264  REAL wghti(len, nd)
    260265  REAL tnk(len), thnk(len), qnk(len), qsnk(len)
    … …  
    263268  REAL plcl(len)
    264269
    265   ! local variables:
     270!local variables:
    266271  INTEGER i, k, iter, niter
    267272  INTEGER ihmin(len)
    … …  
    269274  REAL pup(len), plo(len), pfeed(len)
    270275  REAL plclup(len), plcllo(len), plclfeed(len)
     276  REAL pfeedmin(len)
    271277  REAL posit(len)
    272278  LOGICAL nocond(len)
    273279
    274   ! -------------------------------------------------------------------
    275   ! --- Origin level of ascending parcels for convect3:
    276   ! -------------------------------------------------------------------
     280!jyg20140217<
     281  INTEGER iostat
     282  LOGICAL, SAVE :: first
     283  LOGICAL, SAVE :: ok_new_feed
     284  REAL, SAVE :: dp_lcl_feed
     285!$OMP THREADPRIVATE (first,ok_new_feed,dp_lcl_feed)
     286  DATA first/.TRUE./
     287  DATA dp_lcl_feed/2./
     288
     289  IF (first) THEN
     290!$OMP MASTER
     291    ok_new_feed = ok_conserv_q
     292    OPEN (98, FILE='cv3feed_param.data', STATUS='old', FORM='formatted', IOSTAT=iostat)
     293    IF (iostat==0) THEN
     294      READ (98, *, END=998) ok_new_feed
     295998   CONTINUE
     296      CLOSE (98)
     297    END IF
     298    PRINT *, ' ok_new_feed: ', ok_new_feed
     299    first = .FALSE.
     300!$OMP END MASTER
     301  END IF
     302!jyg>
     303! -------------------------------------------------------------------
     304! --- Origin level of ascending parcels for convect3:
     305! -------------------------------------------------------------------
    277306
    278307  DO i = 1, len
    … …  
    281310  END DO
    282311
    283   ! -------------------------------------------------------------------
    284   ! --- Adjust feeding layer thickness so that lifting up to the top of
    285   ! --- the feeding layer does not induce condensation (i.e. so that
    286   ! --- plcl < p2feed).
    287   ! --- Method : iterative secant method.
    288   ! -------------------------------------------------------------------
    289 
    290   ! 1- First bracketing of the solution : ph(nk+1), p2feed
    291 
    292   ! 1.a- LCL associated to p2feed
     312! -------------------------------------------------------------------
     313! --- Adjust feeding layer thickness so that lifting up to the top of
     314! --- the feeding layer does not induce condensation (i.e. so that
     315! --- plcl < p2feed).
     316! --- Method : iterative secant method.
     317! -------------------------------------------------------------------
     318
     319! 1- First bracketing of the solution : ph(nk+1), p2feed
     320
     321! 1.a- LCL associated with p2feed
    293322  DO i = 1, len
    294323    pup(i) = p2feed(i)
    295324  END DO
    296   CALL cv3_vertmix(len, nd, iflag, p1feed, pup, p, ph, t, q, u, v, wght, &
    297     wghti, nk, tnk, thnk, qnk, qsnk, unk, vnk, plclup)
    298   ! 1.b- LCL associated to ph(nk+1)
     325  CALL cv3_vertmix(len, nd, iflag, p1feed, pup, p, ph, &
     326                   t, q, u, v, wght, &
     327                   wghti, nk, tnk, thnk, qnk, qsnk, unk, vnk, plclup)
     328! 1.b- LCL associated with ph(nk+1)
    299329  DO i = 1, len
    300330    plo(i) = ph(i, nk(i)+1)
    301331  END DO
    302   CALL cv3_vertmix(len, nd, iflag, p1feed, plo, p, ph, t, q, u, v, wght, &
    303     wghti, nk, tnk, thnk, qnk, qsnk, unk, vnk, plcllo)
    304   ! 2- Iterations
     332  CALL cv3_vertmix(len, nd, iflag, p1feed, plo, p, ph, &
     333                   t, q, u, v, wght, &
     334                   wghti, nk, tnk, thnk, qnk, qsnk, unk, vnk, plcllo)
     335! 2- Iterations
    305336  niter = 5
    306337  DO iter = 1, niter
    … …  
    314345        pfeed(i) = pup(i)
    315346      ELSE
    316         pfeed(i) = (pup(i)*(plo(i)-plcllo(i))+plo(i)*(plclup(i)-pup(i)))/ &
    317           (plo(i)-plcllo(i)+plclup(i)-pup(i))
     347!JYG20140217<
     348        IF (ok_new_feed) THEN
     349          pfeed(i) = (pup(i)*(plo(i)-plcllo(i)-dp_lcl_feed)+  &
     350                      plo(i)*(plclup(i)-pup(i)+dp_lcl_feed))/ &
     351                     (plo(i)-plcllo(i)+plclup(i)-pup(i))
     352        ELSE
     353          pfeed(i) = (pup(i)*(plo(i)-plcllo(i))+  &
     354                      plo(i)*(plclup(i)-pup(i)))/ &
     355                     (plo(i)-plcllo(i)+plclup(i)-pup(i))
     356        END IF
     357!JYG>
    318358      END IF
    319359    END DO
    320     CALL cv3_vertmix(len, nd, iflag, p1feed, pfeed, p, ph, t, q, u, v, wght, &
    321       wghti, nk, tnk, thnk, qnk, qsnk, unk, vnk, plclfeed)
     360!jyg20140217<
     361! For the last iteration, make sure that the top of the feeding layer
     362! and LCL are not in the same layer:
     363    IF (ok_new_feed) THEN
     364      IF (iter==niter) THEN
     365        DO k = minorig, nd
     366          DO i = 1, len
     367            IF (ph(i,k)>=plclfeed(i)) pfeedmin(i) = ph(i, k)
     368          END DO
     369        END DO
     370        DO i = 1, len
     371          pfeed(i) = max(pfeedmin(i), pfeed(i))
     372        END DO
     373      END IF
     374    END IF
     375!jyg>
     376
     377    CALL cv3_vertmix(len, nd, iflag, p1feed, pfeed, p, ph, &
     378                   t, q, u, v, wght, &
     379                   wghti, nk, tnk, thnk, qnk, qsnk, unk, vnk, plclfeed)
     380!jyg20140217<
     381    IF (ok_new_feed) THEN
     382      DO i = 1, len
     383        posit(i) = (sign(1.,plclfeed(i)-pfeed(i)+dp_lcl_feed)+1.)*0.5
     384        IF (plclfeed(i)-pfeed(i)+dp_lcl_feed==0.) posit(i) = 1.
     385      END DO
     386    ELSE
     387      DO i = 1, len
     388        posit(i) = (sign(1.,plclfeed(i)-pfeed(i))+1.)*0.5
     389        IF (plclfeed(i)==pfeed(i)) posit(i) = 1.
     390      END DO
     391    END IF
     392!jyg>
    322393    DO i = 1, len
    323       posit(i) = (sign(1.,plclfeed(i)-pfeed(i))+1.)*0.5
    324       IF (plclfeed(i)==pfeed(i)) posit(i) = 1.
    325       ! - posit = 1 when lcl is below top of feeding layer (plclfeed>pfeed)
    326       ! -               => pup=pfeed
    327       ! - posit = 0 when lcl is above top of feeding layer (plclfeed<pfeed)
    328       ! -               => plo=pfeed
     394! - posit = 1 when lcl is below top of feeding layer (plclfeed>pfeed)
     395! -               => pup=pfeed
     396! - posit = 0 when lcl is above top of feeding layer (plclfeed<pfeed)
     397! -               => plo=pfeed
    329398      pup(i) = posit(i)*pfeed(i) + (1.-posit(i))*pup(i)
    330399      plo(i) = (1.-posit(i))*pfeed(i) + posit(i)*plo(i)
    … …  
    343412  END DO
    344413
    345   ! -------------------------------------------------------------------
    346   ! --- Check whether parcel level temperature and specific humidity
    347   ! --- are reasonable
    348   ! -------------------------------------------------------------------
     414! -------------------------------------------------------------------
     415! --- Check whether parcel level temperature and specific humidity
     416! --- are reasonable
     417! -------------------------------------------------------------------
    349418  DO i = 1, len
    350419    IF (((tnk(i)<250.0) .OR. (qnk(i)<=0.0)) .AND. (iflag(i)==0)) iflag(i) = 7
    351420  END DO
    352421
    353   ! -------------------------------------------------------------------
    354   ! --- Calculate first level above lcl (=icb)
    355   ! -------------------------------------------------------------------
    356 
    357   ! @      do 270 i=1,len
    358   ! @       icb(i)=nlm
    359   ! @ 270  continue
    360   ! @c
    361   ! @      do 290 k=minorig,nl
    362   ! @        do 280 i=1,len
    363   ! @          if((k.ge.(nk(i)+1)).and.(p(i,k).lt.plcl(i)))
    364   ! @     &    icb(i)=min(icb(i),k)
    365   ! @ 280    continue
    366   ! @ 290  continue
    367   ! @c
    368   ! @      do 300 i=1,len
    369   ! @        if((icb(i).ge.nlm).and.(iflag(i).eq.0))iflag(i)=9
    370   ! @ 300  continue
     422! -------------------------------------------------------------------
     423! --- Calculate first level above lcl (=icb)
     424! -------------------------------------------------------------------
     425
     426!@      do 270 i=1,len
     427!@       icb(i)=nlm
     428!@ 270  continue
     429!@c
     430!@      do 290 k=minorig,nl
     431!@        do 280 i=1,len
     432!@          if((k.ge.(nk(i)+1)).and.(p(i,k).lt.plcl(i)))
     433!@     &    icb(i)=min(icb(i),k)
     434!@ 280    continue
     435!@ 290  continue
     436!@c
     437!@      do 300 i=1,len
     438!@        if((icb(i).ge.nlm).and.(iflag(i).eq.0))iflag(i)=9
     439!@ 300  continue
    371440
    372441  DO i = 1, len
    … …  
    374443  END DO
    375444
    376   ! la modification consiste a comparer plcl a ph et non a p:
    377   ! icb est defini par :  ph(icb)<plcl<ph(icb-1)
    378   ! @      do 290 k=minorig,nl
     445! la modification consiste a comparer plcl a ph et non a p:
     446! icb est defini par :  ph(icb)<plcl<ph(icb-1)
     447!@      do 290 k=minorig,nl
    379448  DO k = 3, nl - 1 ! modif pour que icb soit sup/egal a 2
    380449    DO i = 1, len
    … …  
    384453
    385454
    386   ! print*,'icb dans cv3_feed '
    387   ! write(*,'(64i2)') icb(2:len-1)
    388   ! call dump2d(64,43,'plcl dans cv3_feed ',plcl(2:len-1))
     455! print*,'icb dans cv3_feed '
     456! write(*,'(64i2)') icb(2:len-1)
     457! call dump2d(64,43,'plcl dans cv3_feed ',plcl(2:len-1))
    389458
    390459  DO i = 1, len
    391     ! @        if((icb(i).ge.nlm).and.(iflag(i).eq.0))iflag(i)=9
     460!@        if((icb(i).ge.nlm).and.(iflag(i).eq.0))iflag(i)=9
    392461    IF ((icb(i)==nlm) .AND. (iflag(i)==0)) iflag(i) = 9
    393462  END DO
    … …  
    397466  END DO
    398467
    399   ! Compute icbmax.
     468! Compute icbmax.
    400469
    401470  icbmax = 2
    402471  DO i = 1, len
    403     ! !        icbmax=max(icbmax,icb(i))
    404     IF (iflag(i)<7) icbmax = max(icbmax, icb(i)) ! sb Jun7th02
     472!!        icbmax=max(icbmax,icb(i))
     473    IF (iflag(i)<7) icbmax = max(icbmax, icb(i))     ! sb Jun7th02
    405474  END DO
    406475
    … …  
    409478
    410479SUBROUTINE cv3_undilute1(len, nd, t, qs, gz, plcl, p, icb, tnk, qnk, gznk, &
    411     tp, tvp, clw, icbs)
     480                         tp, tvp, clw, icbs)
    412481  IMPLICIT NONE
    413482
    414   ! ----------------------------------------------------------------
    415   ! Equivalent de TLIFT entre NK et ICB+1 inclus
    416 
    417   ! Differences with convect4:
    418   ! - specify plcl in input
    419   ! - icbs is the first level above LCL (may differ from icb)
    420   ! - in the iterations, used x(icbs) instead x(icb)
    421   ! - many minor differences in the iterations
    422   ! - tvp is computed in only one time
    423   ! - icbs: first level above Plcl (IMIN de TLIFT) in output
    424   ! - if icbs=icb, compute also tp(icb+1),tvp(icb+1) & clw(icb+1)
    425   ! ----------------------------------------------------------------
     483! ----------------------------------------------------------------
     484! Equivalent de TLIFT entre NK et ICB+1 inclus
     485
     486! Differences with convect4:
     487!    - specify plcl in input
     488!    - icbs is the first level above LCL (may differ from icb)
     489!    - in the iterations, used x(icbs) instead x(icb)
     490!    - many minor differences in the iterations
     491!    - tvp is computed in only one time
     492!    - icbs: first level above Plcl (IMIN de TLIFT) in output
     493!    - if icbs=icb, compute also tp(icb+1),tvp(icb+1) & clw(icb+1)
     494! ----------------------------------------------------------------
    426495
    427496  include "cvthermo.h"
    428497  include "cv3param.h"
    429498
    430   ! inputs:
     499! inputs:
    431500  INTEGER len, nd
    432501  INTEGER icb(len)
    … …  
    436505  REAL plcl(len) ! convect3
    437506
    438   ! outputs:
     507! outputs:
    439508  REAL tp(len, nd), tvp(len, nd), clw(len, nd)
    440509
    441   ! local variables:
     510! local variables:
    442511  INTEGER i, k
    443512  INTEGER icb1(len), icbs(len), icbsmax2 ! convect3
    … …  
    448517  REAL cpinv(len) ! convect3
    449518
    450   ! -------------------------------------------------------------------
    451   ! --- Calculates the lifted parcel virtual temperature at nk,
    452   ! --- the actual temperature, and the adiabatic
    453   ! --- liquid water content. The procedure is to solve the equation.
    454   ! cp*tp+L*qp+phi=cp*tnk+L*qnk+gznk.
    455   ! -------------------------------------------------------------------
    456 
    457 
    458   ! ***  Calculate certain parcel quantities, including static energy   ***
     519! -------------------------------------------------------------------
     520! --- Calculates the lifted parcel virtual temperature at nk,
     521! --- the actual temperature, and the adiabatic
     522! --- liquid water content. The procedure is to solve the equation.
     523!    cp*tp+L*qp+phi=cp*tnk+L*qnk+gznk.
     524! -------------------------------------------------------------------
     525
     526
     527! ***  Calculate certain parcel quantities, including static energy   ***
    459528
    460529  DO i = 1, len
    461     ah0(i) = (cpd*(1.-qnk(i))+cl*qnk(i))*tnk(i) + qnk(i)*(lv0-clmcpv*(tnk(i)- &
    462       273.15)) + gznk(i)
     530    ah0(i) = (cpd*(1.-qnk(i))+cl*qnk(i))*tnk(i) + qnk(i)*(lv0-clmcpv*(tnk(i)-273.15)) + gznk(i)
    463531    cpp(i) = cpd*(1.-qnk(i)) + qnk(i)*cpv
    464532    cpinv(i) = 1./cpp(i)
    465533  END DO
    466534
    467   ! ***   Calculate lifted parcel quantities below cloud base   ***
     535! ***   Calculate lifted parcel quantities below cloud base   ***
     536
     537  DO i = 1, len                                           !convect3
     538    icb1(i) = max(icb(i), 2)                              !convect3
     539    icb1(i) = min(icb(i), nl)                             !convect3
     540! if icb is below LCL, start loop at ICB+1:
     541! (icbs est le premier niveau au-dessus du LCL)
     542    icbs(i) = icb1(i)                                     !convect3
     543    IF (plcl(i)<p(i,icb1(i))) THEN
     544      icbs(i) = min(icbs(i)+1, nl)                        !convect3
     545    END IF
     546  END DO                                                  !convect3
    468547
    469548  DO i = 1, len !convect3
    470     icb1(i) = max(icb(i), 2) !convect3
    471     icb1(i) = min(icb(i), nl) !convect3
    472     ! if icb is below LCL, start loop at ICB+1:
    473     ! (icbs est le premier niveau au-dessus du LCL)
    474     icbs(i) = icb1(i) !convect3
    475     IF (plcl(i)<p(i,icb1(i))) THEN
    476       icbs(i) = min(icbs(i)+1, nl) !convect3
    477     END IF
     549    ticb(i) = t(i, icbs(i))                               !convect3
     550    gzicb(i) = gz(i, icbs(i))                             !convect3
     551    qsicb(i) = qs(i, icbs(i))                             !convect3
    478552  END DO !convect3
    479553
    480   DO i = 1, len !convect3
    481     ticb(i) = t(i, icbs(i)) !convect3
    482     gzicb(i) = gz(i, icbs(i)) !convect3
    483     qsicb(i) = qs(i, icbs(i)) !convect3
    484   END DO !convect3
    485 
    486 
    487   ! Re-compute icbsmax (icbsmax2):        !convect3
    488   ! !convect3
    489   icbsmax2 = 2 !convect3
    490   DO i = 1, len !convect3
    491     icbsmax2 = max(icbsmax2, icbs(i)) !convect3
    492   END DO !convect3
    493 
    494   ! initialization outputs:
    495 
    496   DO k = 1, icbsmax2 ! convect3
    497     DO i = 1, len ! convect3
    498       tp(i, k) = 0.0 ! convect3
    499       tvp(i, k) = 0.0 ! convect3
    500       clw(i, k) = 0.0 ! convect3
    501     END DO ! convect3
    502   END DO ! convect3
    503 
    504   ! tp and tvp below cloud base:
     554
     555! Re-compute icbsmax (icbsmax2):                          !convect3
     556!                                                         !convect3
     557  icbsmax2 = 2                                            !convect3
     558  DO i = 1, len                                           !convect3
     559    icbsmax2 = max(icbsmax2, icbs(i))                     !convect3
     560  END DO                                                  !convect3
     561
     562! initialization outputs:
     563
     564  DO k = 1, icbsmax2                                      ! convect3
     565    DO i = 1, len                                         ! convect3
     566      tp(i, k) = 0.0                                      ! convect3
     567      tvp(i, k) = 0.0                                     ! convect3
     568      clw(i, k) = 0.0                                     ! convect3
     569    END DO                                                ! convect3
     570  END DO                                                  ! convect3
     571
     572! tp and tvp below cloud base:
    505573
    506574  DO k = minorig, icbsmax2 - 1
    507575    DO i = 1, len
    508576      tp(i, k) = tnk(i) - (gz(i,k)-gznk(i))*cpinv(i)
    509       tvp(i, k) = tp(i, k)*(1.+qnk(i)/eps-qnk(i)) !whole thing (convect3)
    510     END DO
    511   END DO
    512 
    513   ! ***  Find lifted parcel quantities above cloud base    ***
     577      tvp(i, k) = tp(i, k)*(1.+qnk(i)/eps-qnk(i))        !whole thing (convect3)
     578    END DO
     579  END DO
     580
     581! ***  Find lifted parcel quantities above cloud base    ***
    514582
    515583  DO i = 1, len
    516584    tg = ticb(i)
    517     ! ori         qg=qs(i,icb(i))
     585! ori         qg=qs(i,icb(i))
    518586    qg = qsicb(i) ! convect3
    519     ! debug         alv=lv0-clmcpv*(ticb(i)-t0)
     587! debug         alv=lv0-clmcpv*(ticb(i)-t0)
    520588    alv = lv0 - clmcpv*(ticb(i)-273.15)
    521589
    522     ! First iteration.
    523 
    524     ! ori          s=cpd+alv*alv*qg/(rrv*ticb(i)*ticb(i))
    525     s = cpd*(1.-qnk(i)) + cl*qnk(i) & ! convect3
    526       +alv*alv*qg/(rrv*ticb(i)*ticb(i)) ! convect3
     590! First iteration.
     591
     592! ori          s=cpd+alv*alv*qg/(rrv*ticb(i)*ticb(i))
     593    s = cpd*(1.-qnk(i)) + cl*qnk(i) + &                  ! convect3
     594        alv*alv*qg/(rrv*ticb(i)*ticb(i))                  ! convect3
    527595    s = 1./s
    528     ! ori          ahg=cpd*tg+(cl-cpd)*qnk(i)*ticb(i)+alv*qg+gzicb(i)
     596! ori          ahg=cpd*tg+(cl-cpd)*qnk(i)*ticb(i)+alv*qg+gzicb(i)
    529597    ahg = cpd*tg + (cl-cpd)*qnk(i)*tg + alv*qg + gzicb(i) ! convect3
    530598    tg = tg + s*(ah0(i)-ahg)
    531     ! ori          tg=max(tg,35.0)
    532     ! debug          tc=tg-t0
     599! ori          tg=max(tg,35.0)
     600! debug          tc=tg-t0
    533601    tc = tg - 273.15
    534602    denom = 243.5 + tc
    535603    denom = max(denom, 1.0) ! convect3
    536     ! ori          if(tc.ge.0.0)then
     604! ori          if(tc.ge.0.0)then
    537605    es = 6.112*exp(17.67*tc/denom)
    538     ! ori          else
    539     ! ori           es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
    540     ! ori          endif
    541     ! ori          qg=eps*es/(p(i,icb(i))-es*(1.-eps))
     606! ori          else
     607! ori           es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
     608! ori          endif
     609! ori          qg=eps*es/(p(i,icb(i))-es*(1.-eps))
    542610    qg = eps*es/(p(i,icbs(i))-es*(1.-eps))
    543611
    544     ! Second iteration.
    545 
    546 
    547     ! ori          s=cpd+alv*alv*qg/(rrv*ticb(i)*ticb(i))
    548     ! ori          s=1./s
    549     ! ori          ahg=cpd*tg+(cl-cpd)*qnk(i)*ticb(i)+alv*qg+gzicb(i)
     612! Second iteration.
     613
     614
     615! ori          s=cpd+alv*alv*qg/(rrv*ticb(i)*ticb(i))
     616! ori          s=1./s
     617! ori          ahg=cpd*tg+(cl-cpd)*qnk(i)*ticb(i)+alv*qg+gzicb(i)
    550618    ahg = cpd*tg + (cl-cpd)*qnk(i)*tg + alv*qg + gzicb(i) ! convect3
    551619    tg = tg + s*(ah0(i)-ahg)
    552     ! ori          tg=max(tg,35.0)
    553     ! debug          tc=tg-t0
     620! ori          tg=max(tg,35.0)
     621! debug          tc=tg-t0
    554622    tc = tg - 273.15
    555623    denom = 243.5 + tc
    556     denom = max(denom, 1.0) ! convect3
    557     ! ori          if(tc.ge.0.0)then
     624    denom = max(denom, 1.0)                               ! convect3
     625! ori          if(tc.ge.0.0)then
    558626    es = 6.112*exp(17.67*tc/denom)
    559     ! ori          else
    560     ! ori           es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
    561     ! ori          end if
    562     ! ori          qg=eps*es/(p(i,icb(i))-es*(1.-eps))
     627! ori          else
     628! ori           es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
     629! ori          end if
     630! ori          qg=eps*es/(p(i,icb(i))-es*(1.-eps))
    563631    qg = eps*es/(p(i,icbs(i))-es*(1.-eps))
    564632
    565633    alv = lv0 - clmcpv*(ticb(i)-273.15)
    566634
    567     ! ori c approximation here:
    568     ! ori         tp(i,icb(i))=(ah0(i)-(cl-cpd)*qnk(i)*ticb(i)
    569     ! ori     &   -gz(i,icb(i))-alv*qg)/cpd
    570 
    571     ! convect3: no approximation:
     635! ori c approximation here:
     636! ori         tp(i,icb(i))=(ah0(i)-(cl-cpd)*qnk(i)*ticb(i)
     637! ori     &   -gz(i,icb(i))-alv*qg)/cpd
     638
     639! convect3: no approximation:
    572640    tp(i, icbs(i)) = (ah0(i)-gz(i,icbs(i))-alv*qg)/(cpd+(cl-cpd)*qnk(i))
    573641
    574     ! ori         clw(i,icb(i))=qnk(i)-qg
    575     ! ori         clw(i,icb(i))=max(0.0,clw(i,icb(i)))
     642! ori         clw(i,icb(i))=qnk(i)-qg
     643! ori         clw(i,icb(i))=max(0.0,clw(i,icb(i)))
    576644    clw(i, icbs(i)) = qnk(i) - qg
    577645    clw(i, icbs(i)) = max(0.0, clw(i,icbs(i)))
    578646
    579647    rg = qg/(1.-qnk(i))
    580     ! ori         tvp(i,icb(i))=tp(i,icb(i))*(1.+rg*epsi)
    581     ! convect3: (qg utilise au lieu du vrai mixing ratio rg)
    582     tvp(i, icbs(i)) = tp(i, icbs(i))*(1.+qg/eps-qnk(i)) !whole thing
    583 
    584   END DO
    585 
    586   ! ori      do 380 k=minorig,icbsmax2
    587   ! ori       do 370 i=1,len
    588   ! ori         tvp(i,k)=tvp(i,k)-tp(i,k)*qnk(i)
    589   ! ori 370   continue
    590   ! ori 380  continue
    591 
    592 
    593   ! -- The following is only for convect3:
    594 
    595   ! * icbs is the first level above the LCL:
    596   ! if plcl<p(icb), then icbs=icb+1
    597   ! if plcl>p(icb), then icbs=icb
    598 
    599   ! * the routine above computes tvp from minorig to icbs (included).
    600 
    601   ! * to compute buoybase (in cv3_trigger.F), both tvp(icb) and tvp(icb+1)
    602   ! must be known. This is the case if icbs=icb+1, but not if icbs=icb.
    603 
    604   ! * therefore, in the case icbs=icb, we compute tvp at level icb+1
    605   ! (tvp at other levels will be computed in cv3_undilute2.F)
     648! ori         tvp(i,icb(i))=tp(i,icb(i))*(1.+rg*epsi)
     649! convect3: (qg utilise au lieu du vrai mixing ratio rg)
     650    tvp(i, icbs(i)) = tp(i, icbs(i))*(1.+qg/eps-qnk(i))   !whole thing
     651
     652  END DO
     653
     654! ori      do 380 k=minorig,icbsmax2
     655! ori       do 370 i=1,len
     656! ori         tvp(i,k)=tvp(i,k)-tp(i,k)*qnk(i)
     657! ori 370   continue
     658! ori 380  continue
     659
     660
     661! -- The following is only for convect3:
     662
     663! * icbs is the first level above the LCL:
     664! if plcl<p(icb), then icbs=icb+1
     665! if plcl>p(icb), then icbs=icb
     666
     667! * the routine above computes tvp from minorig to icbs (included).
     668
     669! * to compute buoybase (in cv3_trigger.F), both tvp(icb) and tvp(icb+1)
     670! must be known. This is the case if icbs=icb+1, but not if icbs=icb.
     671
     672! * therefore, in the case icbs=icb, we compute tvp at level icb+1
     673! (tvp at other levels will be computed in cv3_undilute2.F)
    606674
    607675
    … …  
    615683    tg = ticb(i)
    616684    qg = qsicb(i) ! convect3
    617     ! debug         alv=lv0-clmcpv*(ticb(i)-t0)
     685! debug         alv=lv0-clmcpv*(ticb(i)-t0)
    618686    alv = lv0 - clmcpv*(ticb(i)-273.15)
    619687
    620     ! First iteration.
    621 
    622     ! ori          s=cpd+alv*alv*qg/(rrv*ticb(i)*ticb(i))
    623     s = cpd*(1.-qnk(i)) + cl*qnk(i) & ! convect3
    624       +alv*alv*qg/(rrv*ticb(i)*ticb(i)) ! convect3
     688! First iteration.
     689
     690! ori          s=cpd+alv*alv*qg/(rrv*ticb(i)*ticb(i))
     691    s = cpd*(1.-qnk(i)) + cl*qnk(i) &                         ! convect3
     692      +alv*alv*qg/(rrv*ticb(i)*ticb(i))                       ! convect3
    625693    s = 1./s
    626     ! ori          ahg=cpd*tg+(cl-cpd)*qnk(i)*ticb(i)+alv*qg+gzicb(i)
    627     ahg = cpd*tg + (cl-cpd)*qnk(i)*tg + alv*qg + gzicb(i) ! convect3
     694! ori          ahg=cpd*tg+(cl-cpd)*qnk(i)*ticb(i)+alv*qg+gzicb(i)
     695    ahg = cpd*tg + (cl-cpd)*qnk(i)*tg + alv*qg + gzicb(i)     ! convect3
    628696    tg = tg + s*(ah0(i)-ahg)
    629     ! ori          tg=max(tg,35.0)
    630     ! debug          tc=tg-t0
     697! ori          tg=max(tg,35.0)
     698! debug          tc=tg-t0
    631699    tc = tg - 273.15
    632700    denom = 243.5 + tc
    633     denom = max(denom, 1.0) ! convect3
    634     ! ori          if(tc.ge.0.0)then
     701    denom = max(denom, 1.0)                                   ! convect3
     702! ori          if(tc.ge.0.0)then
    635703    es = 6.112*exp(17.67*tc/denom)
    636     ! ori          else
    637     ! ori           es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
    638     ! ori          endif
    639     ! ori          qg=eps*es/(p(i,icb(i))-es*(1.-eps))
     704! ori          else
     705! ori           es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
     706! ori          endif
     707! ori          qg=eps*es/(p(i,icb(i))-es*(1.-eps))
    640708    qg = eps*es/(p(i,icb(i)+1)-es*(1.-eps))
    641709
    642     ! Second iteration.
    643 
    644 
    645     ! ori          s=cpd+alv*alv*qg/(rrv*ticb(i)*ticb(i))
    646     ! ori          s=1./s
    647     ! ori          ahg=cpd*tg+(cl-cpd)*qnk(i)*ticb(i)+alv*qg+gzicb(i)
    648     ahg = cpd*tg + (cl-cpd)*qnk(i)*tg + alv*qg + gzicb(i) ! convect3
     710! Second iteration.
     711
     712
     713! ori          s=cpd+alv*alv*qg/(rrv*ticb(i)*ticb(i))
     714! ori          s=1./s
     715! ori          ahg=cpd*tg+(cl-cpd)*qnk(i)*ticb(i)+alv*qg+gzicb(i)
     716    ahg = cpd*tg + (cl-cpd)*qnk(i)*tg + alv*qg + gzicb(i)     ! convect3
    649717    tg = tg + s*(ah0(i)-ahg)
    650     ! ori          tg=max(tg,35.0)
    651     ! debug          tc=tg-t0
     718! ori          tg=max(tg,35.0)
     719! debug          tc=tg-t0
    652720    tc = tg - 273.15
    653721    denom = 243.5 + tc
    654     denom = max(denom, 1.0) ! convect3
    655     ! ori          if(tc.ge.0.0)then
     722    denom = max(denom, 1.0)                                   ! convect3
     723! ori          if(tc.ge.0.0)then
    656724    es = 6.112*exp(17.67*tc/denom)
    657     ! ori          else
    658     ! ori           es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
    659     ! ori          end if
    660     ! ori          qg=eps*es/(p(i,icb(i))-es*(1.-eps))
     725! ori          else
     726! ori           es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
     727! ori          end if
     728! ori          qg=eps*es/(p(i,icb(i))-es*(1.-eps))
    661729    qg = eps*es/(p(i,icb(i)+1)-es*(1.-eps))
    662730
    663731    alv = lv0 - clmcpv*(ticb(i)-273.15)
    664732
    665     ! ori c approximation here:
    666     ! ori         tp(i,icb(i))=(ah0(i)-(cl-cpd)*qnk(i)*ticb(i)
    667     ! ori     &   -gz(i,icb(i))-alv*qg)/cpd
    668 
    669     ! convect3: no approximation:
     733! ori c approximation here:
     734! ori         tp(i,icb(i))=(ah0(i)-(cl-cpd)*qnk(i)*ticb(i)
     735! ori     &   -gz(i,icb(i))-alv*qg)/cpd
     736
     737! convect3: no approximation:
    670738    tp(i, icb(i)+1) = (ah0(i)-gz(i,icb(i)+1)-alv*qg)/(cpd+(cl-cpd)*qnk(i))
    671739
    672     ! ori         clw(i,icb(i))=qnk(i)-qg
    673     ! ori         clw(i,icb(i))=max(0.0,clw(i,icb(i)))
     740! ori         clw(i,icb(i))=qnk(i)-qg
     741! ori         clw(i,icb(i))=max(0.0,clw(i,icb(i)))
    674742    clw(i, icb(i)+1) = qnk(i) - qg
    675743    clw(i, icb(i)+1) = max(0.0, clw(i,icb(i)+1))
    676744
    677745    rg = qg/(1.-qnk(i))
    678     ! ori         tvp(i,icb(i))=tp(i,icb(i))*(1.+rg*epsi)
    679     ! convect3: (qg utilise au lieu du vrai mixing ratio rg)
    680     tvp(i, icb(i)+1) = tp(i, icb(i)+1)*(1.+qg/eps-qnk(i)) !whole thing
     746! ori         tvp(i,icb(i))=tp(i,icb(i))*(1.+rg*epsi)
     747! convect3: (qg utilise au lieu du vrai mixing ratio rg)
     748    tvp(i, icb(i)+1) = tp(i, icb(i)+1)*(1.+qg/eps-qnk(i))     !whole thing
    681749
    682750  END DO
    … …  
    685753END SUBROUTINE cv3_undilute1
    686754
    687 SUBROUTINE cv3_trigger(len, nd, icb, plcl, p, th, tv, tvp, thnk, pbase, &
    688     buoybase, iflag, sig, w0)
     755SUBROUTINE cv3_trigger(len, nd, icb, plcl, p, th, tv, tvp, thnk, &
     756                       pbase, buoybase, iflag, sig, w0)
    689757  IMPLICIT NONE
    690758
    691   ! -------------------------------------------------------------------
    692   ! --- TRIGGERING
    693 
    694   ! - computes the cloud base
    695   ! - triggering (crude in this version)
    696   ! - relaxation of sig and w0 when no convection
    697 
    698   ! Caution1: if no convection, we set iflag=4
    699   ! (it used to be 0 in convect3)
    700 
    701   ! Caution2: at this stage, tvp (and thus buoy) are know up
    702   ! through icb only!
    703   ! -> the buoyancy below cloud base not (yet) set to the cloud base buoyancy
    704   ! -------------------------------------------------------------------
     759! -------------------------------------------------------------------
     760! --- TRIGGERING
     761
     762! - computes the cloud base
     763! - triggering (crude in this version)
     764! - relaxation of sig and w0 when no convection
     765
     766! Caution1: if no convection, we set iflag=4
     767! (it used to be 0 in convect3)
     768
     769! Caution2: at this stage, tvp (and thus buoy) are know up
     770! through icb only!
     771! -> the buoyancy below cloud base not (yet) set to the cloud base buoyancy
     772! -------------------------------------------------------------------
    705773
    706774  include "cv3param.h"
    707775
    708   ! input:
     776! input:
    709777  INTEGER len, nd
    710778  INTEGER icb(len)
    … …  
    713781  REAL thnk(len)
    714782
    715   ! output:
     783! output:
    716784  REAL pbase(len), buoybase(len)
    717785
    718   ! input AND output:
     786! input AND output:
    719787  INTEGER iflag(len)
    720788  REAL sig(len, nd), w0(len, nd)
    721789
    722   ! local variables:
     790! local variables:
    723791  INTEGER i, k
    724792  REAL tvpbase, tvbase, tdif, ath, ath1
    725793
    726794
    727   ! ***   set cloud base buoyancy at (plcl+dpbase) level buoyancy
     795! ***   set cloud base buoyancy at (plcl+dpbase) level buoyancy
    728796
    729797  DO i = 1, len
    730798    pbase(i) = plcl(i) + dpbase
    731     tvpbase = tvp(i, icb(i))*(pbase(i)-p(i,icb(i)+1))/ &
    732       (p(i,icb(i))-p(i,icb(i)+1)) + tvp(i, icb(i)+1)*(p(i,icb(i))-pbase(i))/( &
    733       p(i,icb(i))-p(i,icb(i)+1))
    734     tvbase = tv(i, icb(i))*(pbase(i)-p(i,icb(i)+1))/ &
    735       (p(i,icb(i))-p(i,icb(i)+1)) + tv(i, icb(i)+1)*(p(i,icb(i))-pbase(i))/(p &
    736       (i,icb(i))-p(i,icb(i)+1))
     799    tvpbase = tvp(i, icb(i))  *(pbase(i)-p(i,icb(i)+1))/(p(i,icb(i))-p(i,icb(i)+1)) + &
     800              tvp(i, icb(i)+1)*(p(i,icb(i))-pbase(i))  /(p(i,icb(i))-p(i,icb(i)+1))
     801    tvbase = tv(i, icb(i))  *(pbase(i)-p(i,icb(i)+1))/(p(i,icb(i))-p(i,icb(i)+1)) + &
     802             tv(i, icb(i)+1)*(p(i,icb(i))-pbase(i))  /(p(i,icb(i))-p(i,icb(i)+1))
    737803    buoybase(i) = tvpbase - tvbase
    738804  END DO
    739805
    740806
    741   ! ***   make sure that column is dry adiabatic between the surface  ***
    742   ! ***    and cloud base, and that lifted air is positively buoyant  ***
    743   ! ***                         at cloud base                         ***
    744   ! ***       if not, return to calling program after resetting       ***
    745   ! ***                        sig(i) and w0(i)                       ***
    746 
    747 
    748   ! oct3      do 200 i=1,len
    749   ! oct3
    750   ! oct3       tdif = buoybase(i)
    751   ! oct3       ath1 = th(i,1)
    752   ! oct3       ath  = th(i,icb(i)-1) - dttrig
    753   ! oct3
    754   ! oct3       if (tdif.lt.dtcrit .or. ath.gt.ath1) then
    755   ! oct3         do 60 k=1,nl
    756   ! oct3            sig(i,k) = beta*sig(i,k) - 2.*alpha*tdif*tdif
    757   ! oct3            sig(i,k) = AMAX1(sig(i,k),0.0)
    758   ! oct3            w0(i,k)  = beta*w0(i,k)
    759   ! oct3   60    continue
    760   ! oct3         iflag(i)=4 ! pour version vectorisee
    761   ! oct3c convect3         iflag(i)=0
    762   ! oct3cccc         return
    763   ! oct3       endif
    764   ! oct3
    765   ! oct3200   continue
    766 
    767   ! -- oct3: on reecrit la boucle 200 (pour la vectorisation)
     807! ***   make sure that column is dry adiabatic between the surface  ***
     808! ***    and cloud base, and that lifted air is positively buoyant  ***
     809! ***                         at cloud base                         ***
     810! ***       if not, return to calling program after resetting       ***
     811! ***                        sig(i) and w0(i)                       ***
     812
     813
     814! oct3      do 200 i=1,len
     815! oct3
     816! oct3       tdif = buoybase(i)
     817! oct3       ath1 = th(i,1)
     818! oct3       ath  = th(i,icb(i)-1) - dttrig
     819! oct3
     820! oct3       if (tdif.lt.dtcrit .or. ath.gt.ath1) then
     821! oct3         do 60 k=1,nl
     822! oct3            sig(i,k) = beta*sig(i,k) - 2.*alpha*tdif*tdif
     823! oct3            sig(i,k) = AMAX1(sig(i,k),0.0)
     824! oct3            w0(i,k)  = beta*w0(i,k)
     825! oct3   60    continue
     826! oct3         iflag(i)=4 ! pour version vectorisee
     827! oct3c convect3         iflag(i)=0
     828! oct3cccc         return
     829! oct3       endif
     830! oct3
     831! oct3200   continue
     832
     833! -- oct3: on reecrit la boucle 200 (pour la vectorisation)
    768834
    769835  DO k = 1, nl
    … …  
    779845        w0(i, k) = beta*w0(i, k)
    780846        iflag(i) = 4 ! pour version vectorisee
    781         ! convect3         iflag(i)=0
     847! convect3         iflag(i)=0
    782848      END IF
    783849
    … …  
    785851  END DO
    786852
    787   ! fin oct3 --
     853! fin oct3 --
    788854
    789855  RETURN
    790856END SUBROUTINE cv3_trigger
    791857
    792 SUBROUTINE cv3_compress(len, nloc, ncum, nd, ntra, iflag1, nk1, icb1, icbs1, &
    793     plcl1, tnk1, qnk1, gznk1, pbase1, buoybase1, t1, q1, qs1, u1, v1, gz1, &
    794     th1, tra1, h1, lv1, cpn1, p1, ph1, tv1, tp1, tvp1, clw1, sig1, w01, &
    795     iflag, nk, icb, icbs, plcl, tnk, qnk, gznk, pbase, buoybase, t, q, qs, u, &
    796     v, gz, th, tra, h, lv, cpn, p, ph, tv, tp, tvp, clw, sig, w0)
     858SUBROUTINE cv3_compress(len, nloc, ncum, nd, ntra, &
     859                        iflag1, nk1, icb1, icbs1, &
     860                        plcl1, tnk1, qnk1, gznk1, pbase1, buoybase1, &
     861                        t1, q1, qs1, u1, v1, gz1, th1, &
     862                        tra1, &
     863                        h1, lv1, cpn1, p1, ph1, tv1, tp1, tvp1, clw1, &
     864                        sig1, w01, &
     865                        iflag, nk, icb, icbs, &
     866                        plcl, tnk, qnk, gznk, pbase, buoybase, &
     867                        t, q, qs, u, v, gz, th, &
     868                        tra, &
     869                        h, lv, cpn, p, ph, tv, tp, tvp, clw, &
     870                        sig, w0)
    797871  IMPLICIT NONE
    798872
    … …  
    800874  include 'iniprint.h'
    801875
    802   ! inputs:
     876!inputs:
    803877  INTEGER len, ncum, nd, ntra, nloc
    804878  INTEGER iflag1(len), nk1(len), icb1(len), icbs1(len)
    … …  
    813887  REAL tra1(len, nd, ntra)
    814888
    815   ! outputs:
    816   ! en fait, on a nloc=len pour l'instant (cf cv_driver)
     889!outputs:
     890! en fait, on a nloc=len pour l'instant (cf cv_driver)
    817891  INTEGER iflag(nloc), nk(nloc), icb(nloc), icbs(nloc)
    818892  REAL plcl(nloc), tnk(nloc), qnk(nloc), gznk(nloc)
    … …  
    826900  REAL tra(nloc, nd, ntra)
    827901
    828   ! local variables:
     902!local variables:
    829903  INTEGER i, k, nn, j
    830904
    … …  
    859933  END DO
    860934
    861   ! AC!      do 121 j=1,ntra
    862   ! AC!ccccc      do 111 k=1,nl+1
    863   ! AC!      do 111 k=1,nd
    864   ! AC!       nn=0
    865   ! AC!      do 101 i=1,len
    866   ! AC!      if(iflag1(i).eq.0)then
    867   ! AC!       nn=nn+1
    868   ! AC!       tra(nn,k,j)=tra1(i,k,j)
    869   ! AC!      endif
    870   ! AC! 101  continue
    871   ! AC! 111  continue
    872   ! AC! 121  continue
     935!AC!      do 121 j=1,ntra
     936!AC!ccccc      do 111 k=1,nl+1
     937!AC!      do 111 k=1,nd
     938!AC!       nn=0
     939!AC!      do 101 i=1,len
     940!AC!      if(iflag1(i).eq.0)then
     941!AC!       nn=nn+1
     942!AC!       tra(nn,k,j)=tra1(i,k,j)
     943!AC!      endif
     944!AC! 101  continue
     945!AC! 111  continue
     946!AC! 121  continue
    873947
    874948  IF (nn/=ncum) THEN
    … …  
    902976
    903977
    904   ! JAM--------------------------------------------------------------------
    905   ! Calcul de la quantité d'eau sous forme de glace
    906   ! --------------------------------------------------------------------
     978!JAM--------------------------------------------------------------------
     979! Calcul de la quantité d'eau sous forme de glace
     980! --------------------------------------------------------------------
    907981  REAL qi(len, nl)
    908982  REAL t(len, nl), clw(len, nl)
    … …  
    922996        END IF
    923997      END IF
    924       ! print*,t(i,k),qi(i,k),'temp,testglace'
     998! print*,t(i,k),qi(i,k),'temp,testglace'
    925999    END DO
    9261000  END DO
    … …  
    9301004END SUBROUTINE icefrac
    9311005
    932 SUBROUTINE cv3_undilute2(nloc, ncum, nd, icb, icbs, nk, tnk, qnk, gznk, hnk, &
    933     t, q, qs, gz, p, h, tv, lv, lf, pbase, buoybase, plcl, inb, tp, tvp, clw, &
    934     hp, ep, sigp, buoy, frac)
     1006SUBROUTINE cv3_undilute2(nloc, ncum, nd, icb, icbs, nk, &
     1007                         tnk, qnk, gznk, hnk, t, q, qs, gz, &
     1008                         p, h, tv, lv, lf, pbase, buoybase, plcl, &
     1009                         inb, tp, tvp, clw, hp, ep, sigp, buoy, frac)
    9351010  IMPLICIT NONE
    9361011
    937   ! ---------------------------------------------------------------------
    938   ! Purpose:
    939   ! FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
    940   ! &
    941   ! COMPUTE THE PRECIPITATION EFFICIENCIES AND THE
    942   ! FRACTION OF PRECIPITATION FALLING OUTSIDE OF CLOUD
    943   ! &
    944   ! FIND THE LEVEL OF NEUTRAL BUOYANCY
    945 
    946   ! Main differences convect3/convect4:
    947   ! - icbs (input) is the first level above LCL (may differ from icb)
    948   ! - many minor differences in the iterations
    949   ! - condensed water not removed from tvp in convect3
    950   ! - vertical profile of buoyancy computed here (use of buoybase)
    951   ! - the determination of inb is different
    952   ! - no inb1, only inb in output
    953   ! ---------------------------------------------------------------------
     1012! ---------------------------------------------------------------------
     1013! Purpose:
     1014! FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
     1015! &
     1016! COMPUTE THE PRECIPITATION EFFICIENCIES AND THE
     1017! FRACTION OF PRECIPITATION FALLING OUTSIDE OF CLOUD
     1018! &
     1019! FIND THE LEVEL OF NEUTRAL BUOYANCY
     1020
     1021! Main differences convect3/convect4:
     1022!  - icbs (input) is the first level above LCL (may differ from icb)
     1023!  - many minor differences in the iterations
     1024!  - condensed water not removed from tvp in convect3
     1025!  - vertical profile of buoyancy computed here (use of buoybase)
     1026!  - the determination of inb is different
     1027!  - no inb1, only inb in output
     1028! ---------------------------------------------------------------------
    9541029
    9551030  include "cvthermo.h"
    … …  
    9581033  include "cvflag.h"
    9591034
    960   ! inputs:
     1035!inputs:
    9611036  INTEGER ncum, nd, nloc, j
    9621037  INTEGER icb(nloc), icbs(nloc), nk(nloc)
    … …  
    9681043  REAL pbase(nloc), buoybase(nloc), plcl(nloc)
    9691044
    970   ! outputs:
     1045!outputs:
    9711046  INTEGER inb(nloc)
    9721047  REAL tp(nloc, nd), tvp(nloc, nd), clw(nloc, nd)
    … …  
    9741049  REAL buoy(nloc, nd)
    9751050
    976   ! local variables:
     1051!local variables:
    9771052  INTEGER i, k
    9781053  REAL tg, qg, ahg, alv, alf, s, tc, es, esi, denom, rg, tca, elacrit
    … …  
    9861061  REAL fracg
    9871062
    988   ! =====================================================================
    989   ! --- SOME INITIALIZATIONS
    990   ! =====================================================================
     1063! =====================================================================
     1064! --- SOME INITIALIZATIONS
     1065! =====================================================================
    9911066
    9921067  DO k = 1, nl
    … …  
    9981073  END DO
    9991074
    1000   ! =====================================================================
    1001   ! --- FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
    1002   ! =====================================================================
    1003 
    1004   ! ---       The procedure is to solve the equation.
    1005   ! cp*tp+L*qp+phi=cp*tnk+L*qnk+gznk.
    1006 
    1007   ! ***  Calculate certain parcel quantities, including static energy   ***
     1075! =====================================================================
     1076! --- FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
     1077! =====================================================================
     1078
     1079! ---       The procedure is to solve the equation.
     1080!                cp*tp+L*qp+phi=cp*tnk+L*qnk+gznk.
     1081
     1082! ***  Calculate certain parcel quantities, including static energy   ***
    10081083
    10091084
    10101085  DO i = 1, ncum
    1011     ah0(i) = (cpd*(1.-qnk(i))+cl*qnk(i))*tnk(i) & ! debug     &
    1012                                                   ! +qnk(i)*(lv0-clmcpv*(tnk(i)-t0))+gznk(i)
    1013       +qnk(i)*(lv0-clmcpv*(tnk(i)-273.15)) + gznk(i)
    1014   END DO
    1015 
    1016 
    1017   ! ***  Find lifted parcel quantities above cloud base    ***
     1086    ah0(i) = (cpd*(1.-qnk(i))+cl*qnk(i))*tnk(i)+ &
     1087! debug          qnk(i)*(lv0-clmcpv*(tnk(i)-t0))+gznk(i)
     1088             qnk(i)*(lv0-clmcpv*(tnk(i)-273.15)) + gznk(i)
     1089  END DO
     1090
     1091
     1092! ***  Find lifted parcel quantities above cloud base    ***
    10181093
    10191094
    10201095  DO k = minorig + 1, nl
    10211096    DO i = 1, ncum
    1022       ! ori         if(k.ge.(icb(i)+1))then
    1023       IF (k>=(icbs(i)+1)) THEN ! convect3
     1097! ori       if(k.ge.(icb(i)+1))then
     1098      IF (k>=(icbs(i)+1)) THEN                                ! convect3
    10241099        tg = t(i, k)
    10251100        qg = qs(i, k)
    1026         ! debug       alv=lv0-clmcpv*(t(i,k)-t0)
     1101! debug       alv=lv0-clmcpv*(t(i,k)-t0)
    10271102        alv = lv0 - clmcpv*(t(i,k)-273.15)
    10281103
    1029         ! First iteration.
    1030 
    1031         ! ori          s=cpd+alv*alv*qg/(rrv*t(i,k)*t(i,k))
    1032         s = cpd*(1.-qnk(i)) + cl*qnk(i) & ! convect3
    1033           +alv*alv*qg/(rrv*t(i,k)*t(i,k)) ! convect3
     1104! First iteration.
     1105
     1106! ori          s=cpd+alv*alv*qg/(rrv*t(i,k)*t(i,k))
     1107        s = cpd*(1.-qnk(i)) + cl*qnk(i) + &                  ! convect3
     1108            alv*alv*qg/(rrv*t(i,k)*t(i,k))                    ! convect3
    10341109        s = 1./s
    1035         ! ori          ahg=cpd*tg+(cl-cpd)*qnk(i)*t(i,k)+alv*qg+gz(i,k)
     1110! ori          ahg=cpd*tg+(cl-cpd)*qnk(i)*t(i,k)+alv*qg+gz(i,k)
    10361111        ahg = cpd*tg + (cl-cpd)*qnk(i)*tg + alv*qg + gz(i, k) ! convect3
    10371112        tg = tg + s*(ah0(i)-ahg)
    1038         ! ori          tg=max(tg,35.0)
    1039         ! debug        tc=tg-t0
     1113! ori          tg=max(tg,35.0)
     1114! debug        tc=tg-t0
    10401115        tc = tg - 273.15
    10411116        denom = 243.5 + tc
    1042         denom = max(denom, 1.0) ! convect3
    1043         ! ori          if(tc.ge.0.0)then
     1117        denom = max(denom, 1.0)                               ! convect3
     1118! ori          if(tc.ge.0.0)then
    10441119        es = 6.112*exp(17.67*tc/denom)
    1045         ! ori          else
    1046         ! ori                   es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
    1047         ! ori          endif
     1120! ori          else
     1121! ori                   es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
     1122! ori          endif
    10481123        qg = eps*es/(p(i,k)-es*(1.-eps))
    10491124
    1050         ! Second iteration.
    1051 
    1052         ! ori          s=cpd+alv*alv*qg/(rrv*t(i,k)*t(i,k))
    1053         ! ori          s=1./s
    1054         ! ori          ahg=cpd*tg+(cl-cpd)*qnk(i)*t(i,k)+alv*qg+gz(i,k)
     1125! Second iteration.
     1126
     1127! ori          s=cpd+alv*alv*qg/(rrv*t(i,k)*t(i,k))
     1128! ori          s=1./s
     1129! ori          ahg=cpd*tg+(cl-cpd)*qnk(i)*t(i,k)+alv*qg+gz(i,k)
    10551130        ahg = cpd*tg + (cl-cpd)*qnk(i)*tg + alv*qg + gz(i, k) ! convect3
    10561131        tg = tg + s*(ah0(i)-ahg)
    1057         ! ori          tg=max(tg,35.0)
    1058         ! debug        tc=tg-t0
     1132! ori          tg=max(tg,35.0)
     1133! debug        tc=tg-t0
    10591134        tc = tg - 273.15
    10601135        denom = 243.5 + tc
    1061         denom = max(denom, 1.0) ! convect3
    1062         ! ori          if(tc.ge.0.0)then
     1136        denom = max(denom, 1.0)                               ! convect3
     1137! ori          if(tc.ge.0.0)then
    10631138        es = 6.112*exp(17.67*tc/denom)
    1064         ! ori          else
    1065         ! ori                   es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
    1066         ! ori          endif
     1139! ori          else
     1140! ori                   es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
     1141! ori          endif
    10671142        qg = eps*es/(p(i,k)-es*(1.-eps))
    10681143
    1069         ! debug        alv=lv0-clmcpv*(t(i,k)-t0)
     1144! debug        alv=lv0-clmcpv*(t(i,k)-t0)
    10701145        alv = lv0 - clmcpv*(t(i,k)-273.15)
    1071         ! print*,'cpd dans convect2 ',cpd
    1072         ! print*,'tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd'
    1073         ! print*,tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd
    1074 
    1075         ! ori c approximation here:
    1076         ! ori
    1077         ! tp(i,k)=(ah0(i)-(cl-cpd)*qnk(i)*t(i,k)-gz(i,k)-alv*qg)/cpd
    1078 
    1079         ! convect3: no approximation:
     1146! print*,'cpd dans convect2 ',cpd
     1147! print*,'tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd'
     1148! print*,tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd
     1149
     1150! ori c approximation here:
     1151! ori        tp(i,k)=(ah0(i)-(cl-cpd)*qnk(i)*t(i,k)-gz(i,k)-alv*qg)/cpd
     1152
     1153! convect3: no approximation:
    10801154        IF (cvflag_ice) THEN
    10811155          tp(i, k) = max(0., (ah0(i)-gz(i,k)-alv*qg)/(cpd+(cl-cpd)*qnk(i)))
    … …  
    10871161        clw(i, k) = max(0.0, clw(i,k))
    10881162        rg = qg/(1.-qnk(i))
    1089         ! ori               tvp(i,k)=tp(i,k)*(1.+rg*epsi)
    1090         ! convect3: (qg utilise au lieu du vrai mixing ratio rg):
     1163! ori               tvp(i,k)=tp(i,k)*(1.+rg*epsi)
     1164! convect3: (qg utilise au lieu du vrai mixing ratio rg):
    10911165        tvp(i, k) = tp(i, k)*(1.+qg/eps-qnk(i)) ! whole thing
    10921166        IF (cvflag_ice) THEN
    … …  
    10991173
    11001174      IF (cvflag_ice) THEN
    1101         ! CR:attention boucle en klon dans Icefrac
    1102         ! Call Icefrac(t,clw,qi,nl,nloc)
     1175!CR:attention boucle en klon dans Icefrac
     1176! Call Icefrac(t,clw,qi,nl,nloc)
    11031177        IF (t(i,k)>263.15) THEN
    11041178          qi(i, k) = 0.
    … …  
    11111185          END IF
    11121186        END IF
    1113         ! CR: fin test
     1187!CR: fin test
    11141188        IF (t(i,k)<263.15) THEN
    1115           ! CR: on commente les calculs d'Arnaud car division par zero
    1116           ! nouveau calcul propose par JYG
    1117           ! alv=lv0-clmcpv*(t(i,k)-273.15)
    1118           ! alf=lf0-clmci*(t(i,k)-273.15)
    1119           ! tg=tp(i,k)
    1120           ! tc=tp(i,k)-273.15
    1121           ! denom=243.5+tc
    1122           ! do j=1,3
    1123           ! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    1124           ! il faudra que esi vienne en argument de la convection
    1125           ! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    1126           ! tbis=t(i,k)+(tp(i,k)-tg)
    1127           ! esi=exp(23.33086-(6111.72784/tbis)
    1128           ! :               +0.15215*log(tbis))
    1129           ! qsat_new=eps*esi/(p(i,k)-esi*(1.-eps))
    1130           ! snew=cpd*(1.-qnk(i))+cl*qnk(i)+alv*alv*qsat_new/
    1131           ! :                               (rrv*tbis*tbis)
    1132           ! snew=1./snew
    1133           ! print*,esi,qsat_new,snew,'esi,qsat,snew'
    1134           ! tp(i,k)=tg+(alf*qi(i,k)+alv*qg*(1.-(esi/es)))*snew
    1135           ! print*,k,tp(i,k),qnk(i),'avec glace'
    1136           ! print*,'tpNAN',tg,alf,qi(i,k),alv,qg,esi,es,snew
    1137           ! enddo
     1189!CR: on commente les calculs d'Arnaud car division par zero
     1190! nouveau calcul propose par JYG
     1191!      alv=lv0-clmcpv*(t(i,k)-273.15)
     1192!      alf=lf0-clmci*(t(i,k)-273.15)
     1193!      tg=tp(i,k)
     1194!      tc=tp(i,k)-273.15
     1195!      denom=243.5+tc
     1196!      do j=1,3
     1197! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
     1198! il faudra que esi vienne en argument de la convection
     1199! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
     1200!        tbis=t(i,k)+(tp(i,k)-tg)
     1201!        esi=exp(23.33086-(6111.72784/tbis) + &
     1202!                       0.15215*log(tbis))
     1203!        qsat_new=eps*esi/(p(i,k)-esi*(1.-eps))
     1204!        snew=cpd*(1.-qnk(i))+cl*qnk(i)+alv*alv*qsat_new/ &
     1205!                                       (rrv*tbis*tbis)
     1206!        snew=1./snew
     1207!        print*,esi,qsat_new,snew,'esi,qsat,snew'
     1208!        tp(i,k)=tg+(alf*qi(i,k)+alv*qg*(1.-(esi/es)))*snew
     1209!        print*,k,tp(i,k),qnk(i),'avec glace'
     1210!        print*,'tpNAN',tg,alf,qi(i,k),alv,qg,esi,es,snew
     1211!      enddo
    11381212
    11391213          alv = lv0 - clmcpv*(t(i,k)-273.15)
    … …  
    11451219            esi = exp(23.33086-(6111.72784/tp(i,k))+0.15215*log(tp(i,k)))
    11461220            qsat_new = eps*esi/(p(i,k)-esi*(1.-eps))
    1147             snew = cpd*(1.-qnk(i)) + cl*qnk(i) + alv*als*qsat_new/(rrv*tp(i,k &
    1148               )*tp(i,k))
     1221            snew = cpd*(1.-qnk(i)) + cl*qnk(i) + alv*als*qsat_new/ &
     1222                                                 (rrv*tp(i,k)*tp(i,k))
    11491223            snew = 1./snew
    1150             ! c             print*,esi,qsat_new,snew,'esi,qsat,snew'
    1151             tp(i, k) = tp(i, k) + ((cpd*(1.-qnk(i))+cl*qnk(i))*(tg-tp(i, &
    1152               k))+alv*(qg-qsat_new)+alf*qi(i,k))*snew
    1153             ! print*,k,tp(i,k),qsat_new,qnk(i),qi(i,k),
    1154             ! :             'k,tp,q,qt,qi avec glace'
     1224! c             print*,esi,qsat_new,snew,'esi,qsat,snew'
     1225            tp(i, k) = tp(i, k) + &
     1226                       ((cpd*(1.-qnk(i))+cl*qnk(i))*(tg-tp(i,k)) + &
     1227                        alv*(qg-qsat_new)+alf*qi(i,k))*snew
     1228! print*,k,tp(i,k),qsat_new,qnk(i),qi(i,k), &
     1229!              'k,tp,q,qt,qi avec glace'
    11551230          END DO
    11561231
    1157           ! CR:reprise du code AJ
     1232!CR:reprise du code AJ
    11581233          clw(i, k) = qnk(i) - qsat_new
    11591234          clw(i, k) = max(0.0, clw(i,k))
    11601235          tvp(i, k) = max(0., tp(i,k)*(1.+qsat_new/eps-qnk(i)))
    1161           ! print*,tvp(i,k),'tvp'
     1236! print*,tvp(i,k),'tvp'
    11621237        END IF
    11631238        IF (clw(i,k)<1.E-11) THEN
    … …  
    11701245  END DO
    11711246
    1172   ! =====================================================================
    1173   ! --- SET THE PRECIPITATION EFFICIENCIES AND THE FRACTION OF
    1174   ! --- PRECIPITATION FALLING OUTSIDE OF CLOUD
    1175   ! --- THESE MAY BE FUNCTIONS OF TP(I), P(I) AND CLW(I)
    1176   ! =====================================================================
     1247! =====================================================================
     1248! --- SET THE PRECIPITATION EFFICIENCIES AND THE FRACTION OF
     1249! --- PRECIPITATION FALLING OUTSIDE OF CLOUD
     1250! --- THESE MAY BE FUNCTIONS OF TP(I), P(I) AND CLW(I)
     1251! =====================================================================
    11771252
    11781253  IF (flag_epkeorig/=1) THEN
    … …  
    12051280    END DO
    12061281  END IF
    1207   ! =====================================================================
    1208   ! --- CALCULATE VIRTUAL TEMPERATURE AND LIFTED PARCEL
    1209   ! --- VIRTUAL TEMPERATURE
    1210   ! =====================================================================
    1211 
    1212   ! dans convect3, tvp est calcule en une seule fois, et sans retirer
    1213   ! l'eau condensee (~> reversible CAPE)
    1214 
    1215   ! ori      do 340 k=minorig+1,nl
    1216   ! ori        do 330 i=1,ncum
    1217   ! ori        if(k.ge.(icb(i)+1))then
    1218   ! ori          tvp(i,k)=tvp(i,k)*(1.0-qnk(i)+ep(i,k)*clw(i,k))
    1219   ! oric         print*,'i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)'
    1220   ! oric         print*, i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)
    1221   ! ori        endif
    1222   ! ori 330    continue
    1223   ! ori 340  continue
    1224 
    1225   ! ori      do 350 i=1,ncum
    1226   ! ori       tvp(i,nlp)=tvp(i,nl)-(gz(i,nlp)-gz(i,nl))/cpd
    1227   ! ori 350  continue
    1228 
    1229   DO i = 1, ncum ! convect3
    1230     tp(i, nlp) = tp(i, nl) ! convect3
    1231   END DO ! convect3
    1232 
    1233   ! =====================================================================
    1234   ! --- EFFECTIVE VERTICAL PROFILE OF BUOYANCY (convect3 only):
    1235   ! =====================================================================
    1236 
    1237   ! -- this is for convect3 only:
    1238 
    1239   ! first estimate of buoyancy:
     1282! =====================================================================
     1283! --- CALCULATE VIRTUAL TEMPERATURE AND LIFTED PARCEL
     1284! --- VIRTUAL TEMPERATURE
     1285! =====================================================================
     1286
     1287! dans convect3, tvp est calcule en une seule fois, et sans retirer
     1288! l'eau condensee (~> reversible CAPE)
     1289
     1290! ori      do 340 k=minorig+1,nl
     1291! ori        do 330 i=1,ncum
     1292! ori        if(k.ge.(icb(i)+1))then
     1293! ori          tvp(i,k)=tvp(i,k)*(1.0-qnk(i)+ep(i,k)*clw(i,k))
     1294! oric         print*,'i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)'
     1295! oric         print*, i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)
     1296! ori        endif
     1297! ori 330    continue
     1298! ori 340  continue
     1299
     1300! ori      do 350 i=1,ncum
     1301! ori       tvp(i,nlp)=tvp(i,nl)-(gz(i,nlp)-gz(i,nl))/cpd
     1302! ori 350  continue
     1303
     1304  DO i = 1, ncum                                           ! convect3
     1305    tp(i, nlp) = tp(i, nl)                                 ! convect3
     1306  END DO                                                   ! convect3
     1307
     1308! =====================================================================
     1309! --- EFFECTIVE VERTICAL PROFILE OF BUOYANCY (convect3 only):
     1310! =====================================================================
     1311
     1312! -- this is for convect3 only:
     1313
     1314! first estimate of buoyancy:
    12401315
    12411316  DO i = 1, ncum
    … …  
    12451320  END DO
    12461321
    1247   ! set buoyancy=buoybase for all levels below base
    1248   ! for safety, set buoy(icb)=buoybase
     1322! set buoyancy=buoybase for all levels below base
     1323! for safety, set buoy(icb)=buoybase
    12491324
    12501325  DO i = 1, ncum
    … …  
    12541329      END IF
    12551330    END DO
    1256     ! buoy(icb(i),k)=buoybase(i)
     1331!    buoy(icb(i),k)=buoybase(i)
    12571332    buoy(i, icb(i)) = buoybase(i)
    12581333  END DO
    12591334
    1260   ! -- end convect3
    1261 
    1262   ! =====================================================================
    1263   ! --- FIND THE FIRST MODEL LEVEL (INB) ABOVE THE PARCEL'S
    1264   ! --- LEVEL OF NEUTRAL BUOYANCY
    1265   ! =====================================================================
    1266 
    1267   ! -- this is for convect3 only:
     1335! -- end convect3
     1336
     1337! =====================================================================
     1338! --- FIND THE FIRST MODEL LEVEL (INB) ABOVE THE PARCEL'S
     1339! --- LEVEL OF NEUTRAL BUOYANCY
     1340! =====================================================================
     1341
     1342! -- this is for convect3 only:
    12681343
    12691344  DO i = 1, ncum
    … …  
    12731348
    12741349
    1275   ! --    iposit(i) = first level, above icb, with positive buoyancy
     1350! --    iposit(i) = first level, above icb, with positive buoyancy
    12761351  DO k = 1, nl - 1
    12771352    DO i = 1, ncum
    … …  
    12961371  END DO
    12971372
    1298   ! -- end convect3
    1299 
    1300   ! ori      do 510 i=1,ncum
    1301   ! ori        cape(i)=0.0
    1302   ! ori        capem(i)=0.0
    1303   ! ori        inb(i)=icb(i)+1
    1304   ! ori        inb1(i)=inb(i)
    1305   ! ori 510  continue
    1306 
    1307   ! Originial Code
    1308 
    1309   ! do 530 k=minorig+1,nl-1
    1310   ! do 520 i=1,ncum
    1311   ! if(k.ge.(icb(i)+1))then
    1312   ! by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
    1313   ! byp=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
    1314   ! cape(i)=cape(i)+by
    1315   ! if(by.ge.0.0)inb1(i)=k+1
    1316   ! if(cape(i).gt.0.0)then
    1317   ! inb(i)=k+1
    1318   ! capem(i)=cape(i)
    1319   ! endif
    1320   ! endif
    1321   ! 520    continue
    1322   ! 530  continue
    1323   ! do 540 i=1,ncum
    1324   ! byp=(tvp(i,nl)-tv(i,nl))*dph(i,nl)/p(i,nl)
    1325   ! cape(i)=capem(i)+byp
    1326   ! defrac=capem(i)-cape(i)
    1327   ! defrac=max(defrac,0.001)
    1328   ! frac(i)=-cape(i)/defrac
    1329   ! frac(i)=min(frac(i),1.0)
    1330   ! frac(i)=max(frac(i),0.0)
    1331   ! 540   continue
    1332 
    1333   ! K Emanuel fix
    1334 
    1335   ! call zilch(byp,ncum)
    1336   ! do 530 k=minorig+1,nl-1
    1337   ! do 520 i=1,ncum
    1338   ! if(k.ge.(icb(i)+1))then
    1339   ! by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
    1340   ! cape(i)=cape(i)+by
    1341   ! if(by.ge.0.0)inb1(i)=k+1
    1342   ! if(cape(i).gt.0.0)then
    1343   ! inb(i)=k+1
    1344   ! capem(i)=cape(i)
    1345   ! byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
    1346   ! endif
    1347   ! endif
    1348   ! 520    continue
    1349   ! 530  continue
    1350   ! do 540 i=1,ncum
    1351   ! inb(i)=max(inb(i),inb1(i))
    1352   ! cape(i)=capem(i)+byp(i)
    1353   ! defrac=capem(i)-cape(i)
    1354   ! defrac=max(defrac,0.001)
    1355   ! frac(i)=-cape(i)/defrac
    1356   ! frac(i)=min(frac(i),1.0)
    1357   ! frac(i)=max(frac(i),0.0)
    1358   ! 540   continue
    1359 
    1360   ! J Teixeira fix
    1361 
    1362   ! ori      call zilch(byp,ncum)
    1363   ! ori      do 515 i=1,ncum
    1364   ! ori        lcape(i)=.true.
    1365   ! ori 515  continue
    1366   ! ori      do 530 k=minorig+1,nl-1
    1367   ! ori        do 520 i=1,ncum
    1368   ! ori          if(cape(i).lt.0.0)lcape(i)=.false.
    1369   ! ori          if((k.ge.(icb(i)+1)).and.lcape(i))then
    1370   ! ori            by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
    1371   ! ori            byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
    1372   ! ori            cape(i)=cape(i)+by
    1373   ! ori            if(by.ge.0.0)inb1(i)=k+1
    1374   ! ori            if(cape(i).gt.0.0)then
    1375   ! ori              inb(i)=k+1
    1376   ! ori              capem(i)=cape(i)
    1377   ! ori            endif
    1378   ! ori          endif
    1379   ! ori 520    continue
    1380   ! ori 530  continue
    1381   ! ori      do 540 i=1,ncum
    1382   ! ori          cape(i)=capem(i)+byp(i)
    1383   ! ori          defrac=capem(i)-cape(i)
    1384   ! ori          defrac=max(defrac,0.001)
    1385   ! ori          frac(i)=-cape(i)/defrac
    1386   ! ori          frac(i)=min(frac(i),1.0)
    1387   ! ori          frac(i)=max(frac(i),0.0)
    1388   ! ori 540  continue
    1389 
    1390   ! =====================================================================
    1391   ! ---   CALCULATE LIQUID WATER STATIC ENERGY OF LIFTED PARCEL
    1392   ! =====================================================================
     1373! -- end convect3
     1374
     1375! ori      do 510 i=1,ncum
     1376! ori        cape(i)=0.0
     1377! ori        capem(i)=0.0
     1378! ori        inb(i)=icb(i)+1
     1379! ori        inb1(i)=inb(i)
     1380! ori 510  continue
     1381
     1382! Originial Code
     1383
     1384!    do 530 k=minorig+1,nl-1
     1385!    do 520 i=1,ncum
     1386!      if(k.ge.(icb(i)+1))then
     1387!      by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
     1388!      byp=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
     1389!      cape(i)=cape(i)+by
     1390!      if(by.ge.0.0)inb1(i)=k+1
     1391!      if(cape(i).gt.0.0)then
     1392!        inb(i)=k+1
     1393!        capem(i)=cape(i)
     1394!      endif
     1395!      endif
     1396!520    continue
     1397!530  continue
     1398!    do 540 i=1,ncum
     1399!    byp=(tvp(i,nl)-tv(i,nl))*dph(i,nl)/p(i,nl)
     1400!    cape(i)=capem(i)+byp
     1401!    defrac=capem(i)-cape(i)
     1402!    defrac=max(defrac,0.001)
     1403!    frac(i)=-cape(i)/defrac
     1404!    frac(i)=min(frac(i),1.0)
     1405!    frac(i)=max(frac(i),0.0)
     1406!540   continue
     1407
     1408!    K Emanuel fix
     1409
     1410!    call zilch(byp,ncum)
     1411!    do 530 k=minorig+1,nl-1
     1412!    do 520 i=1,ncum
     1413!      if(k.ge.(icb(i)+1))then
     1414!      by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
     1415!      cape(i)=cape(i)+by
     1416!      if(by.ge.0.0)inb1(i)=k+1
     1417!      if(cape(i).gt.0.0)then
     1418!        inb(i)=k+1
     1419!        capem(i)=cape(i)
     1420!        byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
     1421!      endif
     1422!      endif
     1423!520    continue
     1424!530  continue
     1425!    do 540 i=1,ncum
     1426!    inb(i)=max(inb(i),inb1(i))
     1427!    cape(i)=capem(i)+byp(i)
     1428!    defrac=capem(i)-cape(i)
     1429!    defrac=max(defrac,0.001)
     1430!    frac(i)=-cape(i)/defrac
     1431!    frac(i)=min(frac(i),1.0)
     1432!    frac(i)=max(frac(i),0.0)
     1433!540   continue
     1434
     1435! J Teixeira fix
     1436
     1437! ori      call zilch(byp,ncum)
     1438! ori      do 515 i=1,ncum
     1439! ori        lcape(i)=.true.
     1440! ori 515  continue
     1441! ori      do 530 k=minorig+1,nl-1
     1442! ori        do 520 i=1,ncum
     1443! ori          if(cape(i).lt.0.0)lcape(i)=.false.
     1444! ori          if((k.ge.(icb(i)+1)).and.lcape(i))then
     1445! ori            by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
     1446! ori            byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
     1447! ori            cape(i)=cape(i)+by
     1448! ori            if(by.ge.0.0)inb1(i)=k+1
     1449! ori            if(cape(i).gt.0.0)then
     1450! ori              inb(i)=k+1
     1451! ori              capem(i)=cape(i)
     1452! ori            endif
     1453! ori          endif
     1454! ori 520    continue
     1455! ori 530  continue
     1456! ori      do 540 i=1,ncum
     1457! ori          cape(i)=capem(i)+byp(i)
     1458! ori          defrac=capem(i)-cape(i)
     1459! ori          defrac=max(defrac,0.001)
     1460! ori          frac(i)=-cape(i)/defrac
     1461! ori          frac(i)=min(frac(i),1.0)
     1462! ori          frac(i)=max(frac(i),0.0)
     1463! ori 540  continue
     1464
     1465! =====================================================================
     1466! ---   CALCULATE LIQUID WATER STATIC ENERGY OF LIFTED PARCEL
     1467! =====================================================================
    13931468
    13941469  DO k = 1, nd
    … …  
    14051480          frac(i, k) = 1. - (t(i,k)-243.15)/(263.15-243.15)
    14061481          frac(i, k) = min(max(frac(i,k),0.0), 1.0)
    1407           hp(i, k) = hnk(i) + (lv(i,k)+(cpd-cpv)*t(i,k)+frac(i,k)*lf(i,k))*ep &
    1408             (i, k)*clw(i, k)
     1482          hp(i, k) = hnk(i) + (lv(i,k)+(cpd-cpv)*t(i,k)+frac(i,k)*lf(i,k))* &
     1483                              ep(i, k)*clw(i, k)
    14091484
    14101485        ELSE
    … …  
    14191494END SUBROUTINE cv3_undilute2
    14201495
    1421 SUBROUTINE cv3_closure(nloc, ncum, nd, icb, inb, pbase, p, ph, tv, buoy, sig, &
    1422     w0, cape, m, iflag)
     1496SUBROUTINE cv3_closure(nloc, ncum, nd, icb, inb, &
     1497                       pbase, p, ph, tv, buoy, &
     1498                       sig, w0, cape, m, iflag)
    14231499  IMPLICIT NONE
    14241500
    1425   ! ===================================================================
    1426   ! ---  CLOSURE OF CONVECT3
    1427 
    1428   ! vectorization: S. Bony
    1429   ! ===================================================================
     1501! ===================================================================
     1502! ---  CLOSURE OF CONVECT3
     1503!
     1504! vectorization: S. Bony
     1505! ===================================================================
    14301506
    14311507  include "cvthermo.h"
    14321508  include "cv3param.h"
    14331509
    1434   ! input:
     1510!input:
    14351511  INTEGER ncum, nd, nloc
    14361512  INTEGER icb(nloc), inb(nloc)
    … …  
    14391515  REAL tv(nloc, nd), buoy(nloc, nd)
    14401516
    1441   ! input/output:
     1517!input/output:
    14421518  REAL sig(nloc, nd), w0(nloc, nd)
    14431519  INTEGER iflag(nloc)
    14441520
    1445   ! output:
     1521!output:
    14461522  REAL cape(nloc)
    14471523  REAL m(nloc, nd)
    14481524
    1449   ! local variables:
     1525!local variables:
    14501526  INTEGER i, j, k, icbmax
    14511527  REAL deltap, fac, w, amu
    … …  
    14541530
    14551531
    1456   ! -------------------------------------------------------
    1457   ! -- Initialization
    1458   ! -------------------------------------------------------
     1532! -------------------------------------------------------
     1533! -- Initialization
     1534! -------------------------------------------------------
    14591535
    14601536  DO k = 1, nl
    … …  
    14641540  END DO
    14651541
    1466   ! -------------------------------------------------------
    1467   ! -- Reset sig(i) and w0(i) for i>inb and i<icb
    1468   ! -------------------------------------------------------
    1469 
    1470   ! update sig and w0 above LNB:
     1542! -------------------------------------------------------
     1543! -- Reset sig(i) and w0(i) for i>inb and i<icb
     1544! -------------------------------------------------------
     1545
     1546! update sig and w0 above LNB:
    14711547
    14721548  DO k = 1, nl - 1
    14731549    DO i = 1, ncum
    14741550      IF ((inb(i)<(nl-1)) .AND. (k>=(inb(i)+1))) THEN
    1475         sig(i, k) = beta*sig(i, k) + 2.*alpha*buoy(i, inb(i))*abs(buoy(i,inb( &
    1476           i)))
     1551        sig(i, k) = beta*sig(i, k) + &
     1552                    2.*alpha*buoy(i, inb(i))*abs(buoy(i,inb(i)))
    14771553        sig(i, k) = amax1(sig(i,k), 0.0)
    14781554        w0(i, k) = beta*w0(i, k)
    … …  
    14811557  END DO
    14821558
    1483   ! compute icbmax:
     1559! compute icbmax:
    14841560
    14851561  icbmax = 2
    … …  
    14881564  END DO
    14891565
    1490   ! update sig and w0 below cloud base:
     1566! update sig and w0 below cloud base:
    14911567
    14921568  DO k = 1, icbmax
    14931569    DO i = 1, ncum
    14941570      IF (k<=icb(i)) THEN
    1495         sig(i, k) = beta*sig(i, k) - 2.*alpha*buoy(i, icb(i))*buoy(i, icb(i))
     1571        sig(i, k) = beta*sig(i, k) - &
     1572                    2.*alpha*buoy(i, icb(i))*buoy(i, icb(i))
    14961573        sig(i, k) = max(sig(i,k), 0.0)
    14971574        w0(i, k) = beta*w0(i, k)
    … …  
    15001577  END DO
    15011578
    1502   ! !      if(inb.lt.(nl-1))then
    1503   ! !         do 85 i=inb+1,nl-1
    1504   ! !            sig(i)=beta*sig(i)+2.*alpha*buoy(inb)*
    1505   ! !     1              abs(buoy(inb))
    1506   ! !            sig(i)=max(sig(i),0.0)
    1507   ! !            w0(i)=beta*w0(i)
    1508   ! !   85    continue
    1509   ! !      end if
    1510 
    1511   ! !      do 87 i=1,icb
    1512   ! !         sig(i)=beta*sig(i)-2.*alpha*buoy(icb)*buoy(icb)
    1513   ! !         sig(i)=max(sig(i),0.0)
    1514   ! !         w0(i)=beta*w0(i)
    1515   ! !   87 continue
    1516 
    1517   ! -------------------------------------------------------------
    1518   ! -- Reset fractional areas of updrafts and w0 at initial time
    1519   ! -- and after 10 time steps of no convection
    1520   ! -------------------------------------------------------------
     1579!!      if(inb.lt.(nl-1))then
     1580!!         do 85 i=inb+1,nl-1
     1581!!            sig(i)=beta*sig(i)+2.*alpha*buoy(inb)*
     1582!!     1              abs(buoy(inb))
     1583!!            sig(i)=max(sig(i),0.0)
     1584!!            w0(i)=beta*w0(i)
     1585!!   85    continue
     1586!!      end if
     1587
     1588!!      do 87 i=1,icb
     1589!!         sig(i)=beta*sig(i)-2.*alpha*buoy(icb)*buoy(icb)
     1590!!         sig(i)=max(sig(i),0.0)
     1591!!         w0(i)=beta*w0(i)
     1592!!   87 continue
     1593
     1594! -------------------------------------------------------------
     1595! -- Reset fractional areas of updrafts and w0 at initial time
     1596! -- and after 10 time steps of no convection
     1597! -------------------------------------------------------------
    15211598
    15221599  DO k = 1, nl - 1
    … …  
    15291606  END DO
    15301607
    1531   ! -------------------------------------------------------------
    1532   ! -- Calculate convective available potential energy (cape),
    1533   ! -- vertical velocity (w), fractional area covered by
    1534   ! -- undilute updraft (sig), and updraft mass flux (m)
    1535   ! -------------------------------------------------------------
     1608! -------------------------------------------------------------
     1609! -- Calculate convective available potential energy (cape),
     1610! -- vertical velocity (w), fractional area covered by
     1611! -- undilute updraft (sig), and updraft mass flux (m)
     1612! -------------------------------------------------------------
    15361613
    15371614  DO i = 1, ncum
    … …  
    15391616  END DO
    15401617
    1541   ! compute dtmin (minimum buoyancy between ICB and given level k):
     1618! compute dtmin (minimum buoyancy between ICB and given level k):
    15421619
    15431620  DO i = 1, ncum
    … …  
    15501627    DO k = 1, nl
    15511628      DO j = minorig, nl
    1552         IF ((k>=(icb(i)+1)) .AND. (k<=inb(i)) .AND. (j>=icb(i)) .AND. (j<=(k- &
    1553             1))) THEN
     1629        IF ((k>=(icb(i)+1)) .AND. (k<=inb(i)) .AND. (j>=icb(i)) .AND. (j<=(k-1))) THEN
    15541630          dtmin(i, k) = amin1(dtmin(i,k), buoy(i,j))
    15551631        END IF
    … …  
    15581634  END DO
    15591635
    1560   ! the interval on which cape is computed starts at pbase :
     1636! the interval on which cape is computed starts at pbase :
    15611637
    15621638  DO k = 1, nl
    … …  
    15701646        sigold(i, k) = sig(i, k)
    15711647
    1572         ! dtmin(i,k)=100.0
    1573         ! do 97 j=icb(i),k-1 ! mauvaise vectorisation
    1574         ! dtmin(i,k)=AMIN1(dtmin(i,k),buoy(i,j))
    1575         ! 97     continue
     1648! dtmin(i,k)=100.0
     1649! do 97 j=icb(i),k-1 ! mauvaise vectorisation
     1650! dtmin(i,k)=AMIN1(dtmin(i,k),buoy(i,j))
     1651! 97     continue
    15761652
    15771653        sig(i, k) = beta*sig(i, k) + alpha*dtmin(i, k)*abs(dtmin(i,k))
    … …  
    15901666  DO i = 1, ncum
    15911667    w0(i, icb(i)) = 0.5*w0(i, icb(i)+1)
    1592     m(i, icb(i)) = 0.5*m(i, icb(i)+1)*(ph(i,icb(i))-ph(i,icb(i)+1))/ &
    1593       (ph(i,icb(i)+1)-ph(i,icb(i)+2))
     1668    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))
    15941669    sig(i, icb(i)) = sig(i, icb(i)+1)
    15951670    sig(i, icb(i)-1) = sig(i, icb(i))
    15961671  END DO
    15971672
    1598   ! ccc 3. Compute final cloud base mass flux and set iflag to 3 if
    1599   ! ccc    cloud base mass flux is exceedingly small and is decreasing (i.e.
    1600   ! if
    1601   ! ccc    the final mass flux (cbmflast) is greater than the target mass
    1602   ! flux
    1603   ! ccc    (cbmf) ??).
    1604   ! cc
    1605   ! c      do i = 1,ncum
    1606   ! c       cbmflast(i) = 0.
    1607   ! c      enddo
    1608   ! cc
    1609   ! c      do k= 1,nl
    1610   ! c       do i = 1,ncum
    1611   ! c        IF (k .ge. icb(i) .and. k .le. inb(i)) THEN
    1612   ! c         cbmflast(i) = cbmflast(i)+M(i,k)
    1613   ! c        ENDIF
    1614   ! c       enddo
    1615   ! c      enddo
    1616   ! cc
    1617   ! c      do i = 1,ncum
    1618   ! c       IF (cbmflast(i) .lt. 1.e-6) THEN
    1619   ! c         iflag(i) = 3
    1620   ! c       ENDIF
    1621   ! c      enddo
    1622   ! cc
    1623   ! c      do k= 1,nl
    1624   ! c       do i = 1,ncum
    1625   ! c        IF (iflag(i) .ge. 3) THEN
    1626   ! c         M(i,k) = 0.
    1627   ! c         sig(i,k) = 0.
    1628   ! c         w0(i,k) = 0.
    1629   ! c        ENDIF
    1630   ! c       enddo
    1631   ! c      enddo
    1632   ! cc
    1633   ! !      cape=0.0
    1634   ! !      do 98 i=icb+1,inb
    1635   ! !         deltap = min(pbase,ph(i-1))-min(pbase,ph(i))
    1636   ! !         cape=cape+rrd*buoy(i-1)*deltap/p(i-1)
    1637   ! !         dcape=rrd*buoy(i-1)*deltap/p(i-1)
    1638   ! !         dlnp=deltap/p(i-1)
    1639   ! !         cape=max(0.0,cape)
    1640   ! !         sigold=sig(i)
    1641 
    1642   ! !         dtmin=100.0
    1643   ! !         do 97 j=icb,i-1
    1644   ! !            dtmin=amin1(dtmin,buoy(j))
    1645   ! !   97    continue
    1646 
    1647   ! !         sig(i)=beta*sig(i)+alpha*dtmin*abs(dtmin)
    1648   ! !         sig(i)=max(sig(i),0.0)
    1649   ! !         sig(i)=amin1(sig(i),0.01)
    1650   ! !         fac=amin1(((dtcrit-dtmin)/dtcrit),1.0)
    1651   ! !         w=(1.-beta)*fac*sqrt(cape)+beta*w0(i)
    1652   ! !         amu=0.5*(sig(i)+sigold)*w
    1653   ! !         m(i)=amu*0.007*p(i)*(ph(i)-ph(i+1))/tv(i)
    1654   ! !         w0(i)=w
    1655   ! !   98 continue
    1656   ! !      w0(icb)=0.5*w0(icb+1)
    1657   ! !      m(icb)=0.5*m(icb+1)*(ph(icb)-ph(icb+1))/(ph(icb+1)-ph(icb+2))
    1658   ! !      sig(icb)=sig(icb+1)
    1659   ! !      sig(icb-1)=sig(icb)
     1673! ccc 3. Compute final cloud base mass flux and set iflag to 3 if
     1674! ccc    cloud base mass flux is exceedingly small and is decreasing (i.e. if
     1675! ccc    the final mass flux (cbmflast) is greater than the target mass flux
     1676! ccc    (cbmf) ??).
     1677! cc
     1678! c      do i = 1,ncum
     1679! c       cbmflast(i) = 0.
     1680! c      enddo
     1681! cc
     1682! c      do k= 1,nl
     1683! c       do i = 1,ncum
     1684! c        IF (k .ge. icb(i) .and. k .le. inb(i)) THEN
     1685! c         cbmflast(i) = cbmflast(i)+M(i,k)
     1686! c        ENDIF
     1687! c       enddo
     1688! c      enddo
     1689! cc
     1690! c      do i = 1,ncum
     1691! c       IF (cbmflast(i) .lt. 1.e-6) THEN
     1692! c         iflag(i) = 3
     1693! c       ENDIF
     1694! c      enddo
     1695! cc
     1696! c      do k= 1,nl
     1697! c       do i = 1,ncum
     1698! c        IF (iflag(i) .ge. 3) THEN
     1699! c         M(i,k) = 0.
     1700! c         sig(i,k) = 0.
     1701! c         w0(i,k) = 0.
     1702! c        ENDIF
     1703! c       enddo
     1704! c      enddo
     1705! cc
     1706!!      cape=0.0
     1707!!      do 98 i=icb+1,inb
     1708!!         deltap = min(pbase,ph(i-1))-min(pbase,ph(i))
     1709!!         cape=cape+rrd*buoy(i-1)*deltap/p(i-1)
     1710!!         dcape=rrd*buoy(i-1)*deltap/p(i-1)
     1711!!         dlnp=deltap/p(i-1)
     1712!!         cape=max(0.0,cape)
     1713!!         sigold=sig(i)
     1714
     1715!!         dtmin=100.0
     1716!!         do 97 j=icb,i-1
     1717!!            dtmin=amin1(dtmin,buoy(j))
     1718!!   97    continue
     1719
     1720!!         sig(i)=beta*sig(i)+alpha*dtmin*abs(dtmin)
     1721!!         sig(i)=max(sig(i),0.0)
     1722!!         sig(i)=amin1(sig(i),0.01)
     1723!!         fac=amin1(((dtcrit-dtmin)/dtcrit),1.0)
     1724!!         w=(1.-beta)*fac*sqrt(cape)+beta*w0(i)
     1725!!         amu=0.5*(sig(i)+sigold)*w
     1726!!         m(i)=amu*0.007*p(i)*(ph(i)-ph(i+1))/tv(i)
     1727!!         w0(i)=w
     1728!!   98 continue
     1729!!      w0(icb)=0.5*w0(icb+1)
     1730!!      m(icb)=0.5*m(icb+1)*(ph(icb)-ph(icb+1))/(ph(icb+1)-ph(icb+2))
     1731!!      sig(icb)=sig(icb+1)
     1732!!      sig(icb-1)=sig(icb)
    16601733
    16611734  RETURN
    16621735END SUBROUTINE cv3_closure
    16631736
    1664 SUBROUTINE cv3_mixing(nloc, ncum, nd, na, ntra, icb, nk, inb, ph, t, rr, rs, &
    1665     u, v, tra, h, lv, lf, frac, qnk, unk, vnk, hp, tv, tvp, ep, clw, m, sig, &
    1666     ment, qent, uent, vent, nent, sij, elij, ments, qents, traent)
     1737SUBROUTINE cv3_mixing(nloc, ncum, nd, na, ntra, icb, nk, inb, &
     1738                      ph, t, rr, rs, u, v, tra, h, lv, lf, frac, qnk, &
     1739                      unk, vnk, hp, tv, tvp, ep, clw, m, sig, &
     1740                      ment, qent, uent, vent, nent, sij, elij, ments, qents, traent)
    16671741  IMPLICIT NONE
    16681742
    1669   ! ---------------------------------------------------------------------
    1670   ! a faire:
    1671   ! - vectorisation de la partie normalisation des flux (do 789...)
    1672   ! ---------------------------------------------------------------------
     1743! ---------------------------------------------------------------------
     1744! a faire:
     1745! - vectorisation de la partie normalisation des flux (do 789...)
     1746! ---------------------------------------------------------------------
    16731747
    16741748  include "cvthermo.h"
    … …  
    16761750  include "cvflag.h"
    16771751
    1678   ! inputs:
     1752!inputs:
    16791753  INTEGER ncum, nd, na, ntra, nloc
    16801754  INTEGER icb(nloc), inb(nloc), nk(nloc)
    … …  
    16901764  REAL m(nloc, na) ! input of convect3
    16911765
    1692   ! outputs:
     1766!outputs:
    16931767  REAL ment(nloc, na, na), qent(nloc, na, na)
    16941768  REAL uent(nloc, na, na), vent(nloc, na, na)
    … …  
    16991773  INTEGER nent(nloc, nd)
    17001774
    1701   ! local variables:
     1775!local variables:
    17021776  INTEGER i, j, k, il, im, jm
    17031777  INTEGER num1, num2
    … …  
    17101784  LOGICAL lwork(nloc)
    17111785
    1712   ! =====================================================================
    1713   ! --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS
    1714   ! =====================================================================
    1715 
    1716   ! ori        do 360 i=1,ncum*nlp
     1786! =====================================================================
     1787! --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS
     1788! =====================================================================
     1789
     1790! ori        do 360 i=1,ncum*nlp
    17171791  DO j = 1, nl
    17181792    DO i = 1, ncum
    17191793      nent(i, j) = 0
    1720       ! in convect3, m is computed in cv3_closure
    1721       ! ori          m(i,1)=0.0
    1722     END DO
    1723   END DO
    1724 
    1725   ! ori      do 400 k=1,nlp
    1726   ! ori       do 390 j=1,nlp
     1794! in convect3, m is computed in cv3_closure
     1795! ori          m(i,1)=0.0
     1796    END DO
     1797  END DO
     1798
     1799! ori      do 400 k=1,nlp
     1800! ori       do 390 j=1,nlp
    17271801  DO j = 1, nl
    17281802    DO k = 1, nl
    … …  
    17321806        vent(i, k, j) = v(i, j)
    17331807        elij(i, k, j) = 0.0
    1734         ! ym            ment(i,k,j)=0.0
    1735         ! ym            sij(i,k,j)=0.0
     1808!ym            ment(i,k,j)=0.0
     1809!ym            sij(i,k,j)=0.0
    17361810      END DO
    17371811    END DO
    17381812  END DO
    17391813
    1740   ! ym
     1814!ym
    17411815  ment(1:ncum, 1:nd, 1:nd) = 0.0
    17421816  sij(1:ncum, 1:nd, 1:nd) = 0.0
    17431817
    1744   ! AC!      do k=1,ntra
    1745   ! AC!       do j=1,nd  ! instead nlp
    1746   ! AC!        do i=1,nd ! instead nlp
    1747   ! AC!         do il=1,ncum
    1748   ! AC!            traent(il,i,j,k)=tra(il,j,k)
    1749   ! AC!         enddo
    1750   ! AC!        enddo
    1751   ! AC!       enddo
    1752   ! AC!      enddo
     1818!AC!      do k=1,ntra
     1819!AC!       do j=1,nd  ! instead nlp
     1820!AC!        do i=1,nd ! instead nlp
     1821!AC!         do il=1,ncum
     1822!AC!            traent(il,i,j,k)=tra(il,j,k)
     1823!AC!         enddo
     1824!AC!        enddo
     1825!AC!       enddo
     1826!AC!      enddo
    17531827  zm(:, :) = 0.
    17541828
    1755   ! =====================================================================
    1756   ! --- CALCULATE ENTRAINED AIR MASS FLUX (ment), TOTAL WATER MIXING
    1757   ! --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
    1758   ! --- FRACTION (sij)
    1759   ! =====================================================================
     1829! =====================================================================
     1830! --- CALCULATE ENTRAINED AIR MASS FLUX (ment), TOTAL WATER MIXING
     1831! --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
     1832! --- FRACTION (sij)
     1833! =====================================================================
    17601834
    17611835  DO i = minorig + 1, nl
    … …  
    17631837    DO j = minorig, nl
    17641838      DO il = 1, ncum
    1765         IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (j>=(icb(il)- &
    1766             1)) .AND. (j<=inb(il))) THEN
     1839        IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (j>=(icb(il)-1)) .AND. (j<=inb(il))) THEN
    17671840
    17681841          rti = qnk(il) - ep(il, i)*clw(il, i)
    … …  
    17711844
    17721845          IF (cvflag_ice) THEN
    1773             ! print*,cvflag_ice,'cvflag_ice dans do 700'
     1846! print*,cvflag_ice,'cvflag_ice dans do 700'
    17741847            IF (t(il,j)<=263.15) THEN
    1775               bf2 = 1. + (lf(il,j)+lv(il,j))*(lv(il,j)+frac(il,j)*lf(il,j))* &
    1776                 rs(il, j)/(rrv*t(il,j)*t(il,j)*cpd)
     1848              bf2 = 1. + (lf(il,j)+lv(il,j))*(lv(il,j)+frac(il,j)* &
     1849                   lf(il,j))*rs(il, j)/(rrv*t(il,j)*t(il,j)*cpd)
    17771850            END IF
    17781851          END IF
    … …  
    17911864
    17921865            IF (cvflag_ice) THEN
    1793               anum = anum - (lv(il,j)+frac(il,j)*lf(il,j))*(rti-rs(il,j)-cwat &
    1794                 *bf2)
     1866              anum = anum - (lv(il,j)+frac(il,j)*lf(il,j))*(rti-rs(il,j)-cwat*bf2)
    17951867              denom = denom + (lv(il,j)+frac(il,j)*lf(il,j))*(rr(il,i)-rti)
    17961868            ELSE
    … …  
    18011873            IF (abs(denom)<0.01) denom = 0.01
    18021874            sij(il, i, j) = anum/denom
    1803             altem = sij(il, i, j)*rr(il, i) + (1.-sij(il,i,j))*rti - &
    1804               rs(il, j)
     1875            altem = sij(il, i, j)*rr(il, i) + (1.-sij(il,i,j))*rti - rs(il, j)
    18051876            altem = altem - (bf2-1.)*cwat
    18061877          END IF
    18071878          IF (sij(il,i,j)>0.0 .AND. sij(il,i,j)<0.95) THEN
    18081879            qent(il, i, j) = sij(il, i, j)*rr(il, i) + (1.-sij(il,i,j))*rti
    1809             uent(il, i, j) = sij(il, i, j)*u(il, i) + &
    1810               (1.-sij(il,i,j))*unk(il)
    1811             vent(il, i, j) = sij(il, i, j)*v(il, i) + &
    1812               (1.-sij(il,i,j))*vnk(il)
    1813             ! !!!      do k=1,ntra
    1814             ! !!!      traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k)
    1815             ! !!!     :      +(1.-sij(il,i,j))*tra(il,nk(il),k)
    1816             ! !!!      end do
     1880            uent(il, i, j) = sij(il, i, j)*u(il, i) + (1.-sij(il,i,j))*unk(il)
     1881            vent(il, i, j) = sij(il, i, j)*v(il, i) + (1.-sij(il,i,j))*vnk(il)
     1882!!!!      do k=1,ntra
     1883!!!!      traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k)
     1884!!!!     :      +(1.-sij(il,i,j))*tra(il,nk(il),k)
     1885!!!!      end do
    18171886            elij(il, i, j) = altem
    18181887            elij(il, i, j) = max(0.0, elij(il,i,j))
    … …  
    18261895    END DO
    18271896
    1828     ! AC!       do k=1,ntra
    1829     ! AC!        do j=minorig,nl
    1830     ! AC!         do il=1,ncum
    1831     ! AC!          if( (i.ge.icb(il)).and.(i.le.inb(il)).and.
    1832     ! AC!     :       (j.ge.(icb(il)-1)).and.(j.le.inb(il)))then
    1833     ! AC!            traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k)
    1834     ! AC!     :            +(1.-sij(il,i,j))*tra(il,nk(il),k)
    1835     ! AC!          endif
    1836     ! AC!         enddo
    1837     ! AC!        enddo
    1838     ! AC!       enddo
    1839 
    1840 
    1841     ! ***   if no air can entrain at level i assume that updraft detrains
    1842     ! ***
    1843     ! ***   at that level and calculate detrained air flux and properties
    1844     ! ***
    1845 
    1846 
    1847     ! @      do 170 i=icb(il),inb(il)
     1897!AC!       do k=1,ntra
     1898!AC!        do j=minorig,nl
     1899!AC!         do il=1,ncum
     1900!AC!          if( (i.ge.icb(il)).and.(i.le.inb(il)).and.
     1901!AC!     :       (j.ge.(icb(il)-1)).and.(j.le.inb(il)))then
     1902!AC!            traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k)
     1903!AC!     :            +(1.-sij(il,i,j))*tra(il,nk(il),k)
     1904!AC!          endif
     1905!AC!         enddo
     1906!AC!        enddo
     1907!AC!       enddo
     1908
     1909
     1910! ***   if no air can entrain at level i assume that updraft detrains  ***
     1911! ***   at that level and calculate detrained air flux and properties  ***
     1912
     1913
     1914! @      do 170 i=icb(il),inb(il)
    18481915
    18491916    DO il = 1, ncum
    18501917      IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (nent(il,i)==0)) THEN
    1851         ! @      if(nent(il,i).eq.0)then
     1918! @      if(nent(il,i).eq.0)then
    18521919        ment(il, i, i) = m(il, i)
    18531920        qent(il, i, i) = qnk(il) - ep(il, i)*clw(il, i)
    … …  
    18551922        vent(il, i, i) = vnk(il)
    18561923        elij(il, i, i) = clw(il, i)
    1857         ! MAF      sij(il,i,i)=1.0
     1924! MAF      sij(il,i,i)=1.0
    18581925        sij(il, i, i) = 0.0
    18591926      END IF
    … …  
    18611928  END DO
    18621929
    1863   ! AC!      do j=1,ntra
    1864   ! AC!       do i=minorig+1,nl
    1865   ! AC!        do il=1,ncum
    1866   ! AC!         if (i.ge.icb(il) .and. i.le.inb(il) .and. nent(il,i).eq.0)
    1867   ! then
    1868   ! AC!          traent(il,i,i,j)=tra(il,nk(il),j)
    1869   ! AC!         endif
    1870   ! AC!        enddo
    1871   ! AC!       enddo
    1872   ! AC!      enddo
     1930!AC!      do j=1,ntra
     1931!AC!       do i=minorig+1,nl
     1932!AC!        do il=1,ncum
     1933!AC!         if (i.ge.icb(il) .and. i.le.inb(il) .and. nent(il,i).eq.0) then
     1934!AC!          traent(il,i,i,j)=tra(il,nk(il),j)
     1935!AC!         endif
     1936!AC!        enddo
     1937!AC!       enddo
     1938!AC!      enddo
    18731939
    18741940  DO j = minorig, nl
    18751941    DO i = minorig, nl
    18761942      DO il = 1, ncum
    1877         IF ((j>=(icb(il)-1)) .AND. (j<=inb(il)) .AND. (i>=icb(il)) .AND. (i<= &
    1878             inb(il))) THEN
     1943        IF ((j>=(icb(il)-1)) .AND. (j<=inb(il)) .AND. (i>=icb(il)) .AND. (i<=inb(il))) THEN
    18791944          sigij(il, i, j) = sij(il, i, j)
    18801945        END IF
    … …  
    18821947    END DO
    18831948  END DO
    1884   ! @      enddo
    1885 
    1886   ! @170   continue
    1887 
    1888   ! =====================================================================
    1889   ! ---  NORMALIZE ENTRAINED AIR MASS FLUXES
    1890   ! ---  TO REPRESENT EQUAL PROBABILITIES OF MIXING
    1891   ! =====================================================================
     1949! @      enddo
     1950
     1951! @170   continue
     1952
     1953! =====================================================================
     1954! ---  NORMALIZE ENTRAINED AIR MASS FLUXES
     1955! ---  TO REPRESENT EQUAL PROBABILITIES OF MIXING
     1956! =====================================================================
    18921957
    18931958  CALL zilch(asum, nloc*nd)
    … …  
    19151980        IF (cvflag_ice) THEN
    19161981
    1917           anum = h(il, i) - hp(il, i) - (lv(il,i)+frac(il,i)*lf(il,i))*(qp-rs &
    1918             (il,i)) + (cpv-cpd)*t(il, i)*(qp-rr(il,i))
    1919           denom = h(il, i) - hp(il, i) + (lv(il,i)+frac(il,i)*lf(il,i))*(rr( &
    1920             il,i)-qp) + (cpd-cpv)*t(il, i)*(rr(il,i)-qp)
     1982          anum = h(il, i) - hp(il, i) - (lv(il,i)+frac(il,i)*lf(il,i))* &
     1983                       (qp-rs(il,i)) + (cpv-cpd)*t(il, i)*(qp-rr(il,i))
     1984          denom = h(il, i) - hp(il, i) + (lv(il,i)+frac(il,i)*lf(il,i))* &
     1985                       (rr(il,i)-qp) + (cpd-cpv)*t(il, i)*(rr(il,i)-qp)
    19211986        ELSE
    19221987
    19231988          anum = h(il, i) - hp(il, i) - lv(il, i)*(qp-rs(il,i)) + &
    1924             (cpv-cpd)*t(il, i)*(qp-rr(il,i))
     1989                       (cpv-cpd)*t(il, i)*(qp-rr(il,i))
    19251990          denom = h(il, i) - hp(il, i) + lv(il, i)*(rr(il,i)-qp) + &
    1926             (cpd-cpv)*t(il, i)*(rr(il,i)-qp)
     1991                       (cpd-cpv)*t(il, i)*(rr(il,i)-qp)
    19271992        END IF
    19281993
    … …  
    19402005      num2 = 0
    19412006      DO il = 1, ncum
    1942         IF (i>=icb(il) .AND. i<=inb(il) .AND. j>=(icb( &
    1943           il)-1) .AND. j<=inb(il) .AND. lwork(il)) num2 = num2 + 1
     2007        IF (i>=icb(il) .AND. i<=inb(il) .AND. &
     2008            j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
     2009            lwork(il)) num2 = num2 + 1
    19442010      END DO
    19452011      IF (num2<=0) GO TO 175
    19462012
    19472013      DO il = 1, ncum
    1948         IF (i>=icb(il) .AND. i<=inb(il) .AND. j>=(icb( &
    1949             il)-1) .AND. j<=inb(il) .AND. lwork(il)) THEN
     2014        IF (i>=icb(il) .AND. i<=inb(il) .AND. &
     2015            j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
     2016            lwork(il)) THEN
    19502017
    19512018          IF (sij(il,i,j)>1.0E-16 .AND. sij(il,i,j)<0.95) THEN
    … …  
    19882055    DO j = minorig, nl
    19892056      DO il = 1, ncum
    1990         IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. j>=(icb( &
    1991             il)-1) .AND. j<=inb(il)) THEN
     2057        IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
     2058            j>=(icb(il)-1) .AND. j<=inb(il)) THEN
    19922059          ment(il, i, j) = ment(il, i, j)*asij(il)
    19932060        END IF
    … …  
    19972064    DO j = minorig, nl
    19982065      DO il = 1, ncum
    1999         IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. j>=(icb( &
    2000             il)-1) .AND. j<=inb(il)) THEN
     2066        IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
     2067            j>=(icb(il)-1) .AND. j<=inb(il)) THEN
    20012068          asum(il, i) = asum(il, i) + ment(il, i, j)
    20022069          ment(il, i, j) = ment(il, i, j)*sig(il, j)
    … …  
    20152082    DO j = minorig, nl
    20162083      DO il = 1, ncum
    2017         IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. j>=(icb( &
    2018             il)-1) .AND. j<=inb(il)) THEN
     2084        IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
     2085            j>=(icb(il)-1) .AND. j<=inb(il)) THEN
    20192086          ment(il, i, j) = ment(il, i, j)*asum(il, i)*bsum(il, i)
    20202087        END IF
    … …  
    20242091    DO j = minorig, nl
    20252092      DO il = 1, ncum
    2026         IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. j>=(icb( &
    2027             il)-1) .AND. j<=inb(il)) THEN
     2093        IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
     2094            j>=(icb(il)-1) .AND. j<=inb(il)) THEN
    20282095          csum(il, i) = csum(il, i) + ment(il, i, j)
    20292096        END IF
    … …  
    20402107        vent(il, i, i) = vnk(il)
    20412108        elij(il, i, i) = clw(il, i)
    2042         ! MAF        sij(il,i,i)=1.0
     2109! MAF        sij(il,i,i)=1.0
    20432110        sij(il, i, i) = 0.0
    20442111      END IF
    20452112    END DO ! il
    20462113
    2047     ! AC!      do j=1,ntra
    2048     ! AC!       do il=1,ncum
    2049     ! AC!        if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il)
    2050     ! AC!     :     .and. csum(il,i).lt.m(il,i) ) then
    2051     ! AC!         traent(il,i,i,j)=tra(il,nk(il),j)
    2052     ! AC!        endif
    2053     ! AC!       enddo
    2054     ! AC!      enddo
     2114!AC!      do j=1,ntra
     2115!AC!       do il=1,ncum
     2116!AC!        if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il)
     2117!AC!     :     .and. csum(il,i).lt.m(il,i) ) then
     2118!AC!         traent(il,i,i,j)=tra(il,nk(il),j)
     2119!AC!        endif
     2120!AC!       enddo
     2121!AC!      enddo
    20552122789 END DO
    20562123
    2057   ! MAF: renormalisation de MENT
     2124! MAF: renormalisation de MENT
    20582125  CALL zilch(zm, nloc*na)
    20592126  DO jm = 1, nd
    … …  
    20872154END SUBROUTINE cv3_mixing
    20882155
    2089 SUBROUTINE cv3_unsat(nloc, ncum, nd, na, ntra, icb, inb, iflag, t, rr, rs, &
    2090     gz, u, v, tra, p, ph, th, tv, lv, lf, cpn, ep, sigp, clw, m, ment, elij, &
    2091     delt, plcl, coef_clos, mp, rp, up, vp, trap, wt, water, evap, fondue, &
    2092     ice, faci, b, sigd, wdtraina, wdtrainm) ! RomP
     2156SUBROUTINE cv3_unsat(nloc, ncum, nd, na, ntra, icb, inb, iflag, &
     2157                     t, rr, rs, gz, u, v, tra, p, ph, &
     2158                     th, tv, lv, lf, cpn, ep, sigp, clw, &
     2159                     m, ment, elij, delt, plcl, coef_clos, &
     2160                     mp, rp, up, vp, trap, wt, water, evap, fondue, ice, &
     2161                     faci, b, sigd, &
     2162                     wdtrainA, wdtrainM)                                      ! RomP
    20932163  IMPLICIT NONE
    20942164
    … …  
    20982168  include "cvflag.h"
    20992169
    2100   ! inputs:
     2170!inputs:
    21012171  INTEGER ncum, nd, na, ntra, nloc
    21022172  INTEGER icb(nloc), inb(nloc)
    … …  
    21122182  REAL coef_clos(nloc)
    21132183
    2114   ! input/output
     2184!input/output
    21152185  INTEGER iflag(nloc)
    21162186
    2117   ! outputs:
     2187!outputs:
    21182188  REAL mp(nloc, na), rp(nloc, na), up(nloc, na), vp(nloc, na)
    21192189  REAL water(nloc, na), evap(nloc, na), wt(nloc, na)
    … …  
    21212191  REAL trap(nloc, na, ntra)
    21222192  REAL b(nloc, na), sigd(nloc)
    2123   ! 25/08/10 - RomP---- ajout des masses precipitantes ejectees
    2124   ! lascendance adiabatique et des flux melanges Pa et Pm.
    2125   ! Distinction des wdtrain
    2126   ! Pa = wdtrainA     Pm = wdtrainM
    2127   REAL wdtraina(nloc, na), wdtrainm(nloc, na)
    2128 
    2129   ! local variables
     2193! 25/08/10 - RomP---- ajout des masses precipitantes ejectees
     2194! de l ascendance adiabatique et des flux melanges Pa et Pm.
     2195! Distinction des wdtrain
     2196! Pa = wdtrainA     Pm = wdtrainM
     2197  REAL wdtrainA(nloc, na), wdtrainM(nloc, na)
     2198
     2199!local variables
    21302200  INTEGER i, j, k, il, num1, ndp1
    21312201  REAL tinv, delti, coef
    … …  
    21432213
    21442214
    2145   ! ------------------------------------------------------
     2215! ------------------------------------------------------
    21462216
    21472217  delti = 1./delt
    … …  
    21702240    END DO
    21712241  END DO
    2172   ! AC!        do k=1,ntra
    2173   ! AC!         do i=1,nd
    2174   ! AC!          do il=1,ncum
    2175   ! AC!           trap(il,i,k)=tra(il,i,k)
    2176   ! AC!          enddo
    2177   ! AC!         enddo
    2178   ! AC!        enddo
    2179   ! ! RomP >>>
     2242!AC!        do k=1,ntra
     2243!AC!         do i=1,nd
     2244!AC!          do il=1,ncum
     2245!AC!           trap(il,i,k)=tra(il,i,k)
     2246!AC!          enddo
     2247!AC!         enddo
     2248!AC!        enddo
     2249!! RomP >>>
    21802250  DO i = 1, nd
    21812251    DO il = 1, ncum
    2182       wdtraina(il, i) = 0.0
    2183       wdtrainm(il, i) = 0.0
    2184     END DO
    2185   END DO
    2186   ! ! RomP <<<
    2187 
    2188   ! ***  check whether ep(inb)=0, if so, skip precipitating    ***
    2189   ! ***             downdraft calculation                      ***
     2252      wdtrainA(il, i) = 0.0
     2253      wdtrainM(il, i) = 0.0
     2254    END DO
     2255  END DO
     2256!! RomP <<<
     2257
     2258! ***  check whether ep(inb)=0, if so, skip precipitating    ***
     2259! ***             downdraft calculation                      ***
    21902260
    21912261
    21922262  DO il = 1, ncum
    2193     ! !          lwork(il)=.TRUE.
    2194     ! !          if(ep(il,inb(il)).lt.0.0001)lwork(il)=.FALSE.
     2263!!          lwork(il)=.TRUE.
     2264!!          if(ep(il,inb(il)).lt.0.0001)lwork(il)=.FALSE.
    21952265    lwork(il) = ep(il, inb(il)) >= 0.0001
    21962266  END DO
    21972267
    2198   ! ***  Set the fractionnal area sigd of precipitating downdraughts
     2268! ***  Set the fractionnal area sigd of precipitating downdraughts
    21992269  DO il = 1, ncum
    22002270    sigd(il) = sigdz*coef_clos(il)
    … …  
    22022272
    22032273
    2204   ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    2205 
    2206   ! ***                    begin downdraft loop                    ***
    2207 
    2208   ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     2274! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     2275!
     2276! ***                    begin downdraft loop                    ***
     2277!
     2278! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    22092279
    22102280  DO i = nl + 1, 1, -1
    … …  
    22192289
    22202290
    2221     ! ***  integrate liquid water equation to find condensed water   ***
    2222     ! ***                and condensed water flux                    ***
    2223 
    2224 
    2225     ! ***              calculate detrained precipitation             ***
     2291! ***  integrate liquid water equation to find condensed water   ***
     2292! ***                and condensed water flux                    ***
     2293!
     2294!
     2295! ***              calculate detrained precipitation             ***
    22262296
    22272297    DO il = 1, ncum
    … …  
    22292299        IF (cvflag_grav) THEN
    22302300          wdtrain(il) = grav*ep(il, i)*m(il, i)*clw(il, i)
    2231           wdtraina(il, i) = wdtrain(il)/grav !   Pa   RomP
     2301          wdtrainA(il, i) = wdtrain(il)/grav                        !   Pa   RomP
    22322302        ELSE
    22332303          wdtrain(il) = 10.0*ep(il, i)*m(il, i)*clw(il, i)
    2234           wdtraina(il, i) = wdtrain(il)/10. !   Pa   RomP
     2304          wdtrainA(il, i) = wdtrain(il)/10.                        !   Pa   RomP
    22352305        END IF
    22362306      END IF
    … …  
    22452315            IF (cvflag_grav) THEN
    22462316              wdtrain(il) = wdtrain(il) + grav*awat*ment(il, j, i)
    2247               wdtrainm(il, i) = wdtrain(il)/grav - wdtraina(il, i) !   Pm  RomP
     2317              wdtrainM(il, i) = wdtrain(il)/grav - wdtrainA(il, i) !   Pm  RomP
    22482318            ELSE
    22492319              wdtrain(il) = wdtrain(il) + 10.0*awat*ment(il, j, i)
    2250               wdtrainm(il, i) = wdtrain(il)/10. - wdtraina(il, i) !   Pm  RomP
     2320              wdtrainM(il, i) = wdtrain(il)/10. - wdtrainA(il, i)  !   Pm  RomP
    22512321            END IF
    22522322          END IF
    … …  
    22562326
    22572327
    2258     ! ***    find rain water and evaporation using provisional   ***
    2259     ! ***              estimates of rp(i)and rp(i-1)             ***
     2328! ***    find rain water and evaporation using provisional   ***
     2329! ***              estimates of rp(i)and rp(i-1)             ***
    22602330
    22612331
    … …  
    22832353          END IF
    22842354
    2285           rp(il, i) = rp(il, i+1) + (cpd*(t(il,i+1)-t(il, &
    2286             i))+gz(il,i+1)-gz(il,i))/lv(il, i)
     2355          rp(il, i) = rp(il, i+1) + &
     2356                      (cpd*(t(il,i+1)-t(il,i))+gz(il,i+1)-gz(il,i))/lv(il, i)
    22872357          rp(il, i) = 0.5*(rp(il,i)+rr(il,i))
    22882358        END IF
    … …  
    22962366          afac = p(il, 1)*(rs(il,1)-rp(il,1))/(1.0E4+2000.0*p(il,1)*rs(il,1))
    22972367          IF (cvflag_ice) THEN
    2298             afac1 = p(il, i)*(rs(il,1)-rp(il,1))/(1.0E4+2000.0*p(il,1)*rs(il, &
    2299               1))
     2368            afac1 = p(il, i)*(rs(il,1)-rp(il,1))/(1.0E4+2000.0*p(il,1)*rs(il,1))
    23002369          END IF
    23012370        ELSE
    2302           rp(il, i-1) = rp(il, i) + (cpd*(t(il,i)-t(il, &
    2303             i-1))+gz(il,i)-gz(il,i-1))/lv(il, i)
     2371          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)
    23042372          rp(il, i-1) = 0.5*(rp(il,i-1)+rr(il,i-1))
    23052373          rp(il, i-1) = amin1(rp(il,i-1), rs(il,i-1))
    23062374          rp(il, i-1) = max(rp(il,i-1), 0.0)
    2307           afac1 = p(il, i)*(rs(il,i)-rp(il,i))/(1.0E4+2000.0*p(il,i)*rs(il,i) &
    2308             )
    2309           afac2 = p(il, i-1)*(rs(il,i-1)-rp(il,i-1))/ &
    2310             (1.0E4+2000.0*p(il,i-1)*rs(il,i-1))
     2375          afac1 = p(il, i)*(rs(il,i)-rp(il,i))/(1.0E4+2000.0*p(il,i)*rs(il,i))
     2376          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))
    23112377          afac = 0.5*(afac1+afac2)
    23122378        END IF
    … …  
    23152381        bfac = 1./(sigd(il)*wt(il,i))
    23162382
    2317         ! jyg1
    2318         ! cc        sigt=1.0
    2319         ! cc        if(i.ge.icb)sigt=sigp(i)
    2320         ! prise en compte de la variation progressive de sigt dans
    2321         ! les couches icb et icb-1:
    2322         ! pour plcl<ph(i+1), pr1=0 & pr2=1
    2323         ! pour plcl>ph(i),   pr1=1 & pr2=0
    2324         ! pour ph(i+1)<plcl<ph(i), pr1 est la proportion a cheval
    2325         ! sur le nuage, et pr2 est la proportion sous la base du
    2326         ! nuage.
     2383!JYG1
     2384! cc        sigt=1.0
     2385! cc        if(i.ge.icb)sigt=sigp(i)
     2386! prise en compte de la variation progressive de sigt dans
     2387! les couches icb et icb-1:
     2388! pour plcl<ph(i+1), pr1=0 & pr2=1
     2389! pour plcl>ph(i),   pr1=1 & pr2=0
     2390! pour ph(i+1)<plcl<ph(i), pr1 est la proportion a cheval
     2391! sur le nuage, et pr2 est la proportion sous la base du
     2392! nuage.
    23272393        pr1 = (plcl(il)-ph(il,i+1))/(ph(il,i)-ph(il,i+1))
    23282394        pr1 = max(0., min(1.,pr1))
    … …  
    23302396        pr2 = max(0., min(1.,pr2))
    23312397        sigt = sigp(il, i)*pr1 + pr2
    2332         ! jyg2
    2333 
    2334         ! jyg----
    2335         ! b6 = bfac*100.*sigd(il)*(ph(il,i)-ph(il,i+1))*sigt*afac
    2336         ! c6 = water(il,i+1) + wdtrain(il)*bfac
    2337         ! c6 = prec(il,i+1) + wdtrain(il)*bfac
    2338         ! revap=0.5*(-b6+sqrt(b6*b6+4.*c6))
    2339         ! evap(il,i)=sigt*afac*revap
    2340         ! water(il,i)=revap*revap
    2341         ! prec(il,i)=revap*revap
    2342         ! c        print *,' i,b6,c6,revap,evap(il,i),water(il,i),wdtrain(il)
    2343         ! ',
    2344         ! c     $            i,b6,c6,revap,evap(il,i),water(il,i),wdtrain(il)
    2345         ! c---end jyg---
    2346 
    2347         ! --------retour à la formulation originale d'Emanuel.
     2398!JYG2
     2399
     2400!JYG----
     2401!    b6 = bfac*100.*sigd(il)*(ph(il,i)-ph(il,i+1))*sigt*afac
     2402!    c6 = water(il,i+1) + wdtrain(il)*bfac
     2403!    c6 = prec(il,i+1) + wdtrain(il)*bfac
     2404!    revap=0.5*(-b6+sqrt(b6*b6+4.*c6))
     2405!    evap(il,i)=sigt*afac*revap
     2406!    water(il,i)=revap*revap
     2407!    prec(il,i)=revap*revap
     2408!!        print *,' i,b6,c6,revap,evap(il,i),water(il,i),wdtrain(il) ', &
     2409!!                 i,b6,c6,revap,evap(il,i),water(il,i),wdtrain(il)
     2410!!---end jyg---
     2411
     2412! --------retour à la formulation originale d'Emanuel.
    23482413        IF (cvflag_ice) THEN
    23492414
    2350           ! b6=bfac*50.*sigd(il)*(ph(il,i)-ph(il,i+1))*sigt*afac
    2351           ! c6=prec(il,i+1)+bfac*wdtrain(il)
    2352           ! :    -50.*sigd(il)*bfac*(ph(il,i)-ph(il,i+1))*evap(il,i+1)
    2353           ! if(c6.gt.0.0)then
    2354           ! revap=0.5*(-b6+sqrt(b6*b6+4.*c6))
    2355 
    2356           ! JAM  Attention: evap=sigt*E
    2357           ! Modification: evap devient l'évaporation en milieu de couche
    2358           ! car nécessaire dans cv3_yield
    2359           ! Du coup, il faut modifier pas mal d'équations...
    2360           ! et l'expression de afac qui devient afac1
    2361           ! revap=sqrt((prec(i+1)+prec(i))/2)
     2415!  b6=bfac*50.*sigd(il)*(ph(il,i)-ph(il,i+1))*sigt*afac
     2416!   c6=prec(il,i+1)+bfac*wdtrain(il) &
     2417!       -50.*sigd(il)*bfac*(ph(il,i)-ph(il,i+1))*evap(il,i+1)
     2418!  if(c6.gt.0.0)then
     2419!  revap=0.5*(-b6+sqrt(b6*b6+4.*c6))
     2420
     2421!JAM  Attention: evap=sigt*E
     2422!    Modification: evap devient l'évaporation en milieu de couche
     2423!    car nécessaire dans cv3_yield
     2424!    Du coup, il faut modifier pas mal d'équations...
     2425!    et l'expression de afac qui devient afac1
     2426!    revap=sqrt((prec(i+1)+prec(i))/2)
    23622427
    23632428          b6 = bfac*50.*sigd(il)*(ph(il,i)-ph(il,i+1))*sigt*afac1
    23642429          c6 = prec(il, i+1) + 0.5*bfac*wdtrain(il)
    2365           ! print *,'bfac,sigd(il),sigt,afac1 ',bfac,sigd(il),sigt,afac1
    2366           ! print *,'prec(il,i+1),wdtrain(il) ',prec(il,i+1),wdtrain(il)
    2367           ! print *,'b6,c6,b6*b6+4.*c6 ',b6,c6,b6*b6+4.*c6
     2430! print *,'bfac,sigd(il),sigt,afac1 ',bfac,sigd(il),sigt,afac1
     2431! print *,'prec(il,i+1),wdtrain(il) ',prec(il,i+1),wdtrain(il)
     2432! print *,'b6,c6,b6*b6+4.*c6 ',b6,c6,b6*b6+4.*c6
    23682433          IF (c6>b6*b6+1.E-20) THEN
    23692434            revap = 2.*c6/(b6+sqrt(b6*b6+4.*c6))
    … …  
    23722437          END IF
    23732438          prec(il, i) = max(0., 2.*revap*revap-prec(il,i+1))
    2374           ! print*,prec(il,i),'neige'
    2375 
    2376           ! jyg    Dans sa formulation originale, Emanuel calcule
    2377           ! l'evaporation par:
    2378           ! c             evap(il,i)=sigt*afac*revap
    2379           ! ce qui n'est pas correct. Dans cv_routines, la formulation a été
    2380           ! modifiee.
    2381           ! Ici,l'evaporation evap est simplement calculee par l'equation de
    2382           ! conservation.
    2383           ! prec(il,i)=revap*revap
    2384           ! else
    2385           ! jyg----   Correction : si c6 <= 0, water(il,i)=0.
    2386           ! prec(il,i)=0.
    2387           ! endif
    2388 
    2389           ! jyg---   Dans tous les cas, evaporation = [tt ce qui entre dans
    2390           ! la couche i]
    2391           ! moins [tt ce qui sort de la couche i]
    2392           ! print *, 'evap avec ice'
    2393           evap(il, i) = (wdtrain(il)+sigd(il)*wt(il,i)*(prec(il,i+1)-prec(il, &
    2394             i)))/(sigd(il)*(ph(il,i)-ph(il,i+1))*100.)
    2395 
    2396           d6 = bfac*wdtrain(il) - 100.*sigd(il)*bfac*(ph(il,i)-ph(il,i+1))* &
    2397             evap(il, i)
     2439! print*,prec(il,i),'neige'
     2440
     2441!JYG    Dans sa formulation originale, Emanuel calcule l'evaporation par:
     2442! c             evap(il,i)=sigt*afac*revap
     2443! ce qui n'est pas correct. Dans cv_routines, la formulation a été modifiee.
     2444! Ici,l'evaporation evap est simplement calculee par l'equation de
     2445! conservation.
     2446! prec(il,i)=revap*revap
     2447! else
     2448!JYG----   Correction : si c6 <= 0, water(il,i)=0.
     2449! prec(il,i)=0.
     2450! endif
     2451
     2452!JYG---   Dans tous les cas, evaporation = [tt ce qui entre dans la couche i]
     2453! moins [tt ce qui sort de la couche i]
     2454! print *, 'evap avec ice'
     2455          evap(il, i) = (wdtrain(il)+sigd(il)*wt(il,i)*(prec(il,i+1)-prec(il,i))) / &
     2456                        (sigd(il)*(ph(il,i)-ph(il,i+1))*100.)
     2457
     2458          d6 = bfac*wdtrain(il) - 100.*sigd(il)*bfac*(ph(il,i)-ph(il,i+1))*evap(il, i)
    23982459          e6 = bfac*wdtrain(il)
    23992460          f6 = -100.*sigd(il)*bfac*(ph(il,i)-ph(il,i+1))*evap(il, i)
    … …  
    24152476          END IF
    24162477
    2417           ! water(il,i)=water(il,i+1)+(1.-fraci(il,i))*e6+(1.-faci(il,i))*f6
    2418           ! water(il,i)=max(water(il,i),0.)
    2419           ! ice(il,i)=ice(il,i+1)+fraci(il,i)*e6+faci(il,i)*f6
    2420           ! ice(il,i)=max(ice(il,i),0.)
    2421           ! fondue(il,i)=ice(il,i)*thaw
    2422           ! water(il,i)=water(il,i)+fondue(il,i)
    2423           ! ice(il,i)=ice(il,i)-fondue(il,i)
    2424 
    2425           ! if((water(il,i)+ice(il,i)).lt.1.e-30)then
    2426           ! faci(il,i)=0.
    2427           ! else
    2428           ! faci(il,i)=ice(il,i)/(water(il,i)+ice(il,i))
    2429           ! endif
     2478!          water(il,i)=water(il,i+1)+(1.-fraci(il,i))*e6+(1.-faci(il,i))*f6
     2479!          water(il,i)=max(water(il,i),0.)
     2480!          ice(il,i)=ice(il,i+1)+fraci(il,i)*e6+faci(il,i)*f6
     2481!          ice(il,i)=max(ice(il,i),0.)
     2482!          fondue(il,i)=ice(il,i)*thaw
     2483!          water(il,i)=water(il,i)+fondue(il,i)
     2484!          ice(il,i)=ice(il,i)-fondue(il,i)
     2485           
     2486!          if((water(il,i)+ice(il,i)).lt.1.e-30)then
     2487!            faci(il,i)=0.
     2488!          else
     2489!            faci(il,i)=ice(il,i)/(water(il,i)+ice(il,i))
     2490!          endif
    24302491
    24312492        ELSE
    24322493          b6 = bfac*50.*sigd(il)*(ph(il,i)-ph(il,i+1))*sigt*afac
    2433           c6 = water(il, i+1) + bfac*wdtrain(il) - 50.*sigd(il)*bfac*(ph(il,i &
    2434             )-ph(il,i+1))*evap(il, i+1)
     2494          c6 = water(il, i+1) + bfac*wdtrain(il) - &
     2495               50.*sigd(il)*bfac*(ph(il,i)-ph(il,i+1))*evap(il, i+1)
    24352496          IF (c6>0.0) THEN
    24362497            revap = 0.5*(-b6+sqrt(b6*b6+4.*c6))
    … …  
    24392500            water(il, i) = 0.
    24402501          END IF
    2441           ! print *, 'evap sans ice'
    2442           evap(il, i) = (wdtrain(il)+sigd(il)*wt(il,i)*(water(il, &
    2443             i+1)-water(il,i)))/(sigd(il)*(ph(il,i)-ph(il,i+1))*100.)
     2502! print *, 'evap sans ice'
     2503          evap(il, i) = (wdtrain(il)+sigd(il)*wt(il,i)*(water(il,i+1)-water(il,i)))/ &
     2504                        (sigd(il)*(ph(il,i)-ph(il,i+1))*100.)
    24442505
    24452506        END IF
    24462507      END IF !(i.le.inb(il) .and. lwork(il))
    24472508    END DO
    2448     ! ----------------------------------------------------------------
    2449 
    2450     ! cc
    2451     ! ***  calculate precipitating downdraft mass flux under     ***
    2452     ! ***              hydrostatic approximation                 ***
     2509! ----------------------------------------------------------------
     2510
     2511! cc
     2512! ***  calculate precipitating downdraft mass flux under     ***
     2513! ***              hydrostatic approximation                 ***
    24532514
    24542515    DO il = 1, ncum
    … …  
    24592520        IF (cvflag_ice) THEN
    24602521          IF (cvflag_grav) THEN
    2461             mp(il, i) = 100.*ginv*(lvcp(il,i)*sigd(il)*tevap(il)*(p(il, &
    2462               i-1)-p(il,i))/delth+lfcp(il,i)*sigd(il)*faci(il,i)*tevap(il)*(p &
    2463               (il,i-1)-p(il,i))/delth+lfcp(il,i)*sigd(il)*wt(il,i)/100.* &
    2464               fondue(il,i)*(p(il,i-1)-p(il,i))/delth/(ph(il,i)-ph(il,i+1)))
     2522            mp(il, i) = 100.*ginv*(lvcp(il,i)*sigd(il)*tevap(il)* &
     2523                                               (p(il,i-1)-p(il,i))/delth + &
     2524                                   lfcp(il,i)*sigd(il)*faci(il,i)*tevap(il)* &
     2525                                               (p(il,i-1)-p(il,i))/delth + &
     2526                                   lfcp(il,i)*sigd(il)*wt(il,i)/100.*fondue(il,i)* &
     2527                                               (p(il,i-1)-p(il,i))/delth/(ph(il,i)-ph(il,i+1)))
    24652528          ELSE
    2466             mp(il, i) = 10.*(lvcp(il,i)*sigd(il)*tevap(il)*(p(il,i-1)-p(il, &
    2467               i))/delth+lfcp(il,i)*sigd(il)*faci(il,i)*tevap(il)*(p(il, &
    2468               i-1)-p(il,i))/delth+lfcp(il,i)*sigd(il)*wt(il,i)/100.*fondue(il &
    2469               ,i)*(p(il,i-1)-p(il,i))/delth/(ph(il,i)-ph(il,i+1)))
     2529            mp(il, i) = 10.*(lvcp(il,i)*sigd(il)*tevap(il)* &
     2530                                                (p(il,i-1)-p(il,i))/delth + &
     2531                             lfcp(il,i)*sigd(il)*faci(il,i)*tevap(il)* &
     2532                                                (p(il,i-1)-p(il,i))/delth + &
     2533                             lfcp(il,i)*sigd(il)*wt(il,i)/100.*fondue(il,i)* &
     2534                                                (p(il,i-1)-p(il,i))/delth/(ph(il,i)-ph(il,i+1)))
    24702535
    24712536          END IF
    … …  
    24732538          IF (cvflag_grav) THEN
    24742539            mp(il, i) = 100.*ginv*lvcp(il, i)*sigd(il)*tevap(il)* &
    2475               (p(il,i-1)-p(il,i))/delth
     2540                                                (p(il,i-1)-p(il,i))/delth
    24762541          ELSE
    24772542            mp(il, i) = 10.*lvcp(il, i)*sigd(il)*tevap(il)* &
    2478               (p(il,i-1)-p(il,i))/delth
     2543                                                (p(il,i-1)-p(il,i))/delth
    24792544          END IF
    24802545
    … …  
    24832548      END IF !(i.le.inb(il) .and. lwork(il) .and. i.ne.1)
    24842549    END DO
    2485     ! ----------------------------------------------------------------
    2486 
    2487     ! ***           if hydrostatic assumption fails,             ***
    2488     ! ***   solve cubic difference equation for downdraft theta  ***
    2489     ! ***  and mass flux from two simultaneous differential eqns ***
     2550! ----------------------------------------------------------------
     2551
     2552! ***           if hydrostatic assumption fails,             ***
     2553! ***   solve cubic difference equation for downdraft theta  ***
     2554! ***  and mass flux from two simultaneous differential eqns ***
    24902555
    24912556    DO il = 1, ncum
    … …  
    24932558
    24942559        amfac = sigd(il)*sigd(il)*70.0*ph(il, i)*(p(il,i-1)-p(il,i))* &
    2495           (th(il,i)-th(il,i-1))/(tv(il,i)*th(il,i))
     2560                         (th(il,i)-th(il,i-1))/(tv(il,i)*th(il,i))
    24962561        amp2 = abs(mp(il,i+1)*mp(il,i+1)-mp(il,i)*mp(il,i))
    24972562
    24982563        IF (amp2>(0.1*amfac)) THEN
    24992564          xf = 100.0*sigd(il)*sigd(il)*sigd(il)*(ph(il,i)-ph(il,i+1))
    2500           tf = b(il, i) - 5.0*(th(il,i)-th(il,i-1))*t(il, i)/(lvcp(il,i)*sigd &
    2501             (il)*th(il,i))
     2565          tf = b(il, i) - 5.0*(th(il,i)-th(il,i-1))*t(il, i) / &
     2566                              (lvcp(il,i)*sigd(il)*th(il,i))
    25022567          af = xf*tf + mp(il, i+1)*mp(il, i+1)*tinv
    25032568
    25042569          IF (cvflag_ice) THEN
    25052570            bf = 2.*(tinv*mp(il,i+1))**3 + tinv*mp(il, i+1)*xf*tf + &
    2506               50.*(p(il,i-1)-p(il,i))*xf*(tevap(il)*(1.+(lf(il,i)/lv(il,i))* &
    2507               faci(il,i))+(lf(il,i)/lv(il,i))*wt(il,i)/100.*fondue(il,i)/(ph( &
    2508               il,i)-ph(il,i+1)))
     2571                 50.*(p(il,i-1)-p(il,i))*xf*(tevap(il)*(1.+(lf(il,i)/lv(il,i))*faci(il,i)) + &
     2572                (lf(il,i)/lv(il,i))*wt(il,i)/100.*fondue(il,i)/(ph(il,i)-ph(il,i+1)))
    25092573          ELSE
    25102574
    25112575            bf = 2.*(tinv*mp(il,i+1))**3 + tinv*mp(il, i+1)*xf*tf + &
    2512               50.*(p(il,i-1)-p(il,i))*xf*tevap(il)
     2576                                           50.*(p(il,i-1)-p(il,i))*xf*tevap(il)
    25132577          END IF
    25142578
    … …  
    25222586            IF ((0.5*bf-sru)<0.0) fac = -1.0
    25232587            mp(il, i) = mp(il, i+1)*tinv + (0.5*bf+sru)**tinv + &
    2524               fac*(abs(0.5*bf-sru))**tinv
     2588                                           fac*(abs(0.5*bf-sru))**tinv
    25252589          ELSE
    25262590            d = atan(2.*sqrt(-ur)/(bf+1.0E-28))
    … …  
    25322596          IF (cvflag_ice) THEN
    25332597            IF (cvflag_grav) THEN
    2534               ! jyg : il y a vraisemblablement une erreur dans la ligne 2
    2535               ! suivante:
    2536               ! il faut diviser par (mp(il,i)*sigd(il)*grav) et non par
    2537               ! (mp(il,i)+sigd(il)*0.1).
    2538               ! Et il faut bien revoir les facteurs 100.
    2539               b(il, i-1) = b(il, i) + 100.0*(p(il,i-1)-p(il,i))*(tevap(il)*( &
    2540                 1.+(lf(il,i)/lv(il,i))*faci(il,i))+(lf(il,i)/lv(il, &
    2541                 i))*wt(il,i)/100.*fondue(il,i)/(ph(il,i)-ph(il, &
    2542                 i+1)))/(mp(il,i)+sigd(il)*0.1) - 10.0*(th(il,i)-th(il,i-1))*t &
    2543                 (il, i)/(lvcp(il,i)*sigd(il)*th(il,i))
     2598!JYG : il y a vraisemblablement une erreur dans la ligne 2 suivante:
     2599! il faut diviser par (mp(il,i)*sigd(il)*grav) et non par (mp(il,i)+sigd(il)*0.1).
     2600! Et il faut bien revoir les facteurs 100.
     2601              b(il, i-1) = b(il, i) + 100.0*(p(il,i-1)-p(il,i))* &
     2602                           (tevap(il)*(1.+(lf(il,i)/lv(il,i))*faci(il,i)) + &
     2603                           (lf(il,i)/lv(il,i))*wt(il,i)/100.*fondue(il,i) / &
     2604                           (ph(il,i)-ph(il,i+1))) / &
     2605                           (mp(il,i)+sigd(il)*0.1) - &
     2606                           10.0*(th(il,i)-th(il,i-1))*t(il, i) / &
     2607                           (lvcp(il,i)*sigd(il)*th(il,i))
    25442608            ELSE
    2545               b(il, i-1) = b(il, i) + 100.0*(p(il,i-1)-p(il,i))*(tevap(il)*( &
    2546                 1.+(lf(il,i)/lv(il,i))*faci(il,i))+(lf(il,i)/lv(il, &
    2547                 i))*wt(il,i)/100.*fondue(il,i)/(ph(il,i)-ph(il, &
    2548                 i+1)))/(mp(il,i)+sigd(il)*0.1) - 10.0*(th(il,i)-th(il,i-1))*t &
    2549                 (il, i)/(lvcp(il,i)*sigd(il)*th(il,i))
     2609              b(il, i-1) = b(il, i) + 100.0*(p(il,i-1)-p(il,i))*&
     2610                           (tevap(il)*(1.+(lf(il,i)/lv(il,i))*faci(il,i)) + &
     2611                           (lf(il,i)/lv(il,i))*wt(il,i)/100.*fondue(il,i) / &
     2612                           (ph(il,i)-ph(il,i+1))) / &
     2613                           (mp(il,i)+sigd(il)*0.1) - &
     2614                           10.0*(th(il,i)-th(il,i-1))*t(il, i) / &
     2615                           (lvcp(il,i)*sigd(il)*th(il,i))
    25502616            END IF
    25512617          ELSE
    25522618            IF (cvflag_grav) THEN
    2553               b(il, i-1) = b(il, i) + 100.0*(p(il,i-1)-p(il,i))*tevap(il)/(mp &
    2554                 (il,i)+sigd(il)*0.1) - 10.0*(th(il,i)-th(il,i-1))*t(il, i)/( &
    2555                 lvcp(il,i)*sigd(il)*th(il,i))
     2619              b(il, i-1) = b(il, i) + 100.0*(p(il,i-1)-p(il,i))*tevap(il) / &
     2620                           (mp(il,i)+sigd(il)*0.1) - &
     2621                           10.0*(th(il,i)-th(il,i-1))*t(il, i) / &
     2622                           (lvcp(il,i)*sigd(il)*th(il,i))
    25562623            ELSE
    2557               b(il, i-1) = b(il, i) + 100.0*(p(il,i-1)-p(il,i))*tevap(il)/(mp &
    2558                 (il,i)+sigd(il)*0.1) - 10.0*(th(il,i)-th(il,i-1))*t(il, i)/( &
    2559                 lvcp(il,i)*sigd(il)*th(il,i))
     2624              b(il, i-1) = b(il, i) + 100.0*(p(il,i-1)-p(il,i))*tevap(il) / &
     2625                           (mp(il,i)+sigd(il)*0.1) - &
     2626                           10.0*(th(il,i)-th(il,i-1))*t(il, i) / &
     2627                           (lvcp(il,i)*sigd(il)*th(il,i))
    25602628            END IF
    25612629          END IF
    … …  
    25642632        END IF !(amp2.gt.(0.1*amfac))
    25652633
    2566         ! ***         limit magnitude of mp(i) to meet cfl condition      ***
     2634! ***         limit magnitude of mp(i) to meet cfl condition      ***
    25672635
    25682636        ampmax = 2.0*(ph(il,i)-ph(il,i+1))*delti
    … …  
    25712639        mp(il, i) = min(mp(il,i), ampmax)
    25722640
    2573         ! ***      force mp to decrease linearly to zero                 ***
    2574         ! ***       between cloud base and the surface                   ***
    2575 
    2576 
    2577         ! c      if(p(il,i).gt.p(il,icb(il)))then
    2578         ! c
    2579         ! mp(il,i)=mp(il,icb(il))*(p(il,1)-p(il,i))/(p(il,1)-p(il,icb(il)))
    2580         ! c      endif
     2641! ***      force mp to decrease linearly to zero                 ***
     2642! ***       between cloud base and the surface                   ***
     2643
     2644
     2645! c      if(p(il,i).gt.p(il,icb(il)))then
     2646! c       mp(il,i)=mp(il,icb(il))*(p(il,1)-p(il,i))/(p(il,1)-p(il,icb(il)))
     2647! c      endif
    25812648        IF (ph(il,i)>0.9*plcl(il)) THEN
    25822649          mp(il, i) = mp(il, i)*(ph(il,1)-ph(il,i))/(ph(il,1)-0.9*plcl(il))
    … …  
    25852652      END IF ! (i.le.inb(il) .and. lwork(il) .and. i.ne.1)
    25862653    END DO
    2587     ! ----------------------------------------------------------------
    2588 
    2589     ! ***       find mixing ratio of precipitating downdraft     ***
     2654! ----------------------------------------------------------------
     2655
     2656! ***       find mixing ratio of precipitating downdraft     ***
    25902657
    25912658    DO il = 1, ncum
    … …  
    26032670
    26042671          IF (cvflag_grav) THEN
    2605             rp(il, i) = rp(il, i+1)*mp(il, i+1) + rr(il, i)*(mp(il,i)-mp(il,i &
    2606               +1)) + 100.*ginv*0.5*sigd(il)*(ph(il,i)-ph(il,i+1))*(evap(il,i+ &
    2607               1)+evap(il,i))
     2672            rp(il, i) = rp(il, i+1)*mp(il, i+1) + rr(il, i)*(mp(il,i)-mp(il,i+1)) + &
     2673              100.*ginv*0.5*sigd(il)*(ph(il,i)-ph(il,i+1))*(evap(il,i+1)+evap(il,i))
    26082674          ELSE
    2609             rp(il, i) = rp(il, i+1)*mp(il, i+1) + rr(il, i)*(mp(il,i)-mp(il,i &
    2610               +1)) + 5.*sigd(il)*(ph(il,i)-ph(il,i+1))*(evap(il,i+1)+evap(il, &
    2611               i))
     2675            rp(il, i) = rp(il, i+1)*mp(il, i+1) + rr(il, i)*(mp(il,i)-mp(il,i+1)) + &
     2676              5.*sigd(il)*(ph(il,i)-ph(il,i+1))*(evap(il,i+1)+evap(il,i))
    26122677          END IF
    26132678          rp(il, i) = rp(il, i)/mp(il, i)
    2614           up(il, i) = up(il, i+1)*mp(il, i+1) + u(il, i)*(mp(il,i)-mp(il,i+1) &
    2615             )
     2679          up(il, i) = up(il, i+1)*mp(il, i+1) + u(il, i)*(mp(il,i)-mp(il,i+1))
    26162680          up(il, i) = up(il, i)/mp(il, i)
    2617           vp(il, i) = vp(il, i+1)*mp(il, i+1) + v(il, i)*(mp(il,i)-mp(il,i+1) &
    2618             )
     2681          vp(il, i) = vp(il, i+1)*mp(il, i+1) + v(il, i)*(mp(il,i)-mp(il,i+1))
    26192682          vp(il, i) = vp(il, i)/mp(il, i)
    26202683
    … …  
    26232686          IF (mp(il,i+1)>1.0E-16) THEN
    26242687            IF (cvflag_grav) THEN
    2625               rp(il, i) = rp(il, i+1) + 100.*ginv*0.5*sigd(il)*(ph(il,i)-ph( &
    2626                 il,i+1))*(evap(il,i+1)+evap(il,i))/mp(il, i+1)
     2688              rp(il, i) = rp(il,i+1) + 100.*ginv*0.5*sigd(il)*(ph(il,i)-ph(il,i+1)) * &
     2689                                       (evap(il,i+1)+evap(il,i))/mp(il,i+1)
    26272690            ELSE
    2628               rp(il, i) = rp(il, i+1) + 5.*sigd(il)*(ph(il,i)-ph(il,i+1))*( &
    2629                 evap(il,i+1)+evap(il,i))/mp(il, i+1)
     2691              rp(il, i) = rp(il,i+1) + 5.*sigd(il)*(ph(il,i)-ph(il,i+1)) * &
     2692                                       (evap(il,i+1)+evap(il,i))/mp(il, i+1)
    26302693            END IF
    26312694            up(il, i) = up(il, i+1)
    … …  
    26392702      END IF ! (i.lt.inb(il) .and. lwork(il))
    26402703    END DO
    2641     ! ----------------------------------------------------------------
    2642 
    2643     ! ***       find tracer concentrations in precipitating downdraft     ***
    2644 
    2645     ! AC!      do j=1,ntra
    2646     ! AC!       do il = 1,ncum
    2647     ! AC!       if (i.lt.inb(il) .and. lwork(il)) then
    2648     ! AC!c
    2649     ! AC!         if(mplus(il))then
    2650     ! AC!          trap(il,i,j)=trap(il,i+1,j)*mp(il,i+1)
    2651     ! AC!     :              +trap(il,i,j)*(mp(il,i)-mp(il,i+1))
    2652     ! AC!          trap(il,i,j)=trap(il,i,j)/mp(il,i)
    2653     ! AC!         else ! if (mplus(il))
    2654     ! AC!          if(mp(il,i+1).gt.1.0e-16)then
    2655     ! AC!           trap(il,i,j)=trap(il,i+1,j)
    2656     ! AC!          endif
    2657     ! AC!         endif ! (mplus(il)) else if (.not.mplus(il))
    2658     ! AC!c
    2659     ! AC!        endif ! (i.lt.inb(il) .and. lwork(il))
    2660     ! AC!       enddo
    2661     ! AC!      end do
     2704! ----------------------------------------------------------------
     2705
     2706! ***       find tracer concentrations in precipitating downdraft     ***
     2707
     2708!AC!      do j=1,ntra
     2709!AC!       do il = 1,ncum
     2710!AC!       if (i.lt.inb(il) .and. lwork(il)) then
     2711!AC!c
     2712!AC!         if(mplus(il))then
     2713!AC!          trap(il,i,j)=trap(il,i+1,j)*mp(il,i+1)
     2714!AC!     :              +trap(il,i,j)*(mp(il,i)-mp(il,i+1))
     2715!AC!          trap(il,i,j)=trap(il,i,j)/mp(il,i)
     2716!AC!         else ! if (mplus(il))
     2717!AC!          if(mp(il,i+1).gt.1.0e-16)then
     2718!AC!           trap(il,i,j)=trap(il,i+1,j)
     2719!AC!          endif
     2720!AC!         endif ! (mplus(il)) else if (.not.mplus(il))
     2721!AC!c
     2722!AC!        endif ! (i.lt.inb(il) .and. lwork(il))
     2723!AC!       enddo
     2724!AC!      end do
    26622725
    26632726400 END DO
    2664   ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    2665 
    2666   ! ***                    end of downdraft loop                    ***
    2667 
    2668   ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     2727! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     2728
     2729! ***                    end of downdraft loop                    ***
     2730
     2731! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    26692732
    26702733
    … …  
    26722735END SUBROUTINE cv3_unsat
    26732736
    2674 SUBROUTINE cv3_yield(nloc, ncum, nd, na, ntra, icb, inb, delt, t, rr, t_wake, &
    2675     rr_wake, s_wake, u, v, tra, gz, p, ph, h, hp, lv, lf, cpn, th, th_wake, &
    2676     ep, clw, m, tp, mp, rp, up, vp, trap, wt, water, ice, evap, fondue, faci, &
    2677     b, sigd, ment, qent, hent, iflag_mix, uent, vent, nent, elij, traent, &
    2678     sig, tv, tvp, wghti, iflag, precip, vprecip, ft, fr, fu, fv, ftra, cbmf, &
    2679     upwd, dnwd, dnwd0, ma, mip, tls, tps, qcondc, wd, ftd, fqd)
     2737SUBROUTINE cv3_yield(nloc, ncum, nd, na, ntra, ok_conserv_q, &
     2738                     icb, inb, delt, &
     2739                     t, rr, t_wake, rr_wake, s_wake, u, v, tra, &
     2740                     gz, p, ph, h, hp, lv, lf, cpn, th, th_wake, &
     2741                     ep, clw, m, tp, mp, rp, up, vp, trap, &
     2742                     wt, water, ice, evap, fondue, faci, b, sigd, &
     2743                     ment, qent, hent, iflag_mix, uent, vent, &
     2744                     nent, elij, traent, sig, &
     2745                     tv, tvp, wghti, &
     2746                     iflag, precip, Vprecip, ft, fr, fu, fv, ftra, &
     2747                     cbmf, upwd, dnwd, dnwd0, ma, mip, &
     2748                     tls, tps, qcondc, wd, &
     2749                     ftd, fqd)
    26802750
    26812751  IMPLICIT NONE
    … …  
    26862756  include "conema3.h"
    26872757
    2688   ! inputs:
    2689   ! print*,'cv3_yield apres include'
    2690   INTEGER iflag_mix
    2691   INTEGER ncum, nd, na, ntra, nloc
    2692   INTEGER icb(nloc), inb(nloc)
    2693   REAL delt
    2694   REAL t(nloc, nd), rr(nloc, nd), u(nloc, nd), v(nloc, nd)
    2695   REAL t_wake(nloc, nd), rr_wake(nloc, nd)
    2696   REAL s_wake(nloc)
    2697   REAL tra(nloc, nd, ntra), sig(nloc, nd)
    2698   REAL gz(nloc, na), ph(nloc, nd+1), h(nloc, na), hp(nloc, na)
    2699   REAL th(nloc, na), p(nloc, nd), tp(nloc, na)
    2700   REAL lv(nloc, na), cpn(nloc, na), ep(nloc, na), clw(nloc, na)
    2701   REAL lf(nloc, na)
    2702   REAL m(nloc, na), mp(nloc, na), rp(nloc, na), up(nloc, na)
    2703   REAL vp(nloc, na), wt(nloc, nd), trap(nloc, nd, ntra)
    2704   REAL water(nloc, na), evap(nloc, na), b(nloc, na), sigd(nloc)
    2705   REAL fondue(nloc, na), faci(nloc, na), ice(nloc, na)
    2706   REAL ment(nloc, na, na), qent(nloc, na, na), uent(nloc, na, na)
    2707   REAL hent(nloc, na, na)
    2708   ! IM bug   real vent(nloc,na,na), nent(nloc,na), elij(nloc,na,na)
    2709   REAL vent(nloc, na, na), elij(nloc, na, na)
    2710   INTEGER nent(nloc, nd)
    2711   REAL traent(nloc, na, na, ntra)
    2712   REAL tv(nloc, nd), tvp(nloc, nd), wghti(nloc, nd)
    2713   ! print*,'cv3_yield declarations 1'
    2714   ! input/output:
    2715   INTEGER iflag(nloc)
    2716 
    2717   ! outputs:
    2718   REAL precip(nloc)
    2719   REAL ft(nloc, nd), fr(nloc, nd), fu(nloc, nd), fv(nloc, nd)
    2720   REAL ftd(nloc, nd), fqd(nloc, nd)
    2721   REAL ftra(nloc, nd, ntra)
    2722   REAL upwd(nloc, nd), dnwd(nloc, nd), ma(nloc, nd)
    2723   REAL dnwd0(nloc, nd), mip(nloc, nd)
    2724   REAL vprecip(nloc, nd+1)
    2725   REAL tls(nloc, nd), tps(nloc, nd)
    2726   REAL qcondc(nloc, nd) ! cld
    2727   REAL wd(nloc) ! gust
    2728   REAL cbmf(nloc)
    2729   ! print*,'cv3_yield declarations 2'
    2730   ! local variables:
    2731   INTEGER i, k, il, n, j, num1
    2732   REAL rat, delti
    2733   REAL ax, bx, cx, dx, ex
    2734   REAL cpinv, rdcp, dpinv
    2735   REAL awat(nloc)
    2736   REAL lvcp(nloc, na), lfcp(nloc, na), mke(nloc, na)
    2737   REAL am(nloc), work(nloc), ad(nloc), amp1(nloc)
    2738   ! !!      real up1(nloc), dn1(nloc)
    2739   REAL up1(nloc, nd, nd), dn1(nloc, nd, nd)
    2740   REAL asum(nloc), bsum(nloc), csum(nloc), dsum(nloc)
    2741   REAL esum(nloc), fsum(nloc), gsum(nloc), hsum(nloc)
    2742   REAL th_wake(nloc, nd)
    2743   REAL alpha_qpos(nloc), alpha_qpos1(nloc)
    2744   REAL qcond(nloc, nd), nqcond(nloc, nd), wa(nloc, nd) ! cld
    2745   REAL siga(nloc, nd), sax(nloc, nd), mac(nloc, nd) ! cld
    2746 
    2747   ! print*,'cv3_yield declarations 3'
    2748   ! -------------------------------------------------------------
    2749 
    2750   ! initialization:
     2758!inputs:
     2759      INTEGER iflag_mix
     2760      INTEGER ncum, nd, na, ntra, nloc
     2761      LOGICAL ok_conserv_q
     2762      INTEGER icb(nloc), inb(nloc)
     2763      REAL delt
     2764      REAL t(nloc, nd), rr(nloc, nd), u(nloc, nd), v(nloc, nd)
     2765      REAL t_wake(nloc, nd), rr_wake(nloc, nd)
     2766      REAL s_wake(nloc)
     2767      REAL tra(nloc, nd, ntra), sig(nloc, nd)
     2768      REAL gz(nloc, na), ph(nloc, nd+1), h(nloc, na), hp(nloc, na)
     2769      REAL th(nloc, na), p(nloc, nd), tp(nloc, na)
     2770      REAL lv(nloc, na), cpn(nloc, na), ep(nloc, na), clw(nloc, na)
     2771      REAL lf(nloc, na)
     2772      REAL m(nloc, na), mp(nloc, na), rp(nloc, na), up(nloc, na)
     2773      REAL vp(nloc, na), wt(nloc, nd), trap(nloc, nd, ntra)
     2774      REAL water(nloc, na), evap(nloc, na), b(nloc, na), sigd(nloc)
     2775      REAL fondue(nloc, na), faci(nloc, na), ice(nloc, na)
     2776      REAL ment(nloc, na, na), qent(nloc, na, na), uent(nloc, na, na)
     2777      REAL hent(nloc, na, na)
     2778!IM bug   real vent(nloc,na,na), nent(nloc,na), elij(nloc,na,na)
     2779      REAL vent(nloc, na, na), elij(nloc, na, na)
     2780      INTEGER nent(nloc, nd)
     2781      REAL traent(nloc, na, na, ntra)
     2782      REAL tv(nloc, nd), tvp(nloc, nd), wghti(nloc, nd)
     2783!
     2784!input/output:
     2785      INTEGER iflag(nloc)
     2786!
     2787!outputs:
     2788      REAL precip(nloc)
     2789      REAL ft(nloc, nd), fr(nloc, nd), fu(nloc, nd), fv(nloc, nd)
     2790      REAL ftd(nloc, nd), fqd(nloc, nd)
     2791      REAL ftra(nloc, nd, ntra)
     2792      REAL upwd(nloc, nd), dnwd(nloc, nd), ma(nloc, nd)
     2793      REAL dnwd0(nloc, nd), mip(nloc, nd)
     2794      REAL Vprecip(nloc, nd+1)
     2795      REAL tls(nloc, nd), tps(nloc, nd)
     2796      REAL qcondc(nloc, nd) ! cld
     2797      REAL wd(nloc) ! gust
     2798      REAL cbmf(nloc)
     2799!
     2800!local variables:
     2801      INTEGER i, k, il, n, j, num1
     2802      REAL rat, delti
     2803      REAL ax, bx, cx, dx, ex
     2804      REAL cpinv, rdcp, dpinv
     2805      REAL awat(nloc)
     2806      REAL lvcp(nloc, na), lfcp(nloc, na), mke(nloc, na)
     2807      REAL am(nloc), work(nloc), ad(nloc), amp1(nloc)
     2808!!      real up1(nloc), dn1(nloc)
     2809      REAL up1(nloc, nd, nd), dn1(nloc, nd, nd)
     2810      REAL asum(nloc), bsum(nloc), csum(nloc), dsum(nloc)
     2811      REAL esum(nloc), fsum(nloc), gsum(nloc), hsum(nloc)
     2812      REAL th_wake(nloc, nd)
     2813      REAL alpha_qpos(nloc), alpha_qpos1(nloc)
     2814      REAL qcond(nloc, nd), nqcond(nloc, nd), wa(nloc, nd) ! cld
     2815      REAL siga(nloc, nd), sax(nloc, nd), mac(nloc, nd) ! cld
     2816
     2817      REAL sumdq !jyg
     2818!
     2819! -------------------------------------------------------------
     2820
     2821! initialization:
    27512822
    27522823  delti = 1.0/delt
    2753   ! print*,'cv3_yield initialisation delt', delt
    2754   ! precip,Vprecip,ft,fr,fu,fv,ftra
    2755   ! :                    ,cbmf,upwd,dnwd,dnwd0,ma,mip
    2756   ! :                    ,tls,tps,qcondc,wd
    2757   ! :                    ,ftd,fqd  )
     2824! print*,'cv3_yield initialisation delt', delt
     2825!
    27582826  DO il = 1, ncum
    27592827    precip(il) = 0.0
    2760     vprecip(il, nd+1) = 0.0
     2828    Vprecip(il, nd+1) = 0.0
    27612829    wd(il) = 0.0 ! gust
    27622830  END DO
    … …  
    27642832  DO i = 1, nd
    27652833    DO il = 1, ncum
    2766       vprecip(il, i) = 0.0
     2834      Vprecip(il, i) = 0.0
    27672835      ft(il, i) = 0.0
    27682836      fr(il, i) = 0.0
    … …  
    27802848    END DO
    27812849  END DO
    2782   ! print*,'cv3_yield initialisation 2'
    2783   ! AC!      do j=1,ntra
    2784   ! AC!       do i=1,nd
    2785   ! AC!        do il=1,ncum
    2786   ! AC!          ftra(il,i,j)=0.0
    2787   ! AC!        enddo
    2788   ! AC!       enddo
    2789   ! AC!      enddo
    2790   ! print*,'cv3_yield initialisation 3'
     2850! print*,'cv3_yield initialisation 2'
     2851!AC!      do j=1,ntra
     2852!AC!       do i=1,nd
     2853!AC!        do il=1,ncum
     2854!AC!          ftra(il,i,j)=0.0
     2855!AC!        enddo
     2856!AC!       enddo
     2857!AC!      enddo
     2858! print*,'cv3_yield initialisation 3'
    27912859  DO i = 1, nl
    27922860    DO il = 1, ncum
    … …  
    27982866
    27992867
    2800   ! ***  calculate surface precipitation in mm/day     ***
     2868! ***  calculate surface precipitation in mm/day     ***
    28012869
    28022870  DO il = 1, ncum
    28032871    IF (ep(il,inb(il))>=0.0001 .AND. iflag(il)<=1) THEN
    28042872      IF (cvflag_ice) THEN
    2805         IF (cvflag_grav) THEN
    2806           precip(il) = wt(il, 1)*sigd(il)*(water(il,1)+ice(il,1))*86400.* &
    2807             1000./(rowl*grav)
    2808         ELSE
    2809           precip(il) = wt(il, 1)*sigd(il)*(water(il,1)+ice(il,1))*8640.
    2810         END IF
     2873        precip(il) = wt(il, 1)*sigd(il)*(water(il,1)+ice(il,1)) &
     2874                              *86400.*1000./(rowl*grav)
    28112875      ELSE
    2812         IF (cvflag_grav) THEN
    2813           precip(il) = wt(il, 1)*sigd(il)*water(il, 1)*86400.*1000./ &
    2814             (rowl*grav)
    2815         ELSE
    2816           precip(il) = wt(il, 1)*sigd(il)*water(il, 1)*8640.
    2817         END IF
     2876        precip(il) = wt(il, 1)*sigd(il)*water(il, 1) &
     2877                              *86400.*1000./(rowl*grav)
    28182878      END IF
    28192879    END IF
    28202880  END DO
    2821   ! print*,'cv3_yield apres calcul precip'
    2822 
    2823 
    2824   ! ===  calculate vertical profile of  precipitation in kg/m2/s  ===
     2881! print*,'cv3_yield apres calcul precip'
     2882
     2883
     2884! ===  calculate vertical profile of  precipitation in kg/m2/s  ===
    28252885
    28262886  DO i = 1, nl
    … …  
    28282888      IF (ep(il,inb(il))>=0.0001 .AND. i<=inb(il) .AND. iflag(il)<=1) THEN
    28292889        IF (cvflag_ice) THEN
    2830           IF (cvflag_grav) THEN
    2831             vprecip(il, i) = wt(il, i)*sigd(il)*(water(il,i)+ice(il,i))/grav
    2832           ELSE
    2833             vprecip(il, i) = wt(il, i)*sigd(il)*(water(il,i)+ice(il,i))/10.
    2834           END IF
     2890          Vprecip(il, i) = wt(il, i)*sigd(il)*(water(il,i)+ice(il,i))/grav
    28352891        ELSE
    2836           IF (cvflag_grav) THEN
    2837             vprecip(il, i) = wt(il, i)*sigd(il)*water(il, i)/grav
    2838           ELSE
    2839             vprecip(il, i) = wt(il, i)*sigd(il)*water(il, i)/10.
    2840           END IF
     2892          Vprecip(il, i) = wt(il, i)*sigd(il)*water(il, i)/grav
    28412893        END IF
    28422894      END IF
    … …  
    28452897
    28462898
    2847   ! ***  Calculate downdraft velocity scale    ***
    2848   ! ***  NE PAS UTILISER POUR L'INSTANT ***
    2849 
    2850   ! !      do il=1,ncum
    2851   ! !        wd(il)=betad*abs(mp(il,icb(il)))*0.01*rrd*t(il,icb(il))
    2852   ! !     :                                  /(sigd(il)*p(il,icb(il)))
    2853   ! !      enddo
    2854 
    2855 
    2856   ! ***  calculate tendencies of lowest level potential temperature  ***
    2857   ! ***                      and mixing ratio                        ***
     2899! ***  Calculate downdraft velocity scale    ***
     2900! ***  NE PAS UTILISER POUR L'INSTANT ***
     2901
     2902!!      do il=1,ncum
     2903!!        wd(il)=betad*abs(mp(il,icb(il)))*0.01*rrd*t(il,icb(il)) &
     2904!!                                       /(sigd(il)*p(il,icb(il)))
     2905!!      enddo
     2906
     2907
     2908! ***  calculate tendencies of lowest level potential temperature  ***
     2909! ***                      and mixing ratio                        ***
    28582910
    28592911  DO il = 1, ncum
    … …  
    28702922  END DO
    28712923
    2872   ! print*,'cv3_yield avant ft'
    2873   ! AM is the part of cbmf taken from the first level
     2924!    print*,'cv3_yield avant ft'
     2925! am is the part of cbmf taken from the first level
    28742926  DO il = 1, ncum
    28752927    am(il) = cbmf(il)*wghti(il, 1)
    … …  
    28782930  DO il = 1, ncum
    28792931    IF (iflag(il)<=1) THEN
    2880       ! convect3      if((0.1*dpinv*am).ge.delti)iflag(il)=4
    2881       ! jyg  Correction pour conserver l'eau
    2882       ! cc       ft(il,1)=-0.5*lvcp(il,1)*sigd(il)*(evap(il,1)+evap(il,2))
    2883       ! !precip
     2932! convect3      if((0.1*dpinv*am).ge.delti)iflag(il)=4
     2933!JYG  Correction pour conserver l'eau
     2934! cc       ft(il,1)=-0.5*lvcp(il,1)*sigd(il)*(evap(il,1)+evap(il,2))          !precip
    28842935      IF (cvflag_ice) THEN
    28852936        ft(il, 1) = -lvcp(il, 1)*sigd(il)*evap(il, 1) - &
    2886           lfcp(il, 1)*sigd(il)*evap(il, 1)*faci(il, 1) - &
    2887           lfcp(il, 1)*sigd(il)*(fondue(il,1)*wt(il,1))/(100.*(ph(il,1)-ph(il, &
    2888           2))) !precip
     2937                     lfcp(il, 1)*sigd(il)*evap(il, 1)*faci(il, 1) - &
     2938                     lfcp(il, 1)*sigd(il)*(fondue(il,1)*wt(il,1)) / &
     2939                       (100.*(ph(il,1)-ph(il,2)))                            !precip
    28892940      ELSE
    28902941        ft(il, 1) = -lvcp(il, 1)*sigd(il)*evap(il, 1)
    28912942      END IF
    28922943
    2893       IF (cvflag_grav) THEN
    2894         ft(il, 1) = ft(il, 1) - 0.009*grav*sigd(il)*mp(il, 2)*t_wake(il, 1)*b &
    2895           (il, 1)*work(il)
     2944      ft(il, 1) = ft(il, 1) - 0.009*grav*sigd(il)*mp(il, 2)*t_wake(il, 1)*b(il, 1)*work(il)
     2945
     2946      IF (cvflag_ice) THEN
     2947        ft(il, 1) = ft(il, 1) + 0.01*sigd(il)*wt(il, 1)*(cl-cpd)*water(il, 2) * &
     2948                                     (t_wake(il,2)-t_wake(il,1))*work(il)/cpn(il, 1) + &
     2949                                0.01*sigd(il)*wt(il, 1)*(ci-cpd)*ice(il, 2) * &
     2950                                     (t_wake(il,2)-t_wake(il,1))*work(il)/cpn(il, 1)
    28962951      ELSE
    2897         ft(il, 1) = ft(il, 1) - 0.09*sigd(il)*mp(il, 2)*t_wake(il, 1)*b(il, 1 &
    2898           )*work(il)
     2952        ft(il, 1) = ft(il, 1) + 0.01*sigd(il)*wt(il, 1)*(cl-cpd)*water(il, 2) * &
     2953                                     (t_wake(il,2)-t_wake(il,1))*work(il)/cpn(il, 1)
    28992954      END IF
    29002955
    2901       IF (cvflag_ice) THEN
    2902         ft(il, 1) = ft(il, 1) + 0.01*sigd(il)*wt(il, 1)*(cl-cpd)*water(il, 2) &
    2903           *(t_wake(il,2)-t_wake(il,1))*work(il)/cpn(il, 1) + &
    2904           0.01*sigd(il)*wt(il, 1)*(ci-cpd)*ice(il, 2)* &
    2905           (t_wake(il,2)-t_wake(il,1))*work(il)/cpn(il, 1)
    2906       ELSE
    2907         ft(il, 1) = ft(il, 1) + 0.01*sigd(il)*wt(il, 1)*(cl-cpd)*water(il, 2) &
    2908           *(t_wake(il,2)-t_wake(il,1))*work(il)/cpn(il, 1)
    2909       END IF
    2910 
    2911       ftd(il, 1) = ft(il, 1) ! fin precip
    2912 
    2913       IF (cvflag_grav) THEN !sature
    2914         IF ((0.01*grav*work(il)*am(il))>=delti) iflag(il) = 1 !consist vect
    2915         ft(il, 1) = ft(il, 1) + 0.01*grav*work(il)*am(il)*(t(il,2)-t(il,1)+( &
    2916           gz(il,2)-gz(il,1))/cpn(il,1))
    2917       ELSE
    2918         IF ((0.1*work(il)*am(il))>=delti) iflag(il) = 1 !consistency vect
    2919         ft(il, 1) = ft(il, 1) + 0.1*work(il)*am(il)*(t(il,2)-t(il,1)+(gz(il, &
    2920           2)-gz(il,1))/cpn(il,1))
    2921       END IF
     2956      ftd(il, 1) = ft(il, 1)                                                  ! fin precip
     2957
     2958      IF ((0.01*grav*work(il)*am(il))>=delti) iflag(il) = 1 !consist vect
     2959      ft(il, 1) = ft(il, 1) + 0.01*grav*work(il)*am(il) * &
     2960                                   (t(il,2)-t(il,1)+(gz(il,2)-gz(il,1))/cpn(il,1))
    29222961    END IF ! iflag
    29232962  END DO
    … …  
    29272966    IF (iflag_mix>0) THEN
    29282967      DO il = 1, ncum
    2929         ! FH WARNING a modifier :
     2968! FH WARNING a modifier :
    29302969        cpinv = 0.
    2931         ! cpinv=1.0/cpn(il,1)
     2970! cpinv=1.0/cpn(il,1)
    29322971        IF (j<=inb(il) .AND. iflag(il)<=1) THEN
    2933           IF (cvflag_grav) THEN
    2934             ft(il, 1) = ft(il, 1) + 0.01*grav*work(il)*ment(il, j, 1)*(hent( &
    2935               il,j,1)-h(il,1)+t(il,1)*(cpv-cpd)*(rr(il,1)-qent(il,j, &
    2936               1)))*cpinv
    2937           ELSE
    2938             ft(il, 1) = ft(il, 1) + 0.1*work(il)*ment(il, j, 1)*(hent(il,j,1) &
    2939               -h(il,1)+t(il,1)*(cpv-cpd)*(rr(il,1)-qent(il,j,1)))*cpinv
    2940           END IF ! cvflag_grav
     2972          ft(il, 1) = ft(il, 1) + 0.01*grav*work(il)*ment(il, j, 1) * &
     2973                     (hent(il,j,1)-h(il,1)+t(il,1)*(cpv-cpd)*(rr(il,1)-qent(il,j,1)))*cpinv
    29412974        END IF ! j
    29422975      END DO
    29432976    END IF
    29442977  END DO
    2945     ! fin sature
     2978! fin sature
    29462979
    29472980
    29482981  DO il = 1, ncum
    29492982    IF (iflag(il)<=1) THEN
    2950       IF (cvflag_grav) THEN
    2951         ! jyg1  Correction pour mieux conserver l'eau (conformite avec
    2952         ! CONVECT4.3)
    2953         fr(il, 1) = 0.01*grav*mp(il, 2)*(rp(il,2)-rr_wake(il,1))*work(il) + &
    2954           sigd(il)*evap(il, 1)
    2955         ! cc     :          +sigd(il)*0.5*(evap(il,1)+evap(il,2))
    2956 
    2957         fqd(il, 1) = fr(il, 1) !precip
    2958 
    2959         fr(il, 1) = fr(il, 1) + 0.01*grav*am(il)*(rr(il,2)-rr(il,1))*work(il) !sature
    2960 
    2961         fu(il, 1) = fu(il, 1) + 0.01*grav*work(il)*(mp(il,2)*(up(il,2)-u(il, &
    2962           1))+am(il)*(u(il,2)-u(il,1)))
    2963         fv(il, 1) = fv(il, 1) + 0.01*grav*work(il)*(mp(il,2)*(vp(il,2)-v(il, &
    2964           1))+am(il)*(v(il,2)-v(il,1)))
    2965       ELSE ! cvflag_grav
    2966         fr(il, 1) = 0.1*mp(il, 2)*(rp(il,2)-rr_wake(il,1))*work(il) + &
    2967           sigd(il)*evap(il, 1)
    2968         ! cc     :          +sigd(il)*0.5*(evap(il,1)+evap(il,2))
    2969         fqd(il, 1) = fr(il, 1) !precip
    2970         fr(il, 1) = fr(il, 1) + 0.1*am(il)*(rr(il,2)-rr(il,1))*work(il)
    2971         fu(il, 1) = fu(il, 1) + 0.1*work(il)*(mp(il,2)*(up(il,2)-u(il, &
    2972           1))+am(il)*(u(il,2)-u(il,1)))
    2973         fv(il, 1) = fv(il, 1) + 0.1*work(il)*(mp(il,2)*(vp(il,2)-v(il, &
    2974           1))+am(il)*(v(il,2)-v(il,1)))
    2975       END IF ! cvflag_grav
     2983!JYG1  Correction pour mieux conserver l'eau (conformite avec CONVECT4.3)
     2984      fr(il, 1) = 0.01*grav*mp(il, 2)*(rp(il,2)-rr_wake(il,1))*work(il) + &
     2985                  sigd(il)*evap(il, 1)
     2986!!!                  sigd(il)*0.5*(evap(il,1)+evap(il,2))
     2987
     2988      fqd(il, 1) = fr(il, 1) !precip
     2989
     2990      fr(il, 1) = fr(il, 1) + 0.01*grav*am(il)*(rr(il,2)-rr(il,1))*work(il)        !sature
     2991
     2992      fu(il, 1) = fu(il, 1) + 0.01*grav*work(il)*(mp(il,2)*(up(il,2)-u(il,1)) + &
     2993                                                  am(il)*(u(il,2)-u(il,1)))
     2994      fv(il, 1) = fv(il, 1) + 0.01*grav*work(il)*(mp(il,2)*(vp(il,2)-v(il,1)) + &
     2995                                                  am(il)*(v(il,2)-v(il,1)))
    29762996    END IF ! iflag
    29772997  END DO ! il
    29782998
    29792999
    2980   ! AC!     do j=1,ntra
    2981   ! AC!      do il=1,ncum
    2982   ! AC!       if (iflag(il) .le. 1) then
    2983   ! AC!       if (cvflag_grav) then
    2984   ! AC!        ftra(il,1,j)=ftra(il,1,j)+0.01*grav*work(il)
    2985   ! AC!    :                     *(mp(il,2)*(trap(il,2,j)-tra(il,1,j))
    2986   ! AC!    :             +am(il)*(tra(il,2,j)-tra(il,1,j)))
    2987   ! AC!       else
    2988   ! AC!        ftra(il,1,j)=ftra(il,1,j)+0.1*work(il)
    2989   ! AC!    :                     *(mp(il,2)*(trap(il,2,j)-tra(il,1,j))
    2990   ! AC!    :             +am(il)*(tra(il,2,j)-tra(il,1,j)))
    2991   ! AC!       endif
    2992   ! AC!       endif  ! iflag
    2993   ! AC!      enddo
    2994   ! AC!     enddo
     3000!AC!     do j=1,ntra
     3001!AC!      do il=1,ncum
     3002!AC!       if (iflag(il) .le. 1) then
     3003!AC!       if (cvflag_grav) then
     3004!AC!        ftra(il,1,j)=ftra(il,1,j)+0.01*grav*work(il)
     3005!AC!    :                     *(mp(il,2)*(trap(il,2,j)-tra(il,1,j))
     3006!AC!    :             +am(il)*(tra(il,2,j)-tra(il,1,j)))
     3007!AC!       else
     3008!AC!        ftra(il,1,j)=ftra(il,1,j)+0.1*work(il)
     3009!AC!    :                     *(mp(il,2)*(trap(il,2,j)-tra(il,1,j))
     3010!AC!    :             +am(il)*(tra(il,2,j)-tra(il,1,j)))
     3011!AC!       endif
     3012!AC!       endif  ! iflag
     3013!AC!      enddo
     3014!AC!     enddo
    29953015
    29963016  DO j = 2, nl
    29973017    DO il = 1, ncum
    29983018      IF (j<=inb(il) .AND. iflag(il)<=1) THEN
    2999         IF (cvflag_grav) THEN
    3000           fr(il, 1) = fr(il, 1) + 0.01*grav*work(il)*ment(il, j, 1)*(qent(il, &
    3001             j,1)-rr(il,1))
    3002           fu(il, 1) = fu(il, 1) + 0.01*grav*work(il)*ment(il, j, 1)*(uent(il, &
    3003             j,1)-u(il,1))
    3004           fv(il, 1) = fv(il, 1) + 0.01*grav*work(il)*ment(il, j, 1)*(vent(il, &
    3005             j,1)-v(il,1))
    3006         ELSE ! cvflag_grav
    3007           fr(il, 1) = fr(il, 1) + 0.1*work(il)*ment(il, j, 1)*(qent(il,j,1)- &
    3008             rr(il,1))
    3009           fu(il, 1) = fu(il, 1) + 0.1*work(il)*ment(il, j, 1)*(uent(il,j,1)-u &
    3010             (il,1))
    3011           fv(il, 1) = fv(il, 1) + 0.1*work(il)*ment(il, j, 1)*(vent(il,j,1)-v &
    3012             (il,1)) ! fin sature
    3013         END IF ! cvflag_grav
     3019        fr(il, 1) = fr(il, 1) + 0.01*grav*work(il)*ment(il, j, 1)*(qent(il,j,1)-rr(il,1))
     3020        fu(il, 1) = fu(il, 1) + 0.01*grav*work(il)*ment(il, j, 1)*(uent(il,j,1)-u(il,1))
     3021        fv(il, 1) = fv(il, 1) + 0.01*grav*work(il)*ment(il, j, 1)*(vent(il,j,1)-v(il,1))
    30143022      END IF ! j
    30153023    END DO
    30163024  END DO
    30173025
    3018   ! AC!      do k=1,ntra
    3019   ! AC!       do j=2,nl
    3020   ! AC!        do il=1,ncum
    3021   ! AC!         if (j.le.inb(il) .and. iflag(il) .le. 1) then
    3022   ! AC!
    3023   ! AC!          if (cvflag_grav) then
    3024   ! AC!           ftra(il,1,k)=ftra(il,1,k)+0.01*grav*work(il)*ment(il,j,1)
    3025   ! AC!     :                *(traent(il,j,1,k)-tra(il,1,k))
    3026   ! AC!          else
    3027   ! AC!           ftra(il,1,k)=ftra(il,1,k)+0.1*work(il)*ment(il,j,1)
    3028   ! AC!     :                *(traent(il,j,1,k)-tra(il,1,k))
    3029   ! AC!          endif
    3030   ! AC!
    3031   ! AC!         endif
    3032   ! AC!        enddo
    3033   ! AC!       enddo
    3034   ! AC!      enddo
    3035   ! print*,'cv3_yield apres ft'
    3036 
    3037   ! ***  calculate tendencies of potential temperature and mixing ratio  ***
    3038   ! ***               at levels above the lowest level                   ***
    3039 
    3040   ! ***  first find the net saturated updraft and downdraft mass fluxes  ***
    3041   ! ***                      through each level                          ***
     3026!AC!      do k=1,ntra
     3027!AC!       do j=2,nl
     3028!AC!        do il=1,ncum
     3029!AC!         if (j.le.inb(il) .and. iflag(il) .le. 1) then
     3030!AC!
     3031!AC!          if (cvflag_grav) then
     3032!AC!           ftra(il,1,k)=ftra(il,1,k)+0.01*grav*work(il)*ment(il,j,1)
     3033!AC!     :                *(traent(il,j,1,k)-tra(il,1,k))
     3034!AC!          else
     3035!AC!           ftra(il,1,k)=ftra(il,1,k)+0.1*work(il)*ment(il,j,1)
     3036!AC!     :                *(traent(il,j,1,k)-tra(il,1,k))
     3037!AC!          endif
     3038!AC!
     3039!AC!         endif
     3040!AC!        enddo
     3041!AC!       enddo
     3042!AC!      enddo
     3043! print*,'cv3_yield apres ft'
     3044
     3045! ***  calculate tendencies of potential temperature and mixing ratio  ***
     3046! ***               at levels above the lowest level                   ***
     3047
     3048! ***  first find the net saturated updraft and downdraft mass fluxes  ***
     3049! ***                      through each level                          ***
    30423050
    30433051
    … …  
    30603068          END IF
    30613069        ELSE
    3062           ! AMP1 is the part of cbmf taken from layers I and lower
     3070! AMP1 is the part of cbmf taken from layers I and lower
    30633071          IF (k<=i) THEN
    30643072            amp1(il) = amp1(il) + cbmf(il)*wghti(il, k)
    … …  
    30933101        cpinv = 1.0/cpn(il, i)
    30943102
    3095         ! convect3      if((0.1*dpinv*amp1).ge.delti)iflag(il)=4
    3096         IF (cvflag_grav) THEN
    3097           IF ((0.01*grav*dpinv*amp1(il))>=delti) iflag(il) = 1 ! vecto
    3098         ELSE
    3099           IF ((0.1*dpinv*amp1(il))>=delti) iflag(il) = 1 ! vecto
    3100         END IF
    3101 
    3102           ! precip
    3103         ! cc       ft(il,i)=
    3104         ! -0.5*sigd(il)*lvcp(il,i)*(evap(il,i)+evap(il,i+1))
     3103! convect3      if((0.1*dpinv*amp1).ge.delti)iflag(il)=4
     3104        IF ((0.01*grav*dpinv*amp1(il))>=delti) iflag(il) = 1 ! vecto
     3105
     3106! precip
     3107! cc       ft(il,i)= -0.5*sigd(il)*lvcp(il,i)*(evap(il,i)+evap(il,i+1))
    31053108        IF (cvflag_ice) THEN
    31063109          ft(il, i) = -sigd(il)*lvcp(il, i)*evap(il, i) - &
    3107             sigd(il)*lfcp(il, i)*evap(il, i)*faci(il, i) - &
    3108             sigd(il)*lfcp(il, i)*fondue(il, i)*wt(il, i)/(100.*(p(il, &
    3109             i-1)-p(il,i)))
     3110                       sigd(il)*lfcp(il, i)*evap(il, i)*faci(il, i) - &
     3111                       sigd(il)*lfcp(il, i)*fondue(il, i)*wt(il, i)/(100.*(p(il,i-1)-p(il,i)))
    31103112        ELSE
    31113113          ft(il, i) = -sigd(il)*lvcp(il, i)*evap(il, i)
    … …  
    31143116        rat = cpn(il, i-1)*cpinv
    31153117
    3116         IF (cvflag_grav) THEN
    3117           ft(il, i) = ft(il, i) - 0.009*grav*sigd(il)*(mp(il,i+1)*t_wake(il,i &
    3118             )*b(il,i)-mp(il,i)*t_wake(il,i-1)*rat*b(il,i-1))*dpinv
    3119           IF (cvflag_ice) THEN
    3120             ft(il, i) = ft(il, i) + 0.01*sigd(il)*wt(il, i)*(cl-cpd)*water(il &
    3121               , i+1)*(t_wake(il,i+1)-t_wake(il,i))*dpinv*cpinv + &
    3122               0.01*sigd(il)*wt(il, i)*(ci-cpd)*ice(il, i+1)* &
    3123               (t_wake(il,i+1)-t_wake(il,i))*dpinv*cpinv
    3124           ELSE
    3125             ft(il, i) = ft(il, i) + 0.01*sigd(il)*wt(il, i)*(cl-cpd)*water(il &
    3126               , i+1)*(t_wake(il,i+1)-t_wake(il,i))*dpinv*cpinv
    3127           END IF
    3128 
    3129           ftd(il, i) = ft(il, i)
    3130             ! fin precip
    3131 
    3132             ! sature
    3133           ft(il, i) = ft(il, i) + 0.01*grav*dpinv*(amp1(il)*(t(il,i+1)-t(il, &
    3134             i)+(gz(il,i+1)-gz(il,i))*cpinv)-ad(il)*(t(il,i)-t(il, &
    3135             i-1)+(gz(il,i)-gz(il,i-1))*cpinv))
    3136 
    3137 
    3138           IF (iflag_mix==0) THEN
    3139             ft(il, i) = ft(il, i) + 0.01*grav*dpinv*ment(il, i, i)*(hp(il,i)- &
    3140               h(il,i)+t(il,i)*(cpv-cpd)*(rr(il,i)-qent(il,i,i)))*cpinv
    3141           END IF
    3142 
    3143         ELSE ! cvflag_grav
    3144           ft(il, i) = ft(il, i) - 0.09*sigd(il)*(mp(il,i+1)*t_wake(il,i)*b(il &
    3145             ,i)-mp(il,i)*t_wake(il,i-1)*rat*b(il,i-1))*dpinv
    3146 
    3147           IF (cvflag_ice) THEN
    3148             ft(il, i) = ft(il, i) + 0.01*sigd(il)*wt(il, i)*(cl-cpd)*water(il &
    3149               , i+1)*(t_wake(il,i+1)-t_wake(il,i))*dpinv*cpinv + &
    3150               0.01*sigd(il)*wt(il, i)*(ci-cpd)*ice(il, i+1)* &
    3151               (t_wake(il,i+1)-t_wake(il,i))*dpinv*cpinv
    3152           ELSE
    3153             ft(il, i) = ft(il, i) + 0.01*sigd(il)*wt(il, i)*(cl-cpd)*water(il &
    3154               , i+1)*(t_wake(il,i+1)-t_wake(il,i))*dpinv*cpinv
    3155           END IF
    3156 
    3157           ftd(il, i) = ft(il, i)
    3158             ! fin precip
    3159 
    3160             ! sature
    3161           ft(il, i) = ft(il, i) + 0.1*dpinv*(amp1(il)*(t(il,i+1)-t(il, &
    3162             i)+(gz(il,i+1)-gz(il,i))*cpinv)-ad(il)*(t(il,i)-t(il, &
    3163             i-1)+(gz(il,i)-gz(il,i-1))*cpinv))
    3164 
    3165 
    3166           IF (iflag_mix==0) THEN
    3167             ft(il, i) = ft(il, i) + 0.1*dpinv*ment(il, i, i)*(hp(il,i)-h(il,i &
    3168               )+t(il,i)*(cpv-cpd)*(rr(il,i)-qent(il,i,i)))*cpinv
    3169           END IF
    3170         END IF ! cvflag_grav
    3171 
    3172 
    3173         IF (cvflag_grav) THEN
    3174           ! sb: on ne fait pas encore la correction permettant de mieux
    3175           ! conserver l'eau:
    3176           ! jyg: correction permettant de mieux conserver l'eau:
    3177           ! cc         fr(il,i)=0.5*sigd(il)*(evap(il,i)+evap(il,i+1))
    3178           fr(il, i) = sigd(il)*evap(il, i) + 0.01*grav*(mp(il,i+1)*(rp(il, &
    3179             i+1)-rr_wake(il,i))-mp(il,i)*(rp(il,i)-rr_wake(il,i-1)))*dpinv
    3180           fqd(il, i) = fr(il, i) ! precip
    3181 
    3182           fu(il, i) = 0.01*grav*(mp(il,i+1)*(up(il,i+1)-u(il, &
    3183             i))-mp(il,i)*(up(il,i)-u(il,i-1)))*dpinv
    3184           fv(il, i) = 0.01*grav*(mp(il,i+1)*(vp(il,i+1)-v(il, &
    3185             i))-mp(il,i)*(vp(il,i)-v(il,i-1)))*dpinv
    3186         ELSE ! cvflag_grav
    3187           ! cc         fr(il,i)=0.5*sigd(il)*(evap(il,i)+evap(il,i+1))
    3188           fr(il, i) = sigd(il)*evap(il, i) + 0.1*(mp(il,i+1)*(rp(il, &
    3189             i+1)-rr_wake(il,i))-mp(il,i)*(rp(il,i)-rr_wake(il,i-1)))*dpinv
    3190           fqd(il, i) = fr(il, i) ! precip
    3191 
    3192           fu(il, i) = 0.1*(mp(il,i+1)*(up(il,i+1)-u(il,i))-mp(il,i)*(up(il, &
    3193             i)-u(il,i-1)))*dpinv
    3194           fv(il, i) = 0.1*(mp(il,i+1)*(vp(il,i+1)-v(il,i))-mp(il,i)*(vp(il, &
    3195             i)-v(il,i-1)))*dpinv
    3196         END IF ! cvflag_grav
    3197 
    3198 
    3199         IF (cvflag_grav) THEN
    3200           fr(il, i) = fr(il, i) + 0.01*grav*dpinv*(amp1(il)*(rr(il, &
    3201             i+1)-rr(il,i))-ad(il)*(rr(il,i)-rr(il,i-1)))
    3202           fu(il, i) = fu(il, i) + 0.01*grav*dpinv*(amp1(il)*(u(il,i+1)-u(il, &
    3203             i))-ad(il)*(u(il,i)-u(il,i-1)))
    3204           fv(il, i) = fv(il, i) + 0.01*grav*dpinv*(amp1(il)*(v(il,i+1)-v(il, &
    3205             i))-ad(il)*(v(il,i)-v(il,i-1)))
    3206         ELSE ! cvflag_grav
    3207           fr(il, i) = fr(il, i) + 0.1*dpinv*(amp1(il)*(rr(il,i+1)-rr(il, &
    3208             i))-ad(il)*(rr(il,i)-rr(il,i-1)))
    3209           fu(il, i) = fu(il, i) + 0.1*dpinv*(amp1(il)*(u(il,i+1)-u(il, &
    3210             i))-ad(il)*(u(il,i)-u(il,i-1)))
    3211           fv(il, i) = fv(il, i) + 0.1*dpinv*(amp1(il)*(v(il,i+1)-v(il, &
    3212             i))-ad(il)*(v(il,i)-v(il,i-1)))
    3213         END IF ! cvflag_grav
     3118        ft(il, i) = ft(il, i) - 0.009*grav*sigd(il) * &
     3119                     (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
     3120        IF (cvflag_ice) THEN
     3121          ft(il, i) = ft(il, i) + 0.01*sigd(il)*wt(il, i)*(cl-cpd)*water(il, i+1) * &
     3122                                       (t_wake(il,i+1)-t_wake(il,i))*dpinv*cpinv + &
     3123                                  0.01*sigd(il)*wt(il, i)*(ci-cpd)*ice(il, i+1) * &
     3124                                       (t_wake(il,i+1)-t_wake(il,i))*dpinv*cpinv
     3125        ELSE
     3126          ft(il, i) = ft(il, i) + 0.01*sigd(il)*wt(il, i)*(cl-cpd)*water(il, i+1) * &
     3127                                       (t_wake(il,i+1)-t_wake(il,i))*dpinv* &
     3128            cpinv
     3129        END IF
     3130
     3131        ftd(il, i) = ft(il, i)
     3132! fin precip
     3133
     3134! sature
     3135        ft(il, i) = ft(il, i) + 0.01*grav*dpinv * &
     3136                     (amp1(il)*(t(il,i+1)-t(il,i) + (gz(il,i+1)-gz(il,i))*cpinv) - &
     3137                      ad(il)*(t(il,i)-t(il,i-1)+(gz(il,i)-gz(il,i-1))*cpinv))
     3138
     3139
     3140        IF (iflag_mix==0) THEN
     3141          ft(il, i) = ft(il, i) + 0.01*grav*dpinv*ment(il, i, i)*(hp(il,i)-h(il,i) + &
     3142                                    t(il,i)*(cpv-cpd)*(rr(il,i)-qent(il,i,i)))*cpinv
     3143        END IF
     3144
     3145
     3146
     3147! sb: on ne fait pas encore la correction permettant de mieux
     3148! conserver l'eau:
     3149!JYG: correction permettant de mieux conserver l'eau:
     3150! cc         fr(il,i)=0.5*sigd(il)*(evap(il,i)+evap(il,i+1))
     3151        fr(il, i) = sigd(il)*evap(il, i) + 0.01*grav*(mp(il,i+1)*(rp(il,i+1)-rr_wake(il,i)) - &
     3152                                                      mp(il,i)*(rp(il,i)-rr_wake(il,i-1)))*dpinv
     3153        fqd(il, i) = fr(il, i)                                                                     ! precip
     3154
     3155        fu(il, i) = 0.01*grav*(mp(il,i+1)*(up(il,i+1)-u(il,i)) - &
     3156                               mp(il,i)*(up(il,i)-u(il,i-1)))*dpinv
     3157        fv(il, i) = 0.01*grav*(mp(il,i+1)*(vp(il,i+1)-v(il,i)) - &
     3158                               mp(il,i)*(vp(il,i)-v(il,i-1)))*dpinv
     3159
     3160
     3161        fr(il, i) = fr(il, i) + 0.01*grav*dpinv*(amp1(il)*(rr(il,i+1)-rr(il,i)) - &
     3162                                                 ad(il)*(rr(il,i)-rr(il,i-1)))
     3163        fu(il, i) = fu(il, i) + 0.01*grav*dpinv*(amp1(il)*(u(il,i+1)-u(il,i)) - &
     3164                                                 ad(il)*(u(il,i)-u(il,i-1)))
     3165        fv(il, i) = fv(il, i) + 0.01*grav*dpinv*(amp1(il)*(v(il,i+1)-v(il,i)) - &
     3166                                                 ad(il)*(v(il,i)-v(il,i-1)))
    32143167
    32153168      END IF ! i
    32163169    END DO
    32173170
    3218     ! AC!      do k=1,ntra
    3219     ! AC!       do il=1,ncum
    3220     ! AC!        if (i.le.inb(il) .and. iflag(il) .le. 1) then
    3221     ! AC!         dpinv=1.0/(ph(il,i)-ph(il,i+1))
    3222     ! AC!         cpinv=1.0/cpn(il,i)
    3223     ! AC!         if (cvflag_grav) then
    3224     ! AC!           ftra(il,i,k)=ftra(il,i,k)+0.01*grav*dpinv
    3225     ! AC!     :         *(amp1(il)*(tra(il,i+1,k)-tra(il,i,k))
    3226     ! AC!     :           -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
    3227     ! AC!         else
    3228     ! AC!           ftra(il,i,k)=ftra(il,i,k)+0.1*dpinv
    3229     ! AC!     :         *(amp1(il)*(tra(il,i+1,k)-tra(il,i,k))
    3230     ! AC!     :           -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
    3231     ! AC!         endif
    3232     ! AC!        endif
    3233     ! AC!       enddo
    3234     ! AC!      enddo
     3171!AC!      do k=1,ntra
     3172!AC!       do il=1,ncum
     3173!AC!        if (i.le.inb(il) .and. iflag(il) .le. 1) then
     3174!AC!         dpinv=1.0/(ph(il,i)-ph(il,i+1))
     3175!AC!         cpinv=1.0/cpn(il,i)
     3176!AC!         if (cvflag_grav) then
     3177!AC!           ftra(il,i,k)=ftra(il,i,k)+0.01*grav*dpinv
     3178!AC!     :         *(amp1(il)*(tra(il,i+1,k)-tra(il,i,k))
     3179!AC!     :           -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
     3180!AC!         else
     3181!AC!           ftra(il,i,k)=ftra(il,i,k)+0.1*dpinv
     3182!AC!     :         *(amp1(il)*(tra(il,i+1,k)-tra(il,i,k))
     3183!AC!     :           -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
     3184!AC!         endif
     3185!AC!        endif
     3186!AC!       enddo
     3187!AC!      enddo
    32353188
    32363189    DO k = 1, i - 1
    … …  
    32463199            dpinv = 1.0/(ph(il,i)-ph(il,i+1))
    32473200            cpinv = 1.0/cpn(il, i)
    3248             IF (cvflag_grav) THEN
    3249               ft(il, i) = ft(il, i) + 0.01*grav*dpinv*ment(il, k, i)*(hent(il &
    3250                 ,k,i)-h(il,i)+t(il,i)*(cpv-cpd)*(rr(il,i)+awat(il)-qent(il,k, &
    3251                 i)))*cpinv
    3252 
    3253 
    3254 
    3255             ELSE
    3256               ft(il, i) = ft(il, i) + 0.1*dpinv*ment(il, k, i)*(hent(il,k,i)- &
    3257                 h(il,i)+t(il,i)*(cpv-cpd)*(rr(il,i)+awat(il)-qent(il,k, &
    3258                 i)))*cpinv
    3259             END IF !cvflag_grav
     3201            ft(il, i) = ft(il, i) + 0.01*grav*dpinv*ment(il, k, i) * &
     3202                 (hent(il,k,i)-h(il,i)+t(il,i)*(cpv-cpd)*(rr(il,i)+awat(il)-qent(il,k,i)))*cpinv
     3203!
     3204!
    32603205          END IF ! i
    32613206        END DO
    … …  
    32663211          dpinv = 1.0/(ph(il,i)-ph(il,i+1))
    32673212          cpinv = 1.0/cpn(il, i)
    3268           IF (cvflag_grav) THEN
    3269             fr(il, i) = fr(il, i) + 0.01*grav*dpinv*ment(il, k, i)*(qent(il,k &
    3270               ,i)-awat(il)-rr(il,i))
    3271             fu(il, i) = fu(il, i) + 0.01*grav*dpinv*ment(il, k, i)*(uent(il,k &
    3272               ,i)-u(il,i))
    3273             fv(il, i) = fv(il, i) + 0.01*grav*dpinv*ment(il, k, i)*(vent(il,k &
    3274               ,i)-v(il,i))
    3275           ELSE ! cvflag_grav
    3276             fr(il, i) = fr(il, i) + 0.1*dpinv*ment(il, k, i)*(qent(il,k,i)- &
    3277               awat(il)-rr(il,i))
    3278             fu(il, i) = fu(il, i) + 0.01*grav*dpinv*ment(il, k, i)*(uent(il,k &
    3279               ,i)-u(il,i))
    3280             fv(il, i) = fv(il, i) + 0.1*dpinv*ment(il, k, i)*(vent(il,k,i)-v( &
    3281               il,i))
    3282           END IF ! cvflag_grav
    3283 
    3284           ! (saturated updrafts resulting from mixing)        ! cld
    3285           qcond(il, i) = qcond(il, i) + (elij(il,k,i)-awat(il)) ! cld
     3213          fr(il, i) = fr(il, i) + 0.01*grav*dpinv*ment(il, k, i) * &
     3214                                                       (qent(il,k,i)-awat(il)-rr(il,i))
     3215          fu(il, i) = fu(il, i) + 0.01*grav*dpinv*ment(il, k, i)*(uent(il,k,i)-u(il,i))
     3216          fv(il, i) = fv(il, i) + 0.01*grav*dpinv*ment(il, k, i)*(vent(il,k,i)-v(il,i))
     3217
     3218! (saturated updrafts resulting from mixing)                                   ! cld
     3219          qcond(il, i) = qcond(il, i) + (elij(il,k,i)-awat(il))                ! cld
    32863220          nqcond(il, i) = nqcond(il, i) + 1. ! cld
    32873221        END IF ! i
    … …  
    32893223    END DO
    32903224
    3291     ! AC!      do j=1,ntra
    3292     ! AC!       do k=1,i-1
    3293     ! AC!        do il=1,ncum
    3294     ! AC!         if (i.le.inb(il) .and. iflag(il) .le. 1) then
    3295     ! AC!          dpinv=1.0/(ph(il,i)-ph(il,i+1))
    3296     ! AC!          cpinv=1.0/cpn(il,i)
    3297     ! AC!          if (cvflag_grav) then
    3298     ! AC!           ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv*ment(il,k,i)
    3299     ! AC!     :        *(traent(il,k,i,j)-tra(il,i,j))
    3300     ! AC!          else
    3301     ! AC!           ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv*ment(il,k,i)
    3302     ! AC!     :        *(traent(il,k,i,j)-tra(il,i,j))
    3303     ! AC!          endif
    3304     ! AC!         endif
    3305     ! AC!        enddo
    3306     ! AC!       enddo
    3307     ! AC!      enddo
     3225!AC!      do j=1,ntra
     3226!AC!       do k=1,i-1
     3227!AC!        do il=1,ncum
     3228!AC!         if (i.le.inb(il) .and. iflag(il) .le. 1) then
     3229!AC!          dpinv=1.0/(ph(il,i)-ph(il,i+1))
     3230!AC!          cpinv=1.0/cpn(il,i)
     3231!AC!          if (cvflag_grav) then
     3232!AC!           ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv*ment(il,k,i)
     3233!AC!     :        *(traent(il,k,i,j)-tra(il,i,j))
     3234!AC!          else
     3235!AC!           ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv*ment(il,k,i)
     3236!AC!     :        *(traent(il,k,i,j)-tra(il,i,j))
     3237!AC!          endif
     3238!AC!         endif
     3239!AC!        enddo
     3240!AC!       enddo
     3241!AC!      enddo
    33083242
    33093243    DO k = i, nl + 1
    … …  
    33143248            dpinv = 1.0/(ph(il,i)-ph(il,i+1))
    33153249            cpinv = 1.0/cpn(il, i)
    3316             IF (cvflag_grav) THEN
    3317               ft(il, i) = ft(il, i) + 0.01*grav*dpinv*ment(il, k, i)*(hent(il &
    3318                 ,k,i)-h(il,i)+t(il,i)*(cpv-cpd)*(rr(il,i)-qent(il,k, &
    3319                 i)))*cpinv
    3320 
    3321 
    3322             ELSE
    3323               ft(il, i) = ft(il, i) + 0.1*dpinv*ment(il, k, i)*(hent(il,k,i)- &
    3324                 h(il,i)+t(il,i)*(cpv-cpd)*(rr(il,i)-qent(il,k,i)))*cpinv
    3325             END IF !cvflag_grav
     3250            ft(il, i) = ft(il, i) + 0.01*grav*dpinv*ment(il, k, i) * &
     3251                  (hent(il,k,i)-h(il,i)+t(il,i)*(cpv-cpd)*(rr(il,i)-qent(il,k,i)))*cpinv
     3252
     3253
    33263254          END IF ! i
    33273255        END DO
    … …  
    33333261          cpinv = 1.0/cpn(il, i)
    33343262
    3335           IF (cvflag_grav) THEN
    3336             fr(il, i) = fr(il, i) + 0.01*grav*dpinv*ment(il, k, i)*(qent(il,k &
    3337               ,i)-rr(il,i))
    3338             fu(il, i) = fu(il, i) + 0.01*grav*dpinv*ment(il, k, i)*(uent(il,k &
    3339               ,i)-u(il,i))
    3340             fv(il, i) = fv(il, i) + 0.01*grav*dpinv*ment(il, k, i)*(vent(il,k &
    3341               ,i)-v(il,i))
    3342           ELSE ! cvflag_grav
    3343             fr(il, i) = fr(il, i) + 0.1*dpinv*ment(il, k, i)*(qent(il,k,i)-rr &
    3344               (il,i))
    3345             fu(il, i) = fu(il, i) + 0.1*dpinv*ment(il, k, i)*(uent(il,k,i)-u( &
    3346               il,i))
    3347             fv(il, i) = fv(il, i) + 0.1*dpinv*ment(il, k, i)*(vent(il,k,i)-v( &
    3348               il,i))
    3349           END IF ! cvflag_grav
     3263          fr(il, i) = fr(il, i) + 0.01*grav*dpinv*ment(il, k, i)*(qent(il,k,i)-rr(il,i))
     3264          fu(il, i) = fu(il, i) + 0.01*grav*dpinv*ment(il, k, i)*(uent(il,k,i)-u(il,i))
     3265          fv(il, i) = fv(il, i) + 0.01*grav*dpinv*ment(il, k, i)*(vent(il,k,i)-v(il,i))
    33503266        END IF ! i and k
    33513267      END DO
    33523268    END DO
    33533269
    3354     ! AC!      do j=1,ntra
    3355     ! AC!       do k=i,nl+1
    3356     ! AC!        do il=1,ncum
    3357     ! AC!         if (i.le.inb(il) .and. k.le.inb(il)
    3358     ! AC!     $                .and. iflag(il) .le. 1) then
    3359     ! AC!          dpinv=1.0/(ph(il,i)-ph(il,i+1))
    3360     ! AC!          cpinv=1.0/cpn(il,i)
    3361     ! AC!          if (cvflag_grav) then
    3362     ! AC!           ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv*ment(il,k,i)
    3363     ! AC!     :         *(traent(il,k,i,j)-tra(il,i,j))
    3364     ! AC!          else
    3365     ! AC!           ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv*ment(il,k,i)
    3366     ! AC!     :             *(traent(il,k,i,j)-tra(il,i,j))
    3367     ! AC!          endif
    3368     ! AC!         endif ! i and k
    3369     ! AC!        enddo
    3370     ! AC!       enddo
    3371     ! AC!      enddo
    3372 
    3373     ! sb: interface with the cloud parameterization:          ! cld
     3270!AC!      do j=1,ntra
     3271!AC!       do k=i,nl+1
     3272!AC!        do il=1,ncum
     3273!AC!         if (i.le.inb(il) .and. k.le.inb(il)
     3274!AC!     $                .and. iflag(il) .le. 1) then
     3275!AC!          dpinv=1.0/(ph(il,i)-ph(il,i+1))
     3276!AC!          cpinv=1.0/cpn(il,i)
     3277!AC!          if (cvflag_grav) then
     3278!AC!           ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv*ment(il,k,i)
     3279!AC!     :         *(traent(il,k,i,j)-tra(il,i,j))
     3280!AC!          else
     3281!AC!           ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv*ment(il,k,i)
     3282!AC!     :             *(traent(il,k,i,j)-tra(il,i,j))
     3283!AC!          endif
     3284!AC!         endif ! i and k
     3285!AC!        enddo
     3286!AC!       enddo
     3287!AC!      enddo
     3288
     3289! sb: interface with the cloud parameterization:                               ! cld
    33743290
    33753291    DO k = i + 1, nl
    33763292      DO il = 1, ncum
    3377         IF (k<=inb(il) .AND. i<=inb(il) .AND. iflag(il)<=1) THEN ! cld
    3378           ! (saturated downdrafts resulting from mixing)            ! cld
    3379           qcond(il, i) = qcond(il, i) + elij(il, k, i) ! cld
     3293        IF (k<=inb(il) .AND. i<=inb(il) .AND. iflag(il)<=1) THEN               ! cld
     3294! (saturated downdrafts resulting from mixing)                                 ! cld
     3295          qcond(il, i) = qcond(il, i) + elij(il, k, i)                         ! cld
    33803296          nqcond(il, i) = nqcond(il, i) + 1. ! cld
    33813297        END IF ! cld
    … …  
    33833299    END DO ! cld
    33843300
    3385     ! (particular case: no detraining level is found)         ! cld
    3386     DO il = 1, ncum ! cld
    3387       IF (i<=inb(il) .AND. nent(il,i)==0 .AND. iflag(il)<=1) THEN ! cld
    3388         qcond(il, i) = qcond(il, i) + (1.-ep(il,i))*clw(il, i) ! cld
    3389         nqcond(il, i) = nqcond(il, i) + 1. ! cld
    3390       END IF ! cld
    3391     END DO ! cld
    3392 
    3393     DO il = 1, ncum ! cld
    3394       IF (i<=inb(il) .AND. nqcond(il,i)/=0 .AND. iflag(il)<=1) THEN ! cld
    3395         qcond(il, i) = qcond(il, i)/nqcond(il, i) ! cld
    3396       END IF ! cld
    3397     END DO
    3398 
    3399     ! AC!      do j=1,ntra
    3400     ! AC!       do il=1,ncum
    3401     ! AC!        if (i.le.inb(il) .and. iflag(il) .le. 1) then
    3402     ! AC!         dpinv=1.0/(ph(il,i)-ph(il,i+1))
    3403     ! AC!         cpinv=1.0/cpn(il,i)
    3404     ! AC!
    3405     ! AC!         if (cvflag_grav) then
    3406     ! AC!          ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv
    3407     ! AC!     :     *(mp(il,i+1)*(trap(il,i+1,j)-tra(il,i,j))
    3408     ! AC!     :     -mp(il,i)*(trap(il,i,j)-trap(il,i-1,j)))
    3409     ! AC!         else
    3410     ! AC!          ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv
    3411     ! AC!     :     *(mp(il,i+1)*(trap(il,i+1,j)-tra(il,i,j))
    3412     ! AC!     :     -mp(il,i)*(trap(il,i,j)-trap(il,i-1,j)))
    3413     ! AC!         endif
    3414     ! AC!        endif ! i
    3415     ! AC!       enddo
    3416     ! AC!      enddo
     3301! (particular case: no detraining level is found)                              ! cld
     3302    DO il = 1, ncum                                                            ! cld
     3303      IF (i<=inb(il) .AND. nent(il,i)==0 .AND. iflag(il)<=1) THEN              ! cld
     3304        qcond(il, i) = qcond(il, i) + (1.-ep(il,i))*clw(il, i)                 ! cld
     3305        nqcond(il, i) = nqcond(il, i) + 1.                                     ! cld
     3306      END IF                                                                   ! cld
     3307    END DO                                                                     ! cld
     3308
     3309    DO il = 1, ncum                                                            ! cld
     3310      IF (i<=inb(il) .AND. nqcond(il,i)/=0 .AND. iflag(il)<=1) THEN            ! cld
     3311        qcond(il, i) = qcond(il, i)/nqcond(il, i)                              ! cld
     3312      END IF                                                                   ! cld
     3313    END DO
     3314
     3315!AC!      do j=1,ntra
     3316!AC!       do il=1,ncum
     3317!AC!        if (i.le.inb(il) .and. iflag(il) .le. 1) then
     3318!AC!         dpinv=1.0/(ph(il,i)-ph(il,i+1))
     3319!AC!         cpinv=1.0/cpn(il,i)
     3320!AC!
     3321!AC!         if (cvflag_grav) then
     3322!AC!          ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv
     3323!AC!     :     *(mp(il,i+1)*(trap(il,i+1,j)-tra(il,i,j))
     3324!AC!     :     -mp(il,i)*(trap(il,i,j)-trap(il,i-1,j)))
     3325!AC!         else
     3326!AC!          ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv
     3327!AC!     :     *(mp(il,i+1)*(trap(il,i+1,j)-tra(il,i,j))
     3328!AC!     :     -mp(il,i)*(trap(il,i,j)-trap(il,i-1,j)))
     3329!AC!         endif
     3330!AC!        endif ! i
     3331!AC!       enddo
     3332!AC!      enddo
    34173333
    34183334
    34193335500 END DO
    34203336
    3421 
    3422   ! ***   move the detrainment at level inb down to level inb-1   ***
    3423   ! ***        in such a way as to preserve the vertically        ***
    3424   ! ***          integrated enthalpy and water tendencies         ***
    3425 
    3426   ! Correction bug le 18-03-09
     3337!JYG<
     3338!Conservation de l'eau
     3339!   sumdq = 0.
     3340!   DO k = 1, nl
     3341!     sumdq = sumdq + fr(1, k)*100.*(ph(1,k)-ph(1,k+1))/grav
     3342!   END DO
     3343!   PRINT *, 'cv3_yield, apres 500, sum(dq), precip, somme ', sumdq, Vprecip(1, 1), sumdq + vprecip(1, 1)
     3344!JYG>
     3345! ***   move the detrainment at level inb down to level inb-1   ***
     3346! ***        in such a way as to preserve the vertically        ***
     3347! ***          integrated enthalpy and water tendencies         ***
     3348
     3349! Correction bug le 18-03-09
    34273350  DO il = 1, ncum
    34283351    IF (iflag(il)<=1) THEN
    3429       IF (cvflag_grav) THEN
    3430         ax = 0.01*grav*ment(il, inb(il), inb(il))*(hp(il,inb(il))-h(il,inb(il &
    3431           ))+t(il,inb(il))*(cpv-cpd)*(rr(il,inb(il))-qent(il,inb(il), &
    3432           inb(il))))/(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, &
    3435           inb(il))-ph(il,inb(il)+1))/(cpn(il,inb(il)-1)*(ph(il, &
    3436           inb(il)-1)-ph(il,inb(il))))
    3437 
    3438         bx = 0.01*grav*ment(il, inb(il), inb(il))*(qent(il,inb(il),inb(il))- &
    3439           rr(il,inb(il)))/(ph(il,inb(il))-ph(il,inb(il)+1))
    3440         fr(il, inb(il)) = fr(il, inb(il)) - bx
    3441         fr(il, inb(il)-1) = fr(il, inb(il)-1) + bx*(ph(il,inb(il))-ph(il,inb( &
    3442           il)+1))/(ph(il,inb(il)-1)-ph(il,inb(il)))
    3443 
    3444         cx = 0.01*grav*ment(il, inb(il), inb(il))*(uent(il,inb(il),inb(il))-u &
    3445           (il,inb(il)))/(ph(il,inb(il))-ph(il,inb(il)+1))
    3446         fu(il, inb(il)) = fu(il, inb(il)) - cx
    3447         fu(il, inb(il)-1) = fu(il, inb(il)-1) + cx*(ph(il,inb(il))-ph(il,inb( &
    3448           il)+1))/(ph(il,inb(il)-1)-ph(il,inb(il)))
    3449 
    3450         dx = 0.01*grav*ment(il, inb(il), inb(il))*(vent(il,inb(il),inb(il))-v &
    3451           (il,inb(il)))/(ph(il,inb(il))-ph(il,inb(il)+1))
    3452         fv(il, inb(il)) = fv(il, inb(il)) - dx
    3453         fv(il, inb(il)-1) = fv(il, inb(il)-1) + dx*(ph(il,inb(il))-ph(il,inb( &
    3454           il)+1))/(ph(il,inb(il)-1)-ph(il,inb(il)))
    3455       ELSE
    3456         ax = 0.1*ment(il, inb(il), inb(il))*(hp(il,inb(il))-h(il,inb(il))+t( &
    3457           il,inb(il))*(cpv-cpd)*(rr(il,inb(il))-qent(il,inb(il), &
    3458           inb(il))))/(cpn(il,inb(il))*(ph(il,inb(il))-ph(il,inb(il)+1)))
    3459         ft(il, inb(il)) = ft(il, inb(il)) - ax
    3460         ft(il, inb(il)-1) = ft(il, inb(il)-1) + ax*cpn(il, inb(il))*(ph(il, &
    3461           inb(il))-ph(il,inb(il)+1))/(cpn(il,inb(il)-1)*(ph(il, &
    3462           inb(il)-1)-ph(il,inb(il))))
    3463 
    3464         bx = 0.1*ment(il, inb(il), inb(il))*(qent(il,inb(il),inb(il))-rr(il, &
    3465           inb(il)))/(ph(il,inb(il))-ph(il,inb(il)+1))
    3466         fr(il, inb(il)) = fr(il, inb(il)) - bx
    3467         fr(il, inb(il)-1) = fr(il, inb(il)-1) + bx*(ph(il,inb(il))-ph(il,inb( &
    3468           il)+1))/(ph(il,inb(il)-1)-ph(il,inb(il)))
    3469 
    3470         cx = 0.1*ment(il, inb(il), inb(il))*(uent(il,inb(il),inb(il))-u(il, &
    3471           inb(il)))/(ph(il,inb(il))-ph(il,inb(il)+1))
    3472         fu(il, inb(il)) = fu(il, inb(il)) - cx
    3473         fu(il, inb(il)-1) = fu(il, inb(il)-1) + cx*(ph(il,inb(il))-ph(il,inb( &
    3474           il)+1))/(ph(il,inb(il)-1)-ph(il,inb(il)))
    3475 
    3476         dx = 0.1*ment(il, inb(il), inb(il))*(vent(il,inb(il),inb(il))-v(il, &
    3477           inb(il)))/(ph(il,inb(il))-ph(il,inb(il)+1))
    3478         fv(il, inb(il)) = fv(il, inb(il)) - dx
    3479         fv(il, inb(il)-1) = fv(il, inb(il)-1) + dx*(ph(il,inb(il))-ph(il,inb( &
    3480           il)+1))/(ph(il,inb(il)-1)-ph(il,inb(il)))
    3481       END IF
     3352      ax = 0.01*grav*ment(il, inb(il), inb(il))* &
     3353           (hp(il,inb(il))-h(il,inb(il))+t(il,inb(il))*(cpv-cpd)*(rr(il,inb(il))-qent(il,inb(il),inb(il))))/ &
     3354                                (cpn(il,inb(il))*(ph(il,inb(il))-ph(il,inb(il)+1)))
     3355      ft(il, inb(il)) = ft(il, inb(il)) - ax
     3356      ft(il, inb(il)-1) = ft(il, inb(il)-1) + ax*cpn(il, inb(il))*(ph(il,inb(il))-ph(il,inb(il)+1))/ &
     3357                              (cpn(il,inb(il)-1)*(ph(il,inb(il)-1)-ph(il,inb(il))))
     3358
     3359      bx = 0.01*grav*ment(il, inb(il), inb(il))*(qent(il,inb(il),inb(il))-rr(il,inb(il)))/ &
     3360                                                 (ph(il,inb(il))-ph(il,inb(il)+1))
     3361      fr(il, inb(il)) = fr(il, inb(il)) - bx
     3362      fr(il, inb(il)-1) = fr(il, inb(il)-1) + bx*(ph(il,inb(il))-ph(il,inb(il)+1))/ &
     3363                                                 (ph(il,inb(il)-1)-ph(il,inb(il)))
     3364
     3365      cx = 0.01*grav*ment(il, inb(il), inb(il))*(uent(il,inb(il),inb(il))-u(il,inb(il)))/ &
     3366                                                 (ph(il,inb(il))-ph(il,inb(il)+1))
     3367      fu(il, inb(il)) = fu(il, inb(il)) - cx
     3368      fu(il, inb(il)-1) = fu(il, inb(il)-1) + cx*(ph(il,inb(il))-ph(il,inb(il)+1))/ &
     3369                                                 (ph(il,inb(il)-1)-ph(il,inb(il)))
     3370
     3371      dx = 0.01*grav*ment(il, inb(il), inb(il))*(vent(il,inb(il),inb(il))-v(il,inb(il)))/ &
     3372                                                 (ph(il,inb(il))-ph(il,inb(il)+1))
     3373      fv(il, inb(il)) = fv(il, inb(il)) - dx
     3374      fv(il, inb(il)-1) = fv(il, inb(il)-1) + dx*(ph(il,inb(il))-ph(il,inb(il)+1))/ &
     3375                                                 (ph(il,inb(il)-1)-ph(il,inb(il)))
    34823376    END IF !iflag
    34833377  END DO
    34843378
    3485   ! AC!      do j=1,ntra
    3486   ! AC!       do il=1,ncum
    3487   ! AC!        IF (iflag(il) .le. 1) THEN
    3488   ! AC! IF (cvflag_grav) then
    3489   ! AC!        ex=0.01*grav*ment(il,inb(il),inb(il))
    3490   ! AC!     :      *(traent(il,inb(il),inb(il),j)-tra(il,inb(il),j))
    3491   ! AC!     :      /(ph(i l,inb(il))-ph(il,inb(il)+1))
    3492   ! AC!        ftra(il,inb(il),j)=ftra(il,inb(il),j)-ex
    3493   ! AC!        ftra(il,inb(il)-1,j)=ftra(il,inb(il)-1,j)
    3494   ! AC!     :       +ex*(ph(il,inb(il))-ph(il,inb(il)+1))
    3495   ! AC!     :          /(ph(il,inb(il)-1)-ph(il,inb(il)))
    3496   ! AC! else
    3497   ! AC!        ex=0.1*ment(il,inb(il),inb(il))
    3498   ! AC!     :      *(traent(il,inb(il),inb(il),j)-tra(il,inb(il),j))
    3499   ! AC!     :      /(ph(i l,inb(il))-ph(il,inb(il)+1))
    3500   ! AC!        ftra(il,inb(il),j)=ftra(il,inb(il),j)-ex
    3501   ! AC!        ftra(il,inb(il)-1,j)=ftra(il,inb(il)-1,j)
    3502   ! AC!     :       +ex*(ph(il,inb(il))-ph(il,inb(il)+1))
    3503   ! AC!     :          /(ph(il,inb(il)-1)-ph(il,inb(il)))
    3504   ! AC!        ENDIF   !cvflag grav
    3505   ! AC!        ENDIF    !iflag
    3506   ! AC!       enddo
    3507   ! AC!      enddo
    3508 
    3509 
    3510   ! ***    homogenize tendencies below cloud base    ***
     3379!JYG<
     3380!Conservation de l'eau
     3381!   sumdq = 0.
     3382!   DO k = 1, nl
     3383!     sumdq = sumdq + fr(1, k)*100.*(ph(1,k)-ph(1,k+1))/grav
     3384!   END DO
     3385!   PRINT *, 'cv3_yield, apres 503, sum(dq), precip, somme ', sumdq, Vprecip(1, 1), sumdq + vprecip(1, 1)
     3386!JYG>
     3387
     3388!AC!      do j=1,ntra
     3389!AC!       do il=1,ncum
     3390!AC!        IF (iflag(il) .le. 1) THEN
     3391!AC!    IF (cvflag_grav) then
     3392!AC!        ex=0.01*grav*ment(il,inb(il),inb(il))
     3393!AC!     :      *(traent(il,inb(il),inb(il),j)-tra(il,inb(il),j))
     3394!AC!     :      /(ph(i l,inb(il))-ph(il,inb(il)+1))
     3395!AC!        ftra(il,inb(il),j)=ftra(il,inb(il),j)-ex
     3396!AC!        ftra(il,inb(il)-1,j)=ftra(il,inb(il)-1,j)
     3397!AC!     :       +ex*(ph(il,inb(il))-ph(il,inb(il)+1))
     3398!AC!     :          /(ph(il,inb(il)-1)-ph(il,inb(il)))
     3399!AC!    else
     3400!AC!        ex=0.1*ment(il,inb(il),inb(il))
     3401!AC!     :      *(traent(il,inb(il),inb(il),j)-tra(il,inb(il),j))
     3402!AC!     :      /(ph(i l,inb(il))-ph(il,inb(il)+1))
     3403!AC!        ftra(il,inb(il),j)=ftra(il,inb(il),j)-ex
     3404!AC!        ftra(il,inb(il)-1,j)=ftra(il,inb(il)-1,j)
     3405!AC!     :       +ex*(ph(il,inb(il))-ph(il,inb(il)+1))
     3406!AC!     :          /(ph(il,inb(il)-1)-ph(il,inb(il)))
     3407!AC!        ENDIF   !cvflag grav
     3408!AC!        ENDIF    !iflag
     3409!AC!       enddo
     3410!AC!      enddo
     3411
     3412
     3413! ***    homogenize tendencies below cloud base    ***
    35113414
    35123415
    … …  
    35223425  END DO
    35233426
    3524   ! do i=1,nl
    3525   ! do il=1,ncum
    3526   ! th_wake(il,i)=t_wake(il,i)*(1000.0/p(il,i))**rdcp
    3527   ! enddo
    3528   ! enddo
     3427!do i=1,nl
     3428!do il=1,ncum
     3429!th_wake(il,i)=t_wake(il,i)*(1000.0/p(il,i))**rdcp
     3430!enddo
     3431!enddo
    35293432
    35303433  DO i = 1, nl
    35313434    DO il = 1, ncum
    35323435      IF (i<=(icb(il)-1) .AND. iflag(il)<=1) THEN
    3533         ! jyg  Saturated part : use T profile
     3436!jyg  Saturated part : use T profile
    35343437        asum(il) = asum(il) + (ft(il,i)-ftd(il,i))*(ph(il,i)-ph(il,i+1))
    3535         bsum(il) = bsum(il) + (fr(il,i)-fqd(il,i))*(lv(il,i)+(cl-cpd)*(t(il, &
    3536           i)-t(il,1)))*(ph(il,i)-ph(il,i+1))
    3537         csum(il) = csum(il) + (lv(il,i)+(cl-cpd)*(t(il,i)-t(il, &
    3538           1)))*(ph(il,i)-ph(il,i+1))
     3438!jyg<20140311
     3439!Correction pour conserver l eau
     3440        IF (ok_conserv_q) THEN
     3441          bsum(il) = bsum(il) + (fr(il,i)-fqd(il,i))*(ph(il,i)-ph(il,i+1))
     3442          csum(il) = csum(il) + (ph(il,i)-ph(il,i+1))
     3443
     3444        ELSE
     3445          bsum(il)=bsum(il)+(fr(il,i)-fqd(il,i))*(lv(il,i)+(cl-cpd)*(t(il,i)-t(il,1)))* &
     3446                            (ph(il,i)-ph(il,i+1))
     3447          csum(il)=csum(il)+(lv(il,i)+(cl-cpd)*(t(il,i)-t(il,1)))* &
     3448                            (ph(il,i)-ph(il,i+1))
     3449        ENDIF ! (ok_conserv_q)
     3450!jyg>
    35393451        dsum(il) = dsum(il) + t(il, i)*(ph(il,i)-ph(il,i+1))/th(il, i)
    3540         ! jyg  Unsaturated part : use T_wake profile
     3452!jyg  Unsaturated part : use T_wake profile
    35413453        esum(il) = esum(il) + ftd(il, i)*(ph(il,i)-ph(il,i+1))
    3542         fsum(il) = fsum(il) + fqd(il, i)*(lv(il,i)+(cl-cpd)*(t_wake(il, &
    3543           i)-t_wake(il,1)))*(ph(il,i)-ph(il,i+1))
    3544         gsum(il) = gsum(il) + (lv(il,i)+(cl-cpd)*(t_wake(il,i)-t_wake(il, &
    3545           1)))*(ph(il,i)-ph(il,i+1))
    3546         hsum(il) = hsum(il) + t_wake(il, i)*(ph(il,i)-ph(il,i+1))/th_wake(il, &
    3547           i)
     3454!jyg<20140311
     3455!Correction pour conserver l eau
     3456        IF (ok_conserv_q) THEN
     3457          fsum(il) = fsum(il) + fqd(il, i)*(ph(il,i)-ph(il,i+1))
     3458          gsum(il) = gsum(il) + (ph(il,i)-ph(il,i+1))
     3459        ELSE
     3460          fsum(il)=fsum(il)+fqd(il,i)*(lv(il,i)+(cl-cpd)*(t_wake(il,i)-t_wake(il,1)))* &
     3461                            (ph(il,i)-ph(il,i+1))
     3462          gsum(il)=gsum(il)+(lv(il,i)+(cl-cpd)*(t_wake(il,i)-t_wake(il,1)))* &
     3463                            (ph(il,i)-ph(il,i+1))
     3464        ENDIF ! (ok_conserv_q)
     3465!jyg>
     3466        hsum(il) = hsum(il) + t_wake(il, i)*(ph(il,i)-ph(il,i+1))/th_wake(il, i)
    35483467      END IF
    35493468    END DO
    35503469  END DO
    35513470
    3552   ! !!!      do 700 i=1,icb(il)-1
     3471!!!!      do 700 i=1,icb(il)-1
    35533472  DO i = 1, nl
    35543473    DO il = 1, ncum
    … …  
    35623481  END DO
    35633482
    3564 
    3565   ! ***   Check that moisture stays positive. If not, scale tendencies
    3566   ! in order to ensure moisture positivity
     3483!jyg<
     3484!Conservation de l'eau
     3485!!  sumdq = 0.
     3486!!  DO k = 1, nl
     3487!!    sumdq = sumdq + fr(1, k)*100.*(ph(1,k)-ph(1,k+1))/grav
     3488!!  END DO
     3489!!  PRINT *, 'cv3_yield, apres hom, sum(dq), precip, somme ', sumdq, Vprecip(1, 1), sumdq + vprecip(1, 1)
     3490!jyg>
     3491
     3492
     3493! ***   Check that moisture stays positive. If not, scale tendencies
     3494! in order to ensure moisture positivity
    35673495  DO il = 1, ncum
    35683496    alpha_qpos(il) = 1.
    35693497    IF (iflag(il)<=1) THEN
    35703498      IF (fr(il,1)<=0.) THEN
    3571         alpha_qpos(il) = max(alpha_qpos(il), (-delt*fr(il, &
    3572           1))/(s_wake(il)*rr_wake(il,1)+(1.-s_wake(il))*rr(il,1)))
     3499        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)))
    35733500      END IF
    35743501    END IF
    … …  
    35783505      IF (iflag(il)<=1) THEN
    35793506        IF (fr(il,i)<=0.) THEN
    3580           alpha_qpos1(il) = max(1., (-delt*fr(il,i))/(s_wake(il)*rr_wake(il, &
    3581             i)+(1.-s_wake(il))*rr(il,i)))
    3582           IF (alpha_qpos1(il)>=alpha_qpos(il)) alpha_qpos(il) &
    3583             = alpha_qpos1(il)
     3507          alpha_qpos1(il) = max(1., (-delt*fr(il,i))/(s_wake(il)*rr_wake(il,i)+(1.-s_wake(il))*rr(il,i)))
     3508          IF (alpha_qpos1(il)>=alpha_qpos(il)) alpha_qpos(il) = alpha_qpos1(il)
    35843509        END IF
    35853510      END IF
    … …  
    36083533        m(il, i) = m(il, i)/alpha_qpos(il)
    36093534        mp(il, i) = mp(il, i)/alpha_qpos(il)
    3610         vprecip(il, i) = vprecip(il, i)/alpha_qpos(il)
     3535        Vprecip(il, i) = vprecip(il, i)/alpha_qpos(il)
    36113536      END IF
    36123537    END DO
    … …  
    36223547  END DO
    36233548
    3624   ! AC!      DO j = 1,ntra
    3625   ! AC!      DO i = 1,nl
    3626   ! AC!       DO il = 1,ncum
    3627   ! AC!        IF (iflag(il) .le. 1) THEN
    3628   ! AC!         ftra(il,i,j) = ftra(il,i,j)/alpha_qpos(il)
    3629   ! AC!        ENDIF
    3630   ! AC!       ENDDO
    3631   ! AC!      ENDDO
    3632   ! AC!      ENDDO
    3633 
    3634 
    3635   ! ***           reset counter and return           ***
     3549!AC!      DO j = 1,ntra
     3550!AC!      DO i = 1,nl
     3551!AC!       DO il = 1,ncum
     3552!AC!        IF (iflag(il) .le. 1) THEN
     3553!AC!         ftra(il,i,j) = ftra(il,i,j)/alpha_qpos(il)
     3554!AC!        ENDIF
     3555!AC!       ENDDO
     3556!AC!      ENDDO
     3557!AC!      ENDDO
     3558
     3559
     3560! ***           reset counter and return           ***
    36363561
    36373562  DO il = 1, ncum
    … …  
    37023627          END IF
    37033628        END IF
    3704         ! c        print *,'cbmf',il,i,k,cbmf(il),wghti(il,k)
     3629! c        print *,'cbmf',il,i,k,cbmf(il),wghti(il,k)
    37053630      END DO
    37063631    END DO
    … …  
    37103635    DO k = i, nl
    37113636      DO il = 1, ncum
    3712         ! test         if (i.ge.icb(il).and.i.le.inb(il).and.k.le.inb(il))
    3713         ! then
     3637! test         if (i.ge.icb(il).and.i.le.inb(il).and.k.le.inb(il)) then
    37143638        IF (i<=inb(il) .AND. k<=inb(il)) THEN
    37153639          upwd(il, i) = upwd(il, i) + up1(il, k, i)
    37163640          dnwd(il, i) = dnwd(il, i) + dn1(il, k, i)
    37173641        END IF
    3718         ! c         print
    3719         ! *,'upwd',il,i,k,inb(il),upwd(il,i),m(il,k),up1(il,k,i)
     3642! c         print *,'upwd',il,i,k,inb(il),upwd(il,i),m(il,k),up1(il,k,i)
    37203643      END DO
    37213644    END DO
    … …  
    37233646
    37243647
    3725   ! !!!      DO il=1,ncum
    3726   ! !!!      do i=icb(il),inb(il)
    3727   ! !!!
    3728   ! !!!      upwd(il,i)=0.0
    3729   ! !!!      dnwd(il,i)=0.0
    3730   ! !!!      do k=i,inb(il)
    3731   ! !!!      up1=0.0
    3732   ! !!!      dn1=0.0
    3733   ! !!!      do n=1,i-1
    3734   ! !!!      up1=up1+ment(il,n,k)
    3735   ! !!!      dn1=dn1-ment(il,k,n)
    3736   ! !!!      enddo
    3737   ! !!!      upwd(il,i)=upwd(il,i)+m(il,k)+up1
    3738   ! !!!      dnwd(il,i)=dnwd(il,i)+dn1
    3739   ! !!!      enddo
    3740   ! !!!      enddo
    3741   ! !!!
    3742   ! !!!      ENDDO
    3743 
    3744   ! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    3745   ! determination de la variation de flux ascendant entre
    3746   ! deux niveau non dilue mip
    3747   ! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
     3648!!!!      DO il=1,ncum
     3649!!!!      do i=icb(il),inb(il)
     3650!!!!
     3651!!!!      upwd(il,i)=0.0
     3652!!!!      dnwd(il,i)=0.0
     3653!!!!      do k=i,inb(il)
     3654!!!!      up1=0.0
     3655!!!!      dn1=0.0
     3656!!!!      do n=1,i-1
     3657!!!!      up1=up1+ment(il,n,k)
     3658!!!!      dn1=dn1-ment(il,k,n)
     3659!!!!      enddo
     3660!!!!      upwd(il,i)=upwd(il,i)+m(il,k)+up1
     3661!!!!      dnwd(il,i)=dnwd(il,i)+dn1
     3662!!!!      enddo
     3663!!!!      enddo
     3664!!!!
     3665!!!!      ENDDO
     3666
     3667! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
     3668! determination de la variation de flux ascendant entre
     3669! deux niveau non dilue mip
     3670! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    37483671
    37493672  DO i = 1, nl
    … …  
    37873710  END DO
    37883711
    3789   ! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    3790   ! icb represente de niveau ou se trouve la
    3791   ! base du nuage , et inb le top du nuage
    3792   ! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
     3712! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
     3713! icb represente de niveau ou se trouve la
     3714! base du nuage , et inb le top du nuage
     3715! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    37933716
    37943717  DO i = 1, nd
    … …  
    38003723  DO i = 1, nd
    38013724    DO il = 1, ncum
    3802       rdcp = (rrd*(1.-rr(il,i))-rr(il,i)*rrv)/(cpd*(1.-rr(il, &
    3803         i))+rr(il,i)*cpv)
     3725      rdcp = (rrd*(1.-rr(il,i))-rr(il,i)*rrv)/(cpd*(1.-rr(il,i))+rr(il,i)*cpv)
    38043726      tls(il, i) = t(il, i)*(1000.0/p(il,i))**rdcp
    38053727      tps(il, i) = tp(il, i)
    … …  
    38083730
    38093731
    3810   ! *** diagnose the in-cloud mixing ratio   ***            ! cld
    3811   ! ***           of condensed water         ***            ! cld
    3812   ! ! cld
    3813 
    3814   DO i = 1, nd ! cld
    3815     DO il = 1, ncum ! cld
    3816       mac(il, i) = 0.0 ! cld
    3817       wa(il, i) = 0.0 ! cld
    3818       siga(il, i) = 0.0 ! cld
    3819       sax(il, i) = 0.0 ! cld
    3820     END DO ! cld
    3821   END DO ! cld
    3822 
    3823   DO i = minorig, nl ! cld
    3824     DO k = i + 1, nl + 1 ! cld
    3825       DO il = 1, ncum ! cld
     3732! *** diagnose the in-cloud mixing ratio   ***                       ! cld
     3733! ***           of condensed water         ***                       ! cld
     3734!! cld                                                               
     3735                                                                     
     3736  DO i = 1, nd                                                       ! cld
     3737    DO il = 1, ncum                                                  ! cld
     3738      mac(il, i) = 0.0                                               ! cld
     3739      wa(il, i) = 0.0                                                ! cld
     3740      siga(il, i) = 0.0                                              ! cld
     3741      sax(il, i) = 0.0                                               ! cld
     3742    END DO                                                           ! cld
     3743  END DO                                                             ! cld
     3744                                                                     
     3745  DO i = minorig, nl                                                 ! cld
     3746    DO k = i + 1, nl + 1                                             ! cld
     3747      DO il = 1, ncum                                                ! cld
    38263748        IF (i<=inb(il) .AND. k<=(inb(il)+1) .AND. iflag(il)<=1) THEN ! cld
    3827           mac(il, i) = mac(il, i) + m(il, k) ! cld
    3828         END IF ! cld
    3829       END DO ! cld
    3830     END DO ! cld
    3831   END DO ! cld
    3832 
    3833   DO i = 1, nl ! cld
    3834     DO j = 1, i ! cld
    3835       DO il = 1, ncum ! cld
    3836         IF (i>=icb(il) .AND. i<=(inb(il)-1) & ! cld
    3837             .AND. j>=icb(il) .AND. iflag(il)<=1) THEN ! cld
    3838           sax(il, i) = sax(il, i) + rrd*(tvp(il,j)-tv(il,j)) & ! cld
    3839             *(ph(il,j)-ph(il,j+1))/p(il, j) ! cld
    3840         END IF ! cld
    3841       END DO ! cld
    3842     END DO ! cld
    3843   END DO ! cld
    3844 
    3845   DO i = 1, nl ! cld
    3846     DO il = 1, ncum ! cld
    3847       IF (i>=icb(il) .AND. i<=(inb(il)-1) & ! cld
    3848           .AND. sax(il,i)>0.0 .AND. iflag(il)<=1) THEN ! cld
    3849         wa(il, i) = sqrt(2.*sax(il,i)) ! cld
    3850       END IF ! cld
    3851     END DO ! cld
    3852   END DO ! cld
    3853 
    3854   DO i = 1, nl ! cld
    3855     DO il = 1, ncum ! cld
    3856       IF (wa(il,i)>0.0 .AND. iflag(il)<=1) & ! cld
    3857         siga(il, i) = mac(il, i)/wa(il, i) & ! cld
    3858         *rrd*tvp(il, i)/p(il, i)/100./delta ! cld
    3859       siga(il, i) = min(siga(il,i), 1.0) ! cld
    3860       ! IM cf. FH
    3861       IF (iflag_clw==0) THEN
    3862         qcondc(il, i) = siga(il, i)*clw(il, i)*(1.-ep(il,i)) & ! cld
    3863           +(1.-siga(il,i))*qcond(il, i) ! cld
    3864       ELSE IF (iflag_clw==1) THEN
    3865         qcondc(il, i) = qcond(il, i) ! cld
    3866       END IF
    3867 
    3868     END DO ! cld
    3869   END DO
    3870   ! print*,'cv3_yield fin'
    3871     ! cld
     3749          mac(il, i) = mac(il, i) + m(il, k)                         ! cld
     3750        END IF                                                       ! cld
     3751      END DO                                                         ! cld
     3752    END DO                                                           ! cld
     3753  END DO                                                             ! cld
     3754
     3755  DO i = 1, nl                                                       ! cld
     3756    DO j = 1, i                                                      ! cld
     3757      DO il = 1, ncum                                                ! cld
     3758        IF (i>=icb(il) .AND. i<=(inb(il)-1) &                        ! cld
     3759            .AND. j>=icb(il) .AND. iflag(il)<=1) THEN                ! cld
     3760          sax(il, i) = sax(il, i) + rrd*(tvp(il,j)-tv(il,j)) &       ! cld
     3761            *(ph(il,j)-ph(il,j+1))/p(il, j)                          ! cld
     3762        END IF                                                       ! cld
     3763      END DO                                                         ! cld
     3764    END DO                                                           ! cld
     3765  END DO                                                             ! cld
     3766
     3767  DO i = 1, nl                                                       ! cld
     3768    DO il = 1, ncum                                                  ! cld
     3769      IF (i>=icb(il) .AND. i<=(inb(il)-1) &                          ! cld
     3770          .AND. sax(il,i)>0.0 .AND. iflag(il)<=1) THEN               ! cld
     3771        wa(il, i) = sqrt(2.*sax(il,i))                               ! cld
     3772      END IF                                                         ! cld
     3773    END DO                                                           ! cld
     3774  END DO                                                             ! cld
     3775
     3776  DO i = 1, nl                                                       ! cld
     3777    DO il = 1, ncum                                                  ! cld
     3778      IF (wa(il,i)>0.0 .AND. iflag(il)<=1) &                         ! cld
     3779        siga(il, i) = mac(il, i)/wa(il, i) &                         ! cld
     3780        *rrd*tvp(il, i)/p(il, i)/100./delta                          ! cld
     3781      siga(il, i) = min(siga(il,i), 1.0)                             ! cld
     3782! IM cf. FH                                                         
     3783      IF (iflag_clw==0) THEN                                         ! cld
     3784        qcondc(il, i) = siga(il, i)*clw(il, i)*(1.-ep(il,i)) &       ! cld
     3785          +(1.-siga(il,i))*qcond(il, i)                              ! cld
     3786      ELSE IF (iflag_clw==1) THEN                                    ! cld
     3787        qcondc(il, i) = qcond(il, i)                                 ! cld
     3788      END IF                                                         ! cld
     3789
     3790    END DO                                                           ! cld
     3791  END DO
     3792! print*,'cv3_yield fin'
     3793
    38723794  RETURN
    38733795END SUBROUTINE cv3_yield
    38743796
    3875 ! AC! et !RomP >>>
    3876 SUBROUTINE cv3_tracer(nloc, len, ncum, nd, na, ment, sigij, da, phi, phi2, &
    3877     d1a, dam, ep, vprecip, elij, clw, epmlmmm, eplamm, icb, inb)
     3797!AC! et !RomP >>>
     3798SUBROUTINE cv3_tracer(nloc, len, ncum, nd, na, &
     3799                      ment, sigij, da, phi, phi2, d1a, dam, &
     3800                      ep, Vprecip, elij, clw, epmlmMm, eplaMm, &
     3801                      icb, inb)
    38783802  IMPLICIT NONE
    38793803
    38803804  include "cv3param.h"
    38813805
    3882   ! inputs:
     3806!inputs:
    38833807  INTEGER ncum, nd, na, nloc, len
    38843808  REAL ment(nloc, na, na), sigij(nloc, na, na)
    … …  
    38863810  REAL ep(nloc, na)
    38873811  INTEGER icb(nloc), inb(nloc)
    3888   REAL vprecip(nloc, nd+1)
    3889   ! ouputs:
     3812  REAL Vprecip(nloc, nd+1)
     3813!ouputs:
    38903814  REAL da(nloc, na), phi(nloc, na, na)
    38913815  REAL phi2(nloc, na, na)
    38923816  REAL d1a(nloc, na), dam(nloc, na)
    3893   REAL epmlmmm(nloc, na, na), eplamm(nloc, na)
    3894   ! variables pour tracer dans precip de l'AA et des mel
    3895   ! local variables:
     3817  REAL epmlmMm(nloc, na, na), eplaMm(nloc, na)
     3818! variables pour tracer dans precip de l'AA et des mel
     3819!local variables:
    38963820  INTEGER i, j, k
    38973821  REAL epm(nloc, na, na)
    38983822
    3899   ! variables d'Emanuel : du second indice au troisieme
    3900   ! --->    tab(i,k,j) -> de l origine k a l arrivee j
    3901   ! ment, sigij, elij
    3902   ! variables personnelles : du troisieme au second indice
    3903   ! --->    tab(i,j,k) -> de k a j
    3904   ! phi, phi2
    3905 
    3906   ! initialisations
     3823! variables d'Emanuel : du second indice au troisieme
     3824! --->    tab(i,k,j) -> de l origine k a l arrivee j
     3825! ment, sigij, elij
     3826! variables personnelles : du troisieme au second indice
     3827! --->    tab(i,j,k) -> de k a j
     3828! phi, phi2
     3829
     3830! initialisations
    39073831
    39083832  da(:, :) = 0.
    … …  
    39103834  dam(:, :) = 0.
    39113835  epm(:, :, :) = 0.
    3912   eplamm(:, :) = 0.
    3913   epmlmmm(:, :, :) = 0.
     3836  eplaMm(:, :) = 0.
     3837  epmlmMm(:, :, :) = 0.
    39143838  phi(:, :, :) = 0.
    39153839  phi2(:, :, :) = 0.
    39163840
    3917   ! fraction deau condensee dans les melanges convertie en precip : epm
    3918   ! et eau condensée précipitée dans masse d'air saturé : l_m*dM_m/dzdz.dzdz
     3841! fraction deau condensee dans les melanges convertie en precip : epm
     3842! et eau condensée précipitée dans masse d'air saturé : l_m*dM_m/dzdz.dzdz
    39193843  DO j = 1, na
    39203844    DO k = 1, na
    39213845      DO i = 1, ncum
    3922         IF (k>=icb(i) .AND. k<=inb(i) .AND. & ! !jyg     &
    3923                                               ! j.ge.k.and.j.le.inb(i)) then
    3924           ! !jyg             epm(i,j,k)=1.-(1.-ep(i,j))*clw(i,j)/elij(i,k,j)
     3846        IF (k>=icb(i) .AND. k<=inb(i) .AND. &
     3847!!jyg              j.ge.k.and.j.le.inb(i)) then
     3848!!jyg             epm(i,j,k)=1.-(1.-ep(i,j))*clw(i,j)/elij(i,k,j)
    39253849            j>k .AND. j<=inb(i)) THEN
    39263850          epm(i, j, k) = 1. - (1.-ep(i,j))*clw(i, j)/max(elij(i,k,j), 1.E-16)
    3927           ! !
     3851!!
    39283852          epm(i, j, k) = max(epm(i,j,k), 0.0)
    39293853        END IF
    … …  
    39373861      DO i = 1, ncum
    39383862        IF (k>=icb(i) .AND. k<=inb(i)) THEN
    3939           eplamm(i, j) = eplamm(i, j) + ep(i, j)*clw(i, j)*ment(i, j, k)*(1.- &
    3940             sigij(i,j,k))
     3863          eplaMm(i, j) = eplamm(i, j) + &
     3864                         ep(i, j)*clw(i, j)*ment(i, j, k)*(1.-sigij(i,j,k))
    39413865        END IF
    39423866      END DO
    … …  
    39483872      DO i = 1, ncum
    39493873        IF (k>=icb(i) .AND. k<=inb(i) .AND. j<=inb(i)) THEN
    3950           epmlmmm(i, j, k) = epm(i, j, k)*elij(i, k, j)*ment(i, k, j)
     3874          epmlmMm(i, j, k) = epm(i, j, k)*elij(i, k, j)*ment(i, k, j)
    39513875        END IF
    39523876      END DO
    … …  
    39543878  END DO
    39553879
    3956   ! matrices pour calculer la tendance des concentrations dans cvltr.F90
     3880! matrices pour calculer la tendance des concentrations dans cvltr.F90
    39573881  DO j = 1, na
    39583882    DO k = 1, na
    … …  
    39623886        d1a(i, j) = d1a(i, j) + ment(i, k, j)*ep(i, k)*(1.-sigij(i,k,j))
    39633887        IF (k<=j) THEN
    3964           dam(i, j) = dam(i, j) + ment(i, k, j)*epm(i, k, j)*(1.-ep(i,k))*(1. &
    3965             -sigij(i,k,j))
    3966 
     3888          dam(i, j) = dam(i, j) + ment(i, k, j)*epm(i, k, j)*(1.-ep(i,k))*(1.-sigij(i,k,j))
    39673889          phi2(i, j, k) = phi(i, j, k)*epm(i, j, k)
    39683890        END IF
    … …  
    39733895  RETURN
    39743896END SUBROUTINE cv3_tracer
    3975 ! AC! et !RomP <<<
    3976 
    3977 SUBROUTINE cv3_uncompress(nloc, len, ncum, nd, ntra, idcum, iflag, precip, &
    3978     sig, w0, ft, fq, fu, fv, ftra, ma, upwd, dnwd, dnwd0, qcondc, wd, cape, &
    3979     iflag1, precip1, sig1, w01, ft1, fq1, fu1, fv1, ftra1, ma1, upwd1, dnwd1, &
    3980     dnwd01, qcondc1, wd1, cape1)
     3897!AC! et !RomP <<<
     3898
     3899SUBROUTINE cv3_uncompress(nloc, len, ncum, nd, ntra, idcum, &
     3900                          iflag, &
     3901                          precip, sig, w0, &
     3902                          ft, fq, fu, fv, ftra, &
     3903                          Ma, upwd, dnwd, dnwd0, qcondc, wd, cape, &
     3904                          iflag1, &
     3905                          precip1, sig1, w01, &
     3906                          ft1, fq1, fu1, fv1, ftra1, &
     3907                          Ma1, upwd1, dnwd1, dnwd01, qcondc1, wd1, cape1)
    39813908  IMPLICIT NONE
    39823909
    39833910  include "cv3param.h"
    39843911
    3985   ! inputs:
     3912!inputs:
    39863913  INTEGER len, ncum, nd, ntra, nloc
    39873914  INTEGER idcum(nloc)
    … …  
    39963923  REAL wd(nloc), cape(nloc)
    39973924
    3998   ! outputs:
     3925!outputs:
    39993926  INTEGER iflag1(len)
    40003927  REAL precip1(len)
    … …  
    40073934  REAL wd1(nloc), cape1(nloc)
    40083935
    4009   ! local variables:
     3936!local variables:
    40103937  INTEGER i, k, j
    40113938
    … …  
    40383965
    40393966
    4040   ! AC!        do 2100 j=1,ntra
    4041   ! AC!c oct3         do 2110 k=1,nl
    4042   ! AC!         do 2110 k=1,nd ! oct3
    4043   ! AC!          do 2120 i=1,ncum
    4044   ! AC!            ftra1(idcum(i),k,j)=ftra(i,k,j)
    4045   ! AC! 2120     continue
    4046   ! AC! 2110    continue
    4047   ! AC! 2100   continue
     3967!AC!        do 2100 j=1,ntra
     3968!AC!c oct3         do 2110 k=1,nl
     3969!AC!         do 2110 k=1,nd ! oct3
     3970!AC!          do 2120 i=1,ncum
     3971!AC!            ftra1(idcum(i),k,j)=ftra(i,k,j)
     3972!AC! 2120     continue
     3973!AC! 2110    continue
     3974!AC! 2100   continue
     3975!
    40483976  RETURN
    40493977END SUBROUTINE cv3_uncompress
Note: See TracChangeset for help on using the changeset viewer.