Ignore:
Timestamp:
Jun 11, 2014, 3:46:46 PM (10 years ago)
Author:
Laurent Fairhead
Message:

Merged trunk changes r1997:2055 into testing branch

Location:
LMDZ5/branches/testing
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • LMDZ5/branches/testing

  • LMDZ5/branches/testing/libf/phylmd/cv3_routines.F90

    r1999 r2056  
    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! -- "microphysical" parameters:
    5959    sigdz = 0.01
    6060    spfac = 0.15
    6161    pbcrit = 150.0
    6262    ptcrit = 500.0
    63     ! IM beg: ajout fis. reglage ep
     63! IM beg: ajout fis. reglage ep
    6464    flag_epkeorig = 1
    6565    elcrit = 0.0003
    6666    tlcrit = -55.0
    67     ! IM lu dans physiq.def via conf_phys.F90     epmax  = 0.993
     67! IM lu dans physiq.def via conf_phys.F90     epmax  = 0.993
    6868
    6969    omtrain = 45.0 ! used also for snow (no disctinction rain/snow)
    7070
    71     ! -- misc:
     71! -- misc:
    7272
    7373    dtovsh = -0.2 ! dT for overshoot
    7474    dpbase = -40. ! definition cloud base (400m above LCL)
    75     ! cc      dttrig = 5.   ! (loose) condition for triggering
     75! cc      dttrig = 5.   ! (loose) condition for triggering
    7676    dttrig = 10. ! (loose) condition for triggering
    7777    flag_wb = 1
    7878    wbmax = 6. ! (m/s) adiab updraught speed at LFC (used in cv3p1_closure)
    7979
    80     ! -- rate of approach to quasi-equilibrium:
     80! -- rate of approach to quasi-equilibrium:
    8181
    8282    dtcrit = -2.0
    8383    tau = 8000.
    8484
    85     ! -- interface cloud parameterization:
     85! -- interface cloud parameterization:
    8686
    8787    delta = 0.01 ! cld
    8888
    89     ! -- interface with boundary-layer (gust factor): (sb)
     89! -- interface with boundary-layer (gust factor): (sb)
    9090
    9191    betad = 10.0 ! original value (from convect 4.3)
    9292
    93     OPEN (99, FILE='conv_param.data', STATUS='old', FORM='formatted', &
    94       ERR=9999)
     93    OPEN (99, FILE='conv_param.data', STATUS='old', FORM='formatted', ERR=9999)
    9594    READ (99, *, END=9998) dpbase
    9695    READ (99, *, END=9998) pbcrit
    … …  
    113112    WRITE (*, *) 'wbmax =', wbmax
    114113
    115     ! IM Lecture du fichier ep_param.data
     114! IM Lecture du fichier ep_param.data
    116115    OPEN (79, FILE='ep_param.data', STATUS='old', FORM='formatted', ERR=7999)
    117116    READ (79, *, END=7998) flag_epkeorig
    … …  
    124123    WRITE (*, *) 'elcrit=', elcrit
    125124    WRITE (*, *) 'tlcrit=', tlcrit
    126     ! IM end: ajout fis. reglage ep
     125! IM end: ajout fis. reglage ep
    127126
    128127    first = .FALSE.
    129   END IF
    130 
     128
     129  END IF ! (first)
     130
     131! print*,'tau=',tau
    131132  beta = 1.0 - delt/tau
    132133  alpha1 = 1.5E-3
    133   ! jyg    Correction bug alpha
     134!JYG    Correction bug alpha
    134135  alpha1 = alpha1*1.5
    135136  alpha = alpha1*delt/tau
    136   ! jyg    Bug
    137   ! cc increase alpha to compensate W decrease:
    138   ! c      alpha  = alpha*1.5
     137!JYG    Bug
     138! cc increase alpha to compensate W decrease:
     139! c      alpha  = alpha*1.5
    139140
    140141  RETURN
    141142END SUBROUTINE cv3_param
    142143
    143 SUBROUTINE cv3_prelim(len, nd, ndp1, t, q, p, ph, lv, lf, cpn, tv, gz, h, hm, &
    144     th)
     144SUBROUTINE cv3_prelim(len, nd, ndp1, t, q, p, ph, &
     145                      lv, lf, cpn, tv, gz, h, hm, th)
    145146  IMPLICIT NONE
    146147
    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:
     148! =====================================================================
     149! --- CALCULATE ARRAYS OF GEOPOTENTIAL, HEAT CAPACITY & STATIC ENERGY
     150! "ori": from convect4.3 (vectorized)
     151! "convect3": to be exactly consistent with convect3
     152! =====================================================================
     153
     154! inputs:
    154155  INTEGER len, nd, ndp1
    155156  REAL t(len, nd), q(len, nd), p(len, nd), ph(len, ndp1)
    156157
    157   ! outputs:
     158! outputs:
    158159  REAL lv(len, nd), lf(len, nd), cpn(len, nd), tv(len, nd)
    159160  REAL gz(len, nd), h(len, nd), hm(len, nd)
    160161  REAL th(len, nd)
    161162
    162   ! local variables:
     163! local variables:
    163164  INTEGER k, i
    164165  REAL rdcp
    … …  
    170171
    171172
    172   ! ori      do 110 k=1,nlp
    173   ! abderr     do 110 k=1,nl ! convect3
     173! ori      do 110 k=1,nlp
     174! abderr     do 110 k=1,nl ! convect3
    174175  DO k = 1, nlp
    175176
    176177    DO i = 1, len
    177       ! debug          lv(i,k)= lv0-clmcpv*(t(i,k)-t0)
     178! debug          lv(i,k)= lv0-clmcpv*(t(i,k)-t0)
    178179      lv(i, k) = lv0 - clmcpv*(t(i,k)-273.15)
    179180      lf(i, k) = lf0 - clmci*(t(i,k)-273.15)
    180181      cpn(i, k) = cpd*(1.0-q(i,k)) + cpv*q(i, k)
    181182      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)
     183! ori          tv(i,k)=t(i,k)*(1.0+q(i,k)*epsim1)
    183184      tv(i, k) = t(i, k)*(1.0+q(i,k)/eps-q(i,k))
    184185      rdcp = (rrd*(1.-q(i,k))+q(i,k)*rrv)/cpn(i, k)
    … …  
    187188  END DO
    188189
    189   ! gz = phi at the full levels (same as p).
     190! gz = phi at the full levels (same as p).
    190191
    191192  DO i = 1, len
    192193    gz(i, 1) = 0.0
    193194  END DO
    194   ! ori      do 140 k=2,nlp
     195! ori      do 140 k=2,nlp
    195196  DO k = 2, nl ! convect3
    196197    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
     198      tvx = t(i, k)*(1.+q(i,k)/eps-q(i,k))         !convect3
     199      tvy = t(i, k-1)*(1.+q(i,k-1)/eps-q(i,k-1))   !convect3
     200      gz(i, k) = gz(i, k-1) + 0.5*rrd*(tvx+tvy)* & !convect3
     201                 (p(i,k-1)-p(i,k))/ph(i, k)        !convect3
     202
     203! c        print *,' gz(',k,')',gz(i,k),' tvx',tvx,' tvy ',tvy
     204
     205! ori         gz(i,k)=gz(i,k-1)+hrd*(tv(i,k-1)+tv(i,k))
     206! ori    &         *(p(i,k-1)-p(i,k))/ph(i,k)
     207    END DO
     208  END DO
     209
     210! h  = phi + cpT (dry static energy).
     211! hm = phi + cp(T-Tbase)+Lq
     212
     213! ori      do 170 k=1,nlp
    213214  DO k = 1, nl ! convect3
    214215    DO i = 1, len
    … …  
    221222END SUBROUTINE cv3_prelim
    222223
    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)
     224SUBROUTINE cv3_feed(len, nd, ok_conserv_q, &
     225                    t, q, u, v, p, ph, hm, gz, &
     226                    p1feed, p2feed, wght, &
     227                    wghti, tnk, thnk, qnk, qsnk, unk, vnk, &
     228                    cpnk, hnk, nk, icb, icbmax, iflag, gznk, plcl)
    226229  IMPLICIT NONE
    227230
    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   ! ================================================================
     231! ================================================================
     232! Purpose: CONVECTIVE FEED
     233
     234! Main differences with cv_feed:
     235! - ph added in input
     236! - here, nk(i)=minorig
     237! - icb defined differently (plcl compared with ph instead of p)
     238
     239! Main differences with convect3:
     240! - we do not compute dplcldt and dplcldr of CLIFT anymore
     241! - values iflag different (but tests identical)
     242! - A,B explicitely defined (!...)
     243! ================================================================
    241244
    242245  include "cv3param.h"
    243246  include "cvthermo.h"
    244247
    245   ! inputs:
     248!inputs:
    246249  INTEGER len, nd
     250  LOGICAL ok_conserv_q
    247251  REAL t(len, nd), q(len, nd), p(len, nd)
    248252  REAL u(len, nd), v(len, nd)
    … …  
    250254  REAL ph(len, nd+1)
    251255  REAL p1feed(len)
    252   ! ,  wght(len)
     256! ,  wght(len)
    253257  REAL wght(nd)
    254   ! input-output
     258!input-output
    255259  REAL p2feed(len)
    256   ! outputs:
     260!outputs:
    257261  INTEGER iflag(len), nk(len), icb(len), icbmax
    258   ! real   wghti(len)
     262!  real   wghti(len)
    259263  REAL wghti(len, nd)
    260264  REAL tnk(len), thnk(len), qnk(len), qsnk(len)
    … …  
    263267  REAL plcl(len)
    264268
    265   ! local variables:
     269!local variables:
    266270  INTEGER i, k, iter, niter
    267271  INTEGER ihmin(len)
    … …  
    269273  REAL pup(len), plo(len), pfeed(len)
    270274  REAL plclup(len), plcllo(len), plclfeed(len)
     275  REAL pfeedmin(len)
    271276  REAL posit(len)
    272277  LOGICAL nocond(len)
    273278
    274   ! -------------------------------------------------------------------
    275   ! --- Origin level of ascending parcels for convect3:
    276   ! -------------------------------------------------------------------
     279!jyg20140217<
     280  INTEGER iostat
     281  LOGICAL, SAVE :: first
     282  LOGICAL, SAVE :: ok_new_feed
     283  REAL, SAVE :: dp_lcl_feed
     284!$OMP THREADPRIVATE (first,ok_new_feed,dp_lcl_feed)
     285  DATA first/.TRUE./
     286  DATA dp_lcl_feed/2./
     287
     288  IF (first) THEN
     289!$OMP MASTER
     290    ok_new_feed = ok_conserv_q
     291    OPEN (98, FILE='cv3feed_param.data', STATUS='old', FORM='formatted', IOSTAT=iostat)
     292    IF (iostat==0) THEN
     293      READ (98, *, END=998) ok_new_feed
     294998   CONTINUE
     295      CLOSE (98)
     296    END IF
     297    PRINT *, ' ok_new_feed: ', ok_new_feed
     298    first = .FALSE.
     299!$OMP END MASTER
     300  END IF
     301!jyg>
     302! -------------------------------------------------------------------
     303! --- Origin level of ascending parcels for convect3:
     304! -------------------------------------------------------------------
    277305
    278306  DO i = 1, len
    … …  
    281309  END DO
    282310
    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
     311! -------------------------------------------------------------------
     312! --- Adjust feeding layer thickness so that lifting up to the top of
     313! --- the feeding layer does not induce condensation (i.e. so that
     314! --- plcl < p2feed).
     315! --- Method : iterative secant method.
     316! -------------------------------------------------------------------
     317
     318! 1- First bracketing of the solution : ph(nk+1), p2feed
     319
     320! 1.a- LCL associated with p2feed
    293321  DO i = 1, len
    294322    pup(i) = p2feed(i)
    295323  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)
     324  CALL cv3_vertmix(len, nd, iflag, p1feed, pup, p, ph, &
     325                   t, q, u, v, wght, &
     326                   wghti, nk, tnk, thnk, qnk, qsnk, unk, vnk, plclup)
     327! 1.b- LCL associated with ph(nk+1)
    299328  DO i = 1, len
    300329    plo(i) = ph(i, nk(i)+1)
    301330  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
     331  CALL cv3_vertmix(len, nd, iflag, p1feed, plo, p, ph, &
     332                   t, q, u, v, wght, &
     333                   wghti, nk, tnk, thnk, qnk, qsnk, unk, vnk, plcllo)
     334! 2- Iterations
    305335  niter = 5
    306336  DO iter = 1, niter
    … …  
    314344        pfeed(i) = pup(i)
    315345      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))
     346!JYG20140217<
     347        IF (ok_new_feed) THEN
     348          pfeed(i) = (pup(i)*(plo(i)-plcllo(i)-dp_lcl_feed)+  &
     349                      plo(i)*(plclup(i)-pup(i)+dp_lcl_feed))/ &
     350                     (plo(i)-plcllo(i)+plclup(i)-pup(i))
     351        ELSE
     352          pfeed(i) = (pup(i)*(plo(i)-plcllo(i))+  &
     353                      plo(i)*(plclup(i)-pup(i)))/ &
     354                     (plo(i)-plcllo(i)+plclup(i)-pup(i))
     355        END IF
     356!JYG>
    318357      END IF
    319358    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)
     359!jyg20140217<
     360! For the last iteration, make sure that the top of the feeding layer
     361! and LCL are not in the same layer:
     362    IF (ok_new_feed) THEN
     363      IF (iter==niter) THEN
     364        DO k = minorig, nd
     365          DO i = 1, len
     366            IF (ph(i,k)>=plclfeed(i)) pfeedmin(i) = ph(i, k)
     367          END DO
     368        END DO
     369        DO i = 1, len
     370          pfeed(i) = max(pfeedmin(i), pfeed(i))
     371        END DO
     372      END IF
     373    END IF
     374!jyg>
     375
     376    CALL cv3_vertmix(len, nd, iflag, p1feed, pfeed, p, ph, &
     377                   t, q, u, v, wght, &
     378                   wghti, nk, tnk, thnk, qnk, qsnk, unk, vnk, plclfeed)
     379!jyg20140217<
     380    IF (ok_new_feed) THEN
     381      DO i = 1, len
     382        posit(i) = (sign(1.,plclfeed(i)-pfeed(i)+dp_lcl_feed)+1.)*0.5
     383        IF (plclfeed(i)-pfeed(i)+dp_lcl_feed==0.) posit(i) = 1.
     384      END DO
     385    ELSE
     386      DO i = 1, len
     387        posit(i) = (sign(1.,plclfeed(i)-pfeed(i))+1.)*0.5
     388        IF (plclfeed(i)==pfeed(i)) posit(i) = 1.
     389      END DO
     390    END IF
     391!jyg>
    322392    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
     393! - posit = 1 when lcl is below top of feeding layer (plclfeed>pfeed)
     394! -               => pup=pfeed
     395! - posit = 0 when lcl is above top of feeding layer (plclfeed<pfeed)
     396! -               => plo=pfeed
    329397      pup(i) = posit(i)*pfeed(i) + (1.-posit(i))*pup(i)
    330398      plo(i) = (1.-posit(i))*pfeed(i) + posit(i)*plo(i)
    … …  
    343411  END DO
    344412
    345   ! -------------------------------------------------------------------
    346   ! --- Check whether parcel level temperature and specific humidity
    347   ! --- are reasonable
    348   ! -------------------------------------------------------------------
     413! -------------------------------------------------------------------
     414! --- Check whether parcel level temperature and specific humidity
     415! --- are reasonable
     416! -------------------------------------------------------------------
    349417  DO i = 1, len
    350418    IF (((tnk(i)<250.0) .OR. (qnk(i)<=0.0)) .AND. (iflag(i)==0)) iflag(i) = 7
    351419  END DO
    352420
    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
     421! -------------------------------------------------------------------
     422! --- Calculate first level above lcl (=icb)
     423! -------------------------------------------------------------------
     424
     425!@      do 270 i=1,len
     426!@       icb(i)=nlm
     427!@ 270  continue
     428!@c
     429!@      do 290 k=minorig,nl
     430!@        do 280 i=1,len
     431!@          if((k.ge.(nk(i)+1)).and.(p(i,k).lt.plcl(i)))
     432!@     &    icb(i)=min(icb(i),k)
     433!@ 280    continue
     434!@ 290  continue
     435!@c
     436!@      do 300 i=1,len
     437!@        if((icb(i).ge.nlm).and.(iflag(i).eq.0))iflag(i)=9
     438!@ 300  continue
    371439
    372440  DO i = 1, len
    … …  
    374442  END DO
    375443
    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
     444! la modification consiste a comparer plcl a ph et non a p:
     445! icb est defini par :  ph(icb)<plcl<ph(icb-1)
     446!@      do 290 k=minorig,nl
    379447  DO k = 3, nl - 1 ! modif pour que icb soit sup/egal a 2
    380448    DO i = 1, len
    … …  
    384452
    385453
    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))
     454! print*,'icb dans cv3_feed '
     455! write(*,'(64i2)') icb(2:len-1)
     456! call dump2d(64,43,'plcl dans cv3_feed ',plcl(2:len-1))
    389457
    390458  DO i = 1, len
    391     ! @        if((icb(i).ge.nlm).and.(iflag(i).eq.0))iflag(i)=9
     459!@        if((icb(i).ge.nlm).and.(iflag(i).eq.0))iflag(i)=9
    392460    IF ((icb(i)==nlm) .AND. (iflag(i)==0)) iflag(i) = 9
    393461  END DO
    … …  
    397465  END DO
    398466
    399   ! Compute icbmax.
     467! Compute icbmax.
    400468
    401469  icbmax = 2
    402470  DO i = 1, len
    403     ! !        icbmax=max(icbmax,icb(i))
    404     IF (iflag(i)<7) icbmax = max(icbmax, icb(i)) ! sb Jun7th02
     471!!        icbmax=max(icbmax,icb(i))
     472    IF (iflag(i)<7) icbmax = max(icbmax, icb(i))     ! sb Jun7th02
    405473  END DO
    406474
    … …  
    409477
    410478SUBROUTINE cv3_undilute1(len, nd, t, qs, gz, plcl, p, icb, tnk, qnk, gznk, &
    411     tp, tvp, clw, icbs)
     479                         tp, tvp, clw, icbs)
    412480  IMPLICIT NONE
    413481
    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   ! ----------------------------------------------------------------
     482! ----------------------------------------------------------------
     483! Equivalent de TLIFT entre NK et ICB+1 inclus
     484
     485! Differences with convect4:
     486!    - specify plcl in input
     487!    - icbs is the first level above LCL (may differ from icb)
     488!    - in the iterations, used x(icbs) instead x(icb)
     489!    - many minor differences in the iterations
     490!    - tvp is computed in only one time
     491!    - icbs: first level above Plcl (IMIN de TLIFT) in output
     492!    - if icbs=icb, compute also tp(icb+1),tvp(icb+1) & clw(icb+1)
     493! ----------------------------------------------------------------
    426494
    427495  include "cvthermo.h"
    428496  include "cv3param.h"
    429497
    430   ! inputs:
     498! inputs:
    431499  INTEGER len, nd
    432500  INTEGER icb(len)
    … …  
    436504  REAL plcl(len) ! convect3
    437505
    438   ! outputs:
     506! outputs:
    439507  REAL tp(len, nd), tvp(len, nd), clw(len, nd)
    440508
    441   ! local variables:
     509! local variables:
    442510  INTEGER i, k
    443511  INTEGER icb1(len), icbs(len), icbsmax2 ! convect3
    … …  
    448516  REAL cpinv(len) ! convect3
    449517
    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   ***
     518! -------------------------------------------------------------------
     519! --- Calculates the lifted parcel virtual temperature at nk,
     520! --- the actual temperature, and the adiabatic
     521! --- liquid water content. The procedure is to solve the equation.
     522!    cp*tp+L*qp+phi=cp*tnk+L*qnk+gznk.
     523! -------------------------------------------------------------------
     524
     525
     526! ***  Calculate certain parcel quantities, including static energy   ***
    459527
    460528  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)
     529    ah0(i) = (cpd*(1.-qnk(i))+cl*qnk(i))*tnk(i) + qnk(i)*(lv0-clmcpv*(tnk(i)-273.15)) + gznk(i)
    463530    cpp(i) = cpd*(1.-qnk(i)) + qnk(i)*cpv
    464531    cpinv(i) = 1./cpp(i)
    465532  END DO
    466533
    467   ! ***   Calculate lifted parcel quantities below cloud base   ***
     534! ***   Calculate lifted parcel quantities below cloud base   ***
     535
     536  DO i = 1, len                                           !convect3
     537    icb1(i) = max(icb(i), 2)                              !convect3
     538    icb1(i) = min(icb(i), nl)                             !convect3
     539! if icb is below LCL, start loop at ICB+1:
     540! (icbs est le premier niveau au-dessus du LCL)
     541    icbs(i) = icb1(i)                                     !convect3
     542    IF (plcl(i)<p(i,icb1(i))) THEN
     543      icbs(i) = min(icbs(i)+1, nl)                        !convect3
     544    END IF
     545  END DO                                                  !convect3
    468546
    469547  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
     548    ticb(i) = t(i, icbs(i))                               !convect3
     549    gzicb(i) = gz(i, icbs(i))                             !convect3
     550    qsicb(i) = qs(i, icbs(i))                             !convect3
    478551  END DO !convect3
    479552
    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:
     553
     554! Re-compute icbsmax (icbsmax2):                          !convect3
     555!                                                         !convect3
     556  icbsmax2 = 2                                            !convect3
     557  DO i = 1, len                                           !convect3
     558    icbsmax2 = max(icbsmax2, icbs(i))                     !convect3
     559  END DO                                                  !convect3
     560
     561! initialization outputs:
     562
     563  DO k = 1, icbsmax2                                      ! convect3
     564    DO i = 1, len                                         ! convect3
     565      tp(i, k) = 0.0                                      ! convect3
     566      tvp(i, k) = 0.0                                     ! convect3
     567      clw(i, k) = 0.0                                     ! convect3
     568    END DO                                                ! convect3
     569  END DO                                                  ! convect3
     570
     571! tp and tvp below cloud base:
    505572
    506573  DO k = minorig, icbsmax2 - 1
    507574    DO i = 1, len
    508575      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    ***
     576      tvp(i, k) = tp(i, k)*(1.+qnk(i)/eps-qnk(i))        !whole thing (convect3)
     577    END DO
     578  END DO
     579
     580! ***  Find lifted parcel quantities above cloud base    ***
    514581
    515582  DO i = 1, len
    516583    tg = ticb(i)
    517     ! ori         qg=qs(i,icb(i))
     584! ori         qg=qs(i,icb(i))
    518585    qg = qsicb(i) ! convect3
    519     ! debug         alv=lv0-clmcpv*(ticb(i)-t0)
     586! debug         alv=lv0-clmcpv*(ticb(i)-t0)
    520587    alv = lv0 - clmcpv*(ticb(i)-273.15)
    521588
    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
     589! First iteration.
     590
     591! ori          s=cpd+alv*alv*qg/(rrv*ticb(i)*ticb(i))
     592    s = cpd*(1.-qnk(i)) + cl*qnk(i) + &                  ! convect3
     593        alv*alv*qg/(rrv*ticb(i)*ticb(i))                  ! convect3
    527594    s = 1./s
    528     ! ori          ahg=cpd*tg+(cl-cpd)*qnk(i)*ticb(i)+alv*qg+gzicb(i)
     595! ori          ahg=cpd*tg+(cl-cpd)*qnk(i)*ticb(i)+alv*qg+gzicb(i)
    529596    ahg = cpd*tg + (cl-cpd)*qnk(i)*tg + alv*qg + gzicb(i) ! convect3
    530597    tg = tg + s*(ah0(i)-ahg)
    531     ! ori          tg=max(tg,35.0)
    532     ! debug          tc=tg-t0
     598! ori          tg=max(tg,35.0)
     599! debug          tc=tg-t0
    533600    tc = tg - 273.15
    534601    denom = 243.5 + tc
    535602    denom = max(denom, 1.0) ! convect3
    536     ! ori          if(tc.ge.0.0)then
     603! ori          if(tc.ge.0.0)then
    537604    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))
     605! ori          else
     606! ori           es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
     607! ori          endif
     608! ori          qg=eps*es/(p(i,icb(i))-es*(1.-eps))
    542609    qg = eps*es/(p(i,icbs(i))-es*(1.-eps))
    543610
    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)
     611! Second iteration.
     612
     613
     614! ori          s=cpd+alv*alv*qg/(rrv*ticb(i)*ticb(i))
     615! ori          s=1./s
     616! ori          ahg=cpd*tg+(cl-cpd)*qnk(i)*ticb(i)+alv*qg+gzicb(i)
    550617    ahg = cpd*tg + (cl-cpd)*qnk(i)*tg + alv*qg + gzicb(i) ! convect3
    551618    tg = tg + s*(ah0(i)-ahg)
    552     ! ori          tg=max(tg,35.0)
    553     ! debug          tc=tg-t0
     619! ori          tg=max(tg,35.0)
     620! debug          tc=tg-t0
    554621    tc = tg - 273.15
    555622    denom = 243.5 + tc
    556     denom = max(denom, 1.0) ! convect3
    557     ! ori          if(tc.ge.0.0)then
     623    denom = max(denom, 1.0)                               ! convect3
     624! ori          if(tc.ge.0.0)then
    558625    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))
     626! ori          else
     627! ori           es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
     628! ori          end if
     629! ori          qg=eps*es/(p(i,icb(i))-es*(1.-eps))
    563630    qg = eps*es/(p(i,icbs(i))-es*(1.-eps))
    564631
    565632    alv = lv0 - clmcpv*(ticb(i)-273.15)
    566633
    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:
     634! ori c approximation here:
     635! ori         tp(i,icb(i))=(ah0(i)-(cl-cpd)*qnk(i)*ticb(i)
     636! ori     &   -gz(i,icb(i))-alv*qg)/cpd
     637
     638! convect3: no approximation:
    572639    tp(i, icbs(i)) = (ah0(i)-gz(i,icbs(i))-alv*qg)/(cpd+(cl-cpd)*qnk(i))
    573640
    574     ! ori         clw(i,icb(i))=qnk(i)-qg
    575     ! ori         clw(i,icb(i))=max(0.0,clw(i,icb(i)))
     641! ori         clw(i,icb(i))=qnk(i)-qg
     642! ori         clw(i,icb(i))=max(0.0,clw(i,icb(i)))
    576643    clw(i, icbs(i)) = qnk(i) - qg
    577644    clw(i, icbs(i)) = max(0.0, clw(i,icbs(i)))
    578645
    579646    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)
     647! ori         tvp(i,icb(i))=tp(i,icb(i))*(1.+rg*epsi)
     648! convect3: (qg utilise au lieu du vrai mixing ratio rg)
     649    tvp(i, icbs(i)) = tp(i, icbs(i))*(1.+qg/eps-qnk(i))   !whole thing
     650
     651  END DO
     652
     653! ori      do 380 k=minorig,icbsmax2
     654! ori       do 370 i=1,len
     655! ori         tvp(i,k)=tvp(i,k)-tp(i,k)*qnk(i)
     656! ori 370   continue
     657! ori 380  continue
     658
     659
     660! -- The following is only for convect3:
     661
     662! * icbs is the first level above the LCL:
     663! if plcl<p(icb), then icbs=icb+1
     664! if plcl>p(icb), then icbs=icb
     665
     666! * the routine above computes tvp from minorig to icbs (included).
     667
     668! * to compute buoybase (in cv3_trigger.F), both tvp(icb) and tvp(icb+1)
     669! must be known. This is the case if icbs=icb+1, but not if icbs=icb.
     670
     671! * therefore, in the case icbs=icb, we compute tvp at level icb+1
     672! (tvp at other levels will be computed in cv3_undilute2.F)
    606673
    607674
    … …  
    615682    tg = ticb(i)
    616683    qg = qsicb(i) ! convect3
    617     ! debug         alv=lv0-clmcpv*(ticb(i)-t0)
     684! debug         alv=lv0-clmcpv*(ticb(i)-t0)
    618685    alv = lv0 - clmcpv*(ticb(i)-273.15)
    619686
    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
     687! First iteration.
     688
     689! ori          s=cpd+alv*alv*qg/(rrv*ticb(i)*ticb(i))
     690    s = cpd*(1.-qnk(i)) + cl*qnk(i) &                         ! convect3
     691      +alv*alv*qg/(rrv*ticb(i)*ticb(i))                       ! convect3
    625692    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
     693! ori          ahg=cpd*tg+(cl-cpd)*qnk(i)*ticb(i)+alv*qg+gzicb(i)
     694    ahg = cpd*tg + (cl-cpd)*qnk(i)*tg + alv*qg + gzicb(i)     ! convect3
    628695    tg = tg + s*(ah0(i)-ahg)
    629     ! ori          tg=max(tg,35.0)
    630     ! debug          tc=tg-t0
     696! ori          tg=max(tg,35.0)
     697! debug          tc=tg-t0
    631698    tc = tg - 273.15
    632699    denom = 243.5 + tc
    633     denom = max(denom, 1.0) ! convect3
    634     ! ori          if(tc.ge.0.0)then
     700    denom = max(denom, 1.0)                                   ! convect3
     701! ori          if(tc.ge.0.0)then
    635702    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))
     703! ori          else
     704! ori           es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
     705! ori          endif
     706! ori          qg=eps*es/(p(i,icb(i))-es*(1.-eps))
    640707    qg = eps*es/(p(i,icb(i)+1)-es*(1.-eps))
    641708
    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
     709! Second iteration.
     710
     711
     712! ori          s=cpd+alv*alv*qg/(rrv*ticb(i)*ticb(i))
     713! ori          s=1./s
     714! ori          ahg=cpd*tg+(cl-cpd)*qnk(i)*ticb(i)+alv*qg+gzicb(i)
     715    ahg = cpd*tg + (cl-cpd)*qnk(i)*tg + alv*qg + gzicb(i)     ! convect3
    649716    tg = tg + s*(ah0(i)-ahg)
    650     ! ori          tg=max(tg,35.0)
    651     ! debug          tc=tg-t0
     717! ori          tg=max(tg,35.0)
     718! debug          tc=tg-t0
    652719    tc = tg - 273.15
    653720    denom = 243.5 + tc
    654     denom = max(denom, 1.0) ! convect3
    655     ! ori          if(tc.ge.0.0)then
     721    denom = max(denom, 1.0)                                   ! convect3
     722! ori          if(tc.ge.0.0)then
    656723    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))
     724! ori          else
     725! ori           es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
     726! ori          end if
     727! ori          qg=eps*es/(p(i,icb(i))-es*(1.-eps))
    661728    qg = eps*es/(p(i,icb(i)+1)-es*(1.-eps))
    662729
    663730    alv = lv0 - clmcpv*(ticb(i)-273.15)
    664731
    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:
     732! ori c approximation here:
     733! ori         tp(i,icb(i))=(ah0(i)-(cl-cpd)*qnk(i)*ticb(i)
     734! ori     &   -gz(i,icb(i))-alv*qg)/cpd
     735
     736! convect3: no approximation:
    670737    tp(i, icb(i)+1) = (ah0(i)-gz(i,icb(i)+1)-alv*qg)/(cpd+(cl-cpd)*qnk(i))
    671738
    672     ! ori         clw(i,icb(i))=qnk(i)-qg
    673     ! ori         clw(i,icb(i))=max(0.0,clw(i,icb(i)))
     739! ori         clw(i,icb(i))=qnk(i)-qg
     740! ori         clw(i,icb(i))=max(0.0,clw(i,icb(i)))
    674741    clw(i, icb(i)+1) = qnk(i) - qg
    675742    clw(i, icb(i)+1) = max(0.0, clw(i,icb(i)+1))
    676743
    677744    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
     745! ori         tvp(i,icb(i))=tp(i,icb(i))*(1.+rg*epsi)
     746! convect3: (qg utilise au lieu du vrai mixing ratio rg)
     747    tvp(i, icb(i)+1) = tp(i, icb(i)+1)*(1.+qg/eps-qnk(i))     !whole thing
    681748
    682749  END DO
    … …  
    685752END SUBROUTINE cv3_undilute1
    686753
    687 SUBROUTINE cv3_trigger(len, nd, icb, plcl, p, th, tv, tvp, thnk, pbase, &
    688     buoybase, iflag, sig, w0)
     754SUBROUTINE cv3_trigger(len, nd, icb, plcl, p, th, tv, tvp, thnk, &
     755                       pbase, buoybase, iflag, sig, w0)
    689756  IMPLICIT NONE
    690757
    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   ! -------------------------------------------------------------------
     758! -------------------------------------------------------------------
     759! --- TRIGGERING
     760
     761! - computes the cloud base
     762! - triggering (crude in this version)
     763! - relaxation of sig and w0 when no convection
     764
     765! Caution1: if no convection, we set iflag=4
     766! (it used to be 0 in convect3)
     767
     768! Caution2: at this stage, tvp (and thus buoy) are know up
     769! through icb only!
     770! -> the buoyancy below cloud base not (yet) set to the cloud base buoyancy
     771! -------------------------------------------------------------------
    705772
    706773  include "cv3param.h"
    707774
    708   ! input:
     775! input:
    709776  INTEGER len, nd
    710777  INTEGER icb(len)
    … …  
    713780  REAL thnk(len)
    714781
    715   ! output:
     782! output:
    716783  REAL pbase(len), buoybase(len)
    717784
    718   ! input AND output:
     785! input AND output:
    719786  INTEGER iflag(len)
    720787  REAL sig(len, nd), w0(len, nd)
    721788
    722   ! local variables:
     789! local variables:
    723790  INTEGER i, k
    724791  REAL tvpbase, tvbase, tdif, ath, ath1
    725792
    726793
    727   ! ***   set cloud base buoyancy at (plcl+dpbase) level buoyancy
     794! ***   set cloud base buoyancy at (plcl+dpbase) level buoyancy
    728795
    729796  DO i = 1, len
    730797    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))
     798    tvpbase = tvp(i, icb(i))  *(pbase(i)-p(i,icb(i)+1))/(p(i,icb(i))-p(i,icb(i)+1)) + &
     799              tvp(i, icb(i)+1)*(p(i,icb(i))-pbase(i))  /(p(i,icb(i))-p(i,icb(i)+1))
     800    tvbase = tv(i, icb(i))  *(pbase(i)-p(i,icb(i)+1))/(p(i,icb(i))-p(i,icb(i)+1)) + &
     801             tv(i, icb(i)+1)*(p(i,icb(i))-pbase(i))  /(p(i,icb(i))-p(i,icb(i)+1))
    737802    buoybase(i) = tvpbase - tvbase
    738803  END DO
    739804
    740805
    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)
     806! ***   make sure that column is dry adiabatic between the surface  ***
     807! ***    and cloud base, and that lifted air is positively buoyant  ***
     808! ***                         at cloud base                         ***
     809! ***       if not, return to calling program after resetting       ***
     810! ***                        sig(i) and w0(i)                       ***
     811
     812
     813! oct3      do 200 i=1,len
     814! oct3
     815! oct3       tdif = buoybase(i)
     816! oct3       ath1 = th(i,1)
     817! oct3       ath  = th(i,icb(i)-1) - dttrig
     818! oct3
     819! oct3       if (tdif.lt.dtcrit .or. ath.gt.ath1) then
     820! oct3         do 60 k=1,nl
     821! oct3            sig(i,k) = beta*sig(i,k) - 2.*alpha*tdif*tdif
     822! oct3            sig(i,k) = AMAX1(sig(i,k),0.0)
     823! oct3            w0(i,k)  = beta*w0(i,k)
     824! oct3   60    continue
     825! oct3         iflag(i)=4 ! pour version vectorisee
     826! oct3c convect3         iflag(i)=0
     827! oct3cccc         return
     828! oct3       endif
     829! oct3
     830! oct3200   continue
     831
     832! -- oct3: on reecrit la boucle 200 (pour la vectorisation)
    768833
    769834  DO k = 1, nl
    … …  
    779844        w0(i, k) = beta*w0(i, k)
    780845        iflag(i) = 4 ! pour version vectorisee
    781         ! convect3         iflag(i)=0
     846! convect3         iflag(i)=0
    782847      END IF
    783848
    … …  
    785850  END DO
    786851
    787   ! fin oct3 --
     852! fin oct3 --
    788853
    789854  RETURN
    790855END SUBROUTINE cv3_trigger
    791856
    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)
     857SUBROUTINE cv3_compress(len, nloc, ncum, nd, ntra, &
     858                        iflag1, nk1, icb1, icbs1, &
     859                        plcl1, tnk1, qnk1, gznk1, pbase1, buoybase1, &
     860                        t1, q1, qs1, u1, v1, gz1, th1, &
     861                        tra1, &
     862                        h1, lv1, cpn1, p1, ph1, tv1, tp1, tvp1, clw1, &
     863                        sig1, w01, &
     864                        iflag, nk, icb, icbs, &
     865                        plcl, tnk, qnk, gznk, pbase, buoybase, &
     866                        t, q, qs, u, v, gz, th, &
     867                        tra, &
     868                        h, lv, cpn, p, ph, tv, tp, tvp, clw, &
     869                        sig, w0)
    797870  IMPLICIT NONE
    798871
    … …  
    800873  include 'iniprint.h'
    801874
    802   ! inputs:
     875!inputs:
    803876  INTEGER len, ncum, nd, ntra, nloc
    804877  INTEGER iflag1(len), nk1(len), icb1(len), icbs1(len)
    … …  
    813886  REAL tra1(len, nd, ntra)
    814887
    815   ! outputs:
    816   ! en fait, on a nloc=len pour l'instant (cf cv_driver)
     888!outputs:
     889! en fait, on a nloc=len pour l'instant (cf cv_driver)
    817890  INTEGER iflag(nloc), nk(nloc), icb(nloc), icbs(nloc)
    818891  REAL plcl(nloc), tnk(nloc), qnk(nloc), gznk(nloc)
    … …  
    826899  REAL tra(nloc, nd, ntra)
    827900
    828   ! local variables:
     901!local variables:
    829902  INTEGER i, k, nn, j
    830903
    … …  
    859932  END DO
    860933
    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
     934!AC!      do 121 j=1,ntra
     935!AC!ccccc      do 111 k=1,nl+1
     936!AC!      do 111 k=1,nd
     937!AC!       nn=0
     938!AC!      do 101 i=1,len
     939!AC!      if(iflag1(i).eq.0)then
     940!AC!       nn=nn+1
     941!AC!       tra(nn,k,j)=tra1(i,k,j)
     942!AC!      endif
     943!AC! 101  continue
     944!AC! 111  continue
     945!AC! 121  continue
    873946
    874947  IF (nn/=ncum) THEN
    … …  
    902975
    903976
    904   ! JAM--------------------------------------------------------------------
    905   ! Calcul de la quantité d'eau sous forme de glace
    906   ! --------------------------------------------------------------------
     977!JAM--------------------------------------------------------------------
     978! Calcul de la quantité d'eau sous forme de glace
     979! --------------------------------------------------------------------
    907980  REAL qi(len, nl)
    908981  REAL t(len, nl), clw(len, nl)
    … …  
    922995        END IF
    923996      END IF
    924       ! print*,t(i,k),qi(i,k),'temp,testglace'
     997! print*,t(i,k),qi(i,k),'temp,testglace'
    925998    END DO
    926999  END DO
    … …  
    9301003END SUBROUTINE icefrac
    9311004
    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)
     1005SUBROUTINE cv3_undilute2(nloc, ncum, nd, icb, icbs, nk, &
     1006                         tnk, qnk, gznk, hnk, t, q, qs, gz, &
     1007                         p, h, tv, lv, lf, pbase, buoybase, plcl, &
     1008                         inb, tp, tvp, clw, hp, ep, sigp, buoy, frac)
    9351009  IMPLICIT NONE
    9361010
    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   ! ---------------------------------------------------------------------
     1011! ---------------------------------------------------------------------
     1012! Purpose:
     1013! FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
     1014! &
     1015! COMPUTE THE PRECIPITATION EFFICIENCIES AND THE
     1016! FRACTION OF PRECIPITATION FALLING OUTSIDE OF CLOUD
     1017! &
     1018! FIND THE LEVEL OF NEUTRAL BUOYANCY
     1019
     1020! Main differences convect3/convect4:
     1021!  - icbs (input) is the first level above LCL (may differ from icb)
     1022!  - many minor differences in the iterations
     1023!  - condensed water not removed from tvp in convect3
     1024!  - vertical profile of buoyancy computed here (use of buoybase)
     1025!  - the determination of inb is different
     1026!  - no inb1, only inb in output
     1027! ---------------------------------------------------------------------
    9541028
    9551029  include "cvthermo.h"
    … …  
    9581032  include "cvflag.h"
    9591033
    960   ! inputs:
     1034!inputs:
    9611035  INTEGER ncum, nd, nloc, j
    9621036  INTEGER icb(nloc), icbs(nloc), nk(nloc)
    … …  
    9681042  REAL pbase(nloc), buoybase(nloc), plcl(nloc)
    9691043
    970   ! outputs:
     1044!outputs:
    9711045  INTEGER inb(nloc)
    9721046  REAL tp(nloc, nd), tvp(nloc, nd), clw(nloc, nd)
    … …  
    9741048  REAL buoy(nloc, nd)
    9751049
    976   ! local variables:
     1050!local variables:
    9771051  INTEGER i, k
    9781052  REAL tg, qg, ahg, alv, alf, s, tc, es, esi, denom, rg, tca, elacrit
    … …  
    9861060  REAL fracg
    9871061
    988   ! =====================================================================
    989   ! --- SOME INITIALIZATIONS
    990   ! =====================================================================
     1062! =====================================================================
     1063! --- SOME INITIALIZATIONS
     1064! =====================================================================
    9911065
    9921066  DO k = 1, nl
    … …  
    9981072  END DO
    9991073
    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   ***
     1074! =====================================================================
     1075! --- FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
     1076! =====================================================================
     1077
     1078! ---       The procedure is to solve the equation.
     1079!                cp*tp+L*qp+phi=cp*tnk+L*qnk+gznk.
     1080
     1081! ***  Calculate certain parcel quantities, including static energy   ***
    10081082
    10091083
    10101084  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    ***
     1085    ah0(i) = (cpd*(1.-qnk(i))+cl*qnk(i))*tnk(i)+ &
     1086! debug          qnk(i)*(lv0-clmcpv*(tnk(i)-t0))+gznk(i)
     1087             qnk(i)*(lv0-clmcpv*(tnk(i)-273.15)) + gznk(i)
     1088  END DO
     1089
     1090
     1091! ***  Find lifted parcel quantities above cloud base    ***
    10181092
    10191093
    10201094  DO k = minorig + 1, nl
    10211095    DO i = 1, ncum
    1022       ! ori         if(k.ge.(icb(i)+1))then
    1023       IF (k>=(icbs(i)+1)) THEN ! convect3
     1096! ori       if(k.ge.(icb(i)+1))then
     1097      IF (k>=(icbs(i)+1)) THEN                                ! convect3
    10241098        tg = t(i, k)
    10251099        qg = qs(i, k)
    1026         ! debug       alv=lv0-clmcpv*(t(i,k)-t0)
     1100! debug       alv=lv0-clmcpv*(t(i,k)-t0)
    10271101        alv = lv0 - clmcpv*(t(i,k)-273.15)
    10281102
    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
     1103! First iteration.
     1104
     1105! ori          s=cpd+alv*alv*qg/(rrv*t(i,k)*t(i,k))
     1106        s = cpd*(1.-qnk(i)) + cl*qnk(i) + &                  ! convect3
     1107            alv*alv*qg/(rrv*t(i,k)*t(i,k))                    ! convect3
    10341108        s = 1./s
    1035         ! ori          ahg=cpd*tg+(cl-cpd)*qnk(i)*t(i,k)+alv*qg+gz(i,k)
     1109! ori          ahg=cpd*tg+(cl-cpd)*qnk(i)*t(i,k)+alv*qg+gz(i,k)
    10361110        ahg = cpd*tg + (cl-cpd)*qnk(i)*tg + alv*qg + gz(i, k) ! convect3
    10371111        tg = tg + s*(ah0(i)-ahg)
    1038         ! ori          tg=max(tg,35.0)
    1039         ! debug        tc=tg-t0
     1112! ori          tg=max(tg,35.0)
     1113! debug        tc=tg-t0
    10401114        tc = tg - 273.15
    10411115        denom = 243.5 + tc
    1042         denom = max(denom, 1.0) ! convect3
    1043         ! ori          if(tc.ge.0.0)then
     1116        denom = max(denom, 1.0)                               ! convect3
     1117! ori          if(tc.ge.0.0)then
    10441118        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
     1119! ori          else
     1120! ori                   es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
     1121! ori          endif
    10481122        qg = eps*es/(p(i,k)-es*(1.-eps))
    10491123
    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)
     1124! Second iteration.
     1125
     1126! ori          s=cpd+alv*alv*qg/(rrv*t(i,k)*t(i,k))
     1127! ori          s=1./s
     1128! ori          ahg=cpd*tg+(cl-cpd)*qnk(i)*t(i,k)+alv*qg+gz(i,k)
    10551129        ahg = cpd*tg + (cl-cpd)*qnk(i)*tg + alv*qg + gz(i, k) ! convect3
    10561130        tg = tg + s*(ah0(i)-ahg)
    1057         ! ori          tg=max(tg,35.0)
    1058         ! debug        tc=tg-t0
     1131! ori          tg=max(tg,35.0)
     1132! debug        tc=tg-t0
    10591133        tc = tg - 273.15
    10601134        denom = 243.5 + tc
    1061         denom = max(denom, 1.0) ! convect3
    1062         ! ori          if(tc.ge.0.0)then
     1135        denom = max(denom, 1.0)                               ! convect3
     1136! ori          if(tc.ge.0.0)then
    10631137        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
     1138! ori          else
     1139! ori                   es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
     1140! ori          endif
    10671141        qg = eps*es/(p(i,k)-es*(1.-eps))
    10681142
    1069         ! debug        alv=lv0-clmcpv*(t(i,k)-t0)
     1143! debug        alv=lv0-clmcpv*(t(i,k)-t0)
    10701144        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:
     1145! print*,'cpd dans convect2 ',cpd
     1146! print*,'tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd'
     1147! print*,tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd
     1148
     1149! ori c approximation here:
     1150! ori        tp(i,k)=(ah0(i)-(cl-cpd)*qnk(i)*t(i,k)-gz(i,k)-alv*qg)/cpd
     1151
     1152! convect3: no approximation:
    10801153        IF (cvflag_ice) THEN
    10811154          tp(i, k) = max(0., (ah0(i)-gz(i,k)-alv*qg)/(cpd+(cl-cpd)*qnk(i)))
    … …  
    10871160        clw(i, k) = max(0.0, clw(i,k))
    10881161        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):
     1162! ori               tvp(i,k)=tp(i,k)*(1.+rg*epsi)
     1163! convect3: (qg utilise au lieu du vrai mixing ratio rg):
    10911164        tvp(i, k) = tp(i, k)*(1.+qg/eps-qnk(i)) ! whole thing
    10921165        IF (cvflag_ice) THEN
    … …  
    10991172
    11001173      IF (cvflag_ice) THEN
    1101         ! CR:attention boucle en klon dans Icefrac
    1102         ! Call Icefrac(t,clw,qi,nl,nloc)
     1174!CR:attention boucle en klon dans Icefrac
     1175! Call Icefrac(t,clw,qi,nl,nloc)
    11031176        IF (t(i,k)>263.15) THEN
    11041177          qi(i, k) = 0.
    … …  
    11111184          END IF
    11121185        END IF
    1113         ! CR: fin test
     1186!CR: fin test
    11141187        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
     1188!CR: on commente les calculs d'Arnaud car division par zero
     1189! nouveau calcul propose par JYG
     1190!      alv=lv0-clmcpv*(t(i,k)-273.15)
     1191!      alf=lf0-clmci*(t(i,k)-273.15)
     1192!      tg=tp(i,k)
     1193!      tc=tp(i,k)-273.15
     1194!      denom=243.5+tc
     1195!      do j=1,3
     1196! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
     1197! il faudra que esi vienne en argument de la convection
     1198! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
     1199!        tbis=t(i,k)+(tp(i,k)-tg)
     1200!        esi=exp(23.33086-(6111.72784/tbis) + &
     1201!                       0.15215*log(tbis))
     1202!        qsat_new=eps*esi/(p(i,k)-esi*(1.-eps))
     1203!        snew=cpd*(1.-qnk(i))+cl*qnk(i)+alv*alv*qsat_new/ &
     1204!                                       (rrv*tbis*tbis)
     1205!        snew=1./snew
     1206!        print*,esi,qsat_new,snew,'esi,qsat,snew'
     1207!        tp(i,k)=tg+(alf*qi(i,k)+alv*qg*(1.-(esi/es)))*snew
     1208!        print*,k,tp(i,k),qnk(i),'avec glace'
     1209!        print*,'tpNAN',tg,alf,qi(i,k),alv,qg,esi,es,snew
     1210!      enddo
    11381211
    11391212          alv = lv0 - clmcpv*(t(i,k)-273.15)
    … …  
    11451218            esi = exp(23.33086-(6111.72784/tp(i,k))+0.15215*log(tp(i,k)))
    11461219            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))
     1220            snew = cpd*(1.-qnk(i)) + cl*qnk(i) + alv*als*qsat_new/ &
     1221                                                 (rrv*tp(i,k)*tp(i,k))
    11491222            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'
     1223! c             print*,esi,qsat_new,snew,'esi,qsat,snew'
     1224            tp(i, k) = tp(i, k) + &
     1225                       ((cpd*(1.-qnk(i))+cl*qnk(i))*(tg-tp(i,k)) + &
     1226                        alv*(qg-qsat_new)+alf*qi(i,k))*snew
     1227! print*,k,tp(i,k),qsat_new,qnk(i),qi(i,k), &
     1228!              'k,tp,q,qt,qi avec glace'
    11551229          END DO
    11561230
    1157           ! CR:reprise du code AJ
     1231!CR:reprise du code AJ
    11581232          clw(i, k) = qnk(i) - qsat_new
    11591233          clw(i, k) = max(0.0, clw(i,k))
    11601234          tvp(i, k) = max(0., tp(i,k)*(1.+qsat_new/eps-qnk(i)))
    1161           ! print*,tvp(i,k),'tvp'
     1235! print*,tvp(i,k),'tvp'
    11621236        END IF
    11631237        IF (clw(i,k)<1.E-11) THEN
    … …  
    11701244  END DO
    11711245
    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   ! =====================================================================
     1246! =====================================================================
     1247! --- SET THE PRECIPITATION EFFICIENCIES AND THE FRACTION OF
     1248! --- PRECIPITATION FALLING OUTSIDE OF CLOUD
     1249! --- THESE MAY BE FUNCTIONS OF TP(I), P(I) AND CLW(I)
     1250! =====================================================================
    11771251
    11781252  IF (flag_epkeorig/=1) THEN
    … …  
    12051279    END DO
    12061280  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:
     1281! =====================================================================
     1282! --- CALCULATE VIRTUAL TEMPERATURE AND LIFTED PARCEL
     1283! --- VIRTUAL TEMPERATURE
     1284! =====================================================================
     1285
     1286! dans convect3, tvp est calcule en une seule fois, et sans retirer
     1287! l'eau condensee (~> reversible CAPE)
     1288
     1289! ori      do 340 k=minorig+1,nl
     1290! ori        do 330 i=1,ncum
     1291! ori        if(k.ge.(icb(i)+1))then
     1292! ori          tvp(i,k)=tvp(i,k)*(1.0-qnk(i)+ep(i,k)*clw(i,k))
     1293! oric         print*,'i,k,tvp(i,k),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! ori        endif
     1296! ori 330    continue
     1297! ori 340  continue
     1298
     1299! ori      do 350 i=1,ncum
     1300! ori       tvp(i,nlp)=tvp(i,nl)-(gz(i,nlp)-gz(i,nl))/cpd
     1301! ori 350  continue
     1302
     1303  DO i = 1, ncum                                           ! convect3
     1304    tp(i, nlp) = tp(i, nl)                                 ! convect3
     1305  END DO                                                   ! convect3
     1306
     1307! =====================================================================
     1308! --- EFFECTIVE VERTICAL PROFILE OF BUOYANCY (convect3 only):
     1309! =====================================================================
     1310
     1311! -- this is for convect3 only:
     1312
     1313! first estimate of buoyancy:
    12401314
    12411315  DO i = 1, ncum
    … …  
    12451319  END DO
    12461320
    1247   ! set buoyancy=buoybase for all levels below base
    1248   ! for safety, set buoy(icb)=buoybase
     1321! set buoyancy=buoybase for all levels below base
     1322! for safety, set buoy(icb)=buoybase
    12491323
    12501324  DO i = 1, ncum
    … …  
    12541328      END IF
    12551329    END DO
    1256     ! buoy(icb(i),k)=buoybase(i)
     1330!    buoy(icb(i),k)=buoybase(i)
    12571331    buoy(i, icb(i)) = buoybase(i)
    12581332  END DO
    12591333
    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:
     1334! -- end convect3
     1335
     1336! =====================================================================
     1337! --- FIND THE FIRST MODEL LEVEL (INB) ABOVE THE PARCEL'S
     1338! --- LEVEL OF NEUTRAL BUOYANCY
     1339! =====================================================================
     1340
     1341! -- this is for convect3 only:
    12681342
    12691343  DO i = 1, ncum
    … …  
    12731347
    12741348
    1275   ! --    iposit(i) = first level, above icb, with positive buoyancy
     1349! --    iposit(i) = first level, above icb, with positive buoyancy
    12761350  DO k = 1, nl - 1
    12771351    DO i = 1, ncum
    … …  
    12961370  END DO
    12971371
    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   ! =====================================================================
     1372! -- end convect3
     1373
     1374! ori      do 510 i=1,ncum
     1375! ori        cape(i)=0.0
     1376! ori        capem(i)=0.0
     1377! ori        inb(i)=icb(i)+1
     1378! ori        inb1(i)=inb(i)
     1379! ori 510  continue
     1380
     1381! Originial Code
     1382
     1383!    do 530 k=minorig+1,nl-1
     1384!    do 520 i=1,ncum
     1385!      if(k.ge.(icb(i)+1))then
     1386!      by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
     1387!      byp=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
     1388!      cape(i)=cape(i)+by
     1389!      if(by.ge.0.0)inb1(i)=k+1
     1390!      if(cape(i).gt.0.0)then
     1391!        inb(i)=k+1
     1392!        capem(i)=cape(i)
     1393!      endif
     1394!      endif
     1395!520    continue
     1396!530  continue
     1397!    do 540 i=1,ncum
     1398!    byp=(tvp(i,nl)-tv(i,nl))*dph(i,nl)/p(i,nl)
     1399!    cape(i)=capem(i)+byp
     1400!    defrac=capem(i)-cape(i)
     1401!    defrac=max(defrac,0.001)
     1402!    frac(i)=-cape(i)/defrac
     1403!    frac(i)=min(frac(i),1.0)
     1404!    frac(i)=max(frac(i),0.0)
     1405!540   continue
     1406
     1407!    K Emanuel fix
     1408
     1409!    call zilch(byp,ncum)
     1410!    do 530 k=minorig+1,nl-1
     1411!    do 520 i=1,ncum
     1412!      if(k.ge.(icb(i)+1))then
     1413!      by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
     1414!      cape(i)=cape(i)+by
     1415!      if(by.ge.0.0)inb1(i)=k+1
     1416!      if(cape(i).gt.0.0)then
     1417!        inb(i)=k+1
     1418!        capem(i)=cape(i)
     1419!        byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
     1420!      endif
     1421!      endif
     1422!520    continue
     1423!530  continue
     1424!    do 540 i=1,ncum
     1425!    inb(i)=max(inb(i),inb1(i))
     1426!    cape(i)=capem(i)+byp(i)
     1427!    defrac=capem(i)-cape(i)
     1428!    defrac=max(defrac,0.001)
     1429!    frac(i)=-cape(i)/defrac
     1430!    frac(i)=min(frac(i),1.0)
     1431!    frac(i)=max(frac(i),0.0)
     1432!540   continue
     1433
     1434! J Teixeira fix
     1435
     1436! ori      call zilch(byp,ncum)
     1437! ori      do 515 i=1,ncum
     1438! ori        lcape(i)=.true.
     1439! ori 515  continue
     1440! ori      do 530 k=minorig+1,nl-1
     1441! ori        do 520 i=1,ncum
     1442! ori          if(cape(i).lt.0.0)lcape(i)=.false.
     1443! ori          if((k.ge.(icb(i)+1)).and.lcape(i))then
     1444! ori            by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
     1445! ori            byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
     1446! ori            cape(i)=cape(i)+by
     1447! ori            if(by.ge.0.0)inb1(i)=k+1
     1448! ori            if(cape(i).gt.0.0)then
     1449! ori              inb(i)=k+1
     1450! ori              capem(i)=cape(i)
     1451! ori            endif
     1452! ori          endif
     1453! ori 520    continue
     1454! ori 530  continue
     1455! ori      do 540 i=1,ncum
     1456! ori          cape(i)=capem(i)+byp(i)
     1457! ori          defrac=capem(i)-cape(i)
     1458! ori          defrac=max(defrac,0.001)
     1459! ori          frac(i)=-cape(i)/defrac
     1460! ori          frac(i)=min(frac(i),1.0)
     1461! ori          frac(i)=max(frac(i),0.0)
     1462! ori 540  continue
     1463
     1464! =====================================================================
     1465! ---   CALCULATE LIQUID WATER STATIC ENERGY OF LIFTED PARCEL
     1466! =====================================================================
    13931467
    13941468  DO k = 1, nd
    … …  
    14051479          frac(i, k) = 1. - (t(i,k)-243.15)/(263.15-243.15)
    14061480          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)
     1481          hp(i, k) = hnk(i) + (lv(i,k)+(cpd-cpv)*t(i,k)+frac(i,k)*lf(i,k))* &
     1482                              ep(i, k)*clw(i, k)
    14091483
    14101484        ELSE
    … …  
    14191493END SUBROUTINE cv3_undilute2
    14201494
    1421 SUBROUTINE cv3_closure(nloc, ncum, nd, icb, inb, pbase, p, ph, tv, buoy, sig, &
    1422     w0, cape, m, iflag)
     1495SUBROUTINE cv3_closure(nloc, ncum, nd, icb, inb, &
     1496                       pbase, p, ph, tv, buoy, &
     1497                       sig, w0, cape, m, iflag)
    14231498  IMPLICIT NONE
    14241499
    1425   ! ===================================================================
    1426   ! ---  CLOSURE OF CONVECT3
    1427 
    1428   ! vectorization: S. Bony
    1429   ! ===================================================================
     1500! ===================================================================
     1501! ---  CLOSURE OF CONVECT3
     1502!
     1503! vectorization: S. Bony
     1504! ===================================================================
    14301505
    14311506  include "cvthermo.h"
    14321507  include "cv3param.h"
    14331508
    1434   ! input:
     1509!input:
    14351510  INTEGER ncum, nd, nloc
    14361511  INTEGER icb(nloc), inb(nloc)
    … …  
    14391514  REAL tv(nloc, nd), buoy(nloc, nd)
    14401515
    1441   ! input/output:
     1516!input/output:
    14421517  REAL sig(nloc, nd), w0(nloc, nd)
    14431518  INTEGER iflag(nloc)
    14441519
    1445   ! output:
     1520!output:
    14461521  REAL cape(nloc)
    14471522  REAL m(nloc, nd)
    14481523
    1449   ! local variables:
     1524!local variables:
    14501525  INTEGER i, j, k, icbmax
    14511526  REAL deltap, fac, w, amu
    … …  
    14541529
    14551530
    1456   ! -------------------------------------------------------
    1457   ! -- Initialization
    1458   ! -------------------------------------------------------
     1531! -------------------------------------------------------
     1532! -- Initialization
     1533! -------------------------------------------------------
    14591534
    14601535  DO k = 1, nl
    … …  
    14641539  END DO
    14651540
    1466   ! -------------------------------------------------------
    1467   ! -- Reset sig(i) and w0(i) for i>inb and i<icb
    1468   ! -------------------------------------------------------
    1469 
    1470   ! update sig and w0 above LNB:
     1541! -------------------------------------------------------
     1542! -- Reset sig(i) and w0(i) for i>inb and i<icb
     1543! -------------------------------------------------------
     1544
     1545! update sig and w0 above LNB:
    14711546
    14721547  DO k = 1, nl - 1
    14731548    DO i = 1, ncum
    14741549      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)))
     1550        sig(i, k) = beta*sig(i, k) + &
     1551                    2.*alpha*buoy(i, inb(i))*abs(buoy(i,inb(i)))
    14771552        sig(i, k) = amax1(sig(i,k), 0.0)
    14781553        w0(i, k) = beta*w0(i, k)
    … …  
    14811556  END DO
    14821557
    1483   ! compute icbmax:
     1558! compute icbmax:
    14841559
    14851560  icbmax = 2
    … …  
    14881563  END DO
    14891564
    1490   ! update sig and w0 below cloud base:
     1565! update sig and w0 below cloud base:
    14911566
    14921567  DO k = 1, icbmax
    14931568    DO i = 1, ncum
    14941569      IF (k<=icb(i)) THEN
    1495         sig(i, k) = beta*sig(i, k) - 2.*alpha*buoy(i, icb(i))*buoy(i, icb(i))
     1570        sig(i, k) = beta*sig(i, k) - &
     1571                    2.*alpha*buoy(i, icb(i))*buoy(i, icb(i))
    14961572        sig(i, k) = max(sig(i,k), 0.0)
    14971573        w0(i, k) = beta*w0(i, k)
    … …  
    15001576  END DO
    15011577
    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   ! -------------------------------------------------------------
     1578!!      if(inb.lt.(nl-1))then
     1579!!         do 85 i=inb+1,nl-1
     1580!!            sig(i)=beta*sig(i)+2.*alpha*buoy(inb)*
     1581!!     1              abs(buoy(inb))
     1582!!            sig(i)=max(sig(i),0.0)
     1583!!            w0(i)=beta*w0(i)
     1584!!   85    continue
     1585!!      end if
     1586
     1587!!      do 87 i=1,icb
     1588!!         sig(i)=beta*sig(i)-2.*alpha*buoy(icb)*buoy(icb)
     1589!!         sig(i)=max(sig(i),0.0)
     1590!!         w0(i)=beta*w0(i)
     1591!!   87 continue
     1592
     1593! -------------------------------------------------------------
     1594! -- Reset fractional areas of updrafts and w0 at initial time
     1595! -- and after 10 time steps of no convection
     1596! -------------------------------------------------------------
    15211597
    15221598  DO k = 1, nl - 1
    … …  
    15291605  END DO
    15301606
    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   ! -------------------------------------------------------------
     1607! -------------------------------------------------------------
     1608! -- Calculate convective available potential energy (cape),
     1609! -- vertical velocity (w), fractional area covered by
     1610! -- undilute updraft (sig), and updraft mass flux (m)
     1611! -------------------------------------------------------------
    15361612
    15371613  DO i = 1, ncum
    … …  
    15391615  END DO
    15401616
    1541   ! compute dtmin (minimum buoyancy between ICB and given level k):
     1617! compute dtmin (minimum buoyancy between ICB and given level k):
    15421618
    15431619  DO i = 1, ncum
    … …  
    15501626    DO k = 1, nl
    15511627      DO j = minorig, nl
    1552         IF ((k>=(icb(i)+1)) .AND. (k<=inb(i)) .AND. (j>=icb(i)) .AND. (j<=(k- &
    1553             1))) THEN
     1628        IF ((k>=(icb(i)+1)) .AND. (k<=inb(i)) .AND. (j>=icb(i)) .AND. (j<=(k-1))) THEN
    15541629          dtmin(i, k) = amin1(dtmin(i,k), buoy(i,j))
    15551630        END IF
    … …  
    15581633  END DO
    15591634
    1560   ! the interval on which cape is computed starts at pbase :
     1635! the interval on which cape is computed starts at pbase :
    15611636
    15621637  DO k = 1, nl
    … …  
    15701645        sigold(i, k) = sig(i, k)
    15711646
    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
     1647! dtmin(i,k)=100.0
     1648! do 97 j=icb(i),k-1 ! mauvaise vectorisation
     1649! dtmin(i,k)=AMIN1(dtmin(i,k),buoy(i,j))
     1650! 97     continue
    15761651
    15771652        sig(i, k) = beta*sig(i, k) + alpha*dtmin(i, k)*abs(dtmin(i,k))
    … …  
    15901665  DO i = 1, ncum
    15911666    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))
     1667    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))
    15941668    sig(i, icb(i)) = sig(i, icb(i)+1)
    15951669    sig(i, icb(i)-1) = sig(i, icb(i))
    15961670  END DO
    15971671
    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)
     1672! ccc 3. Compute final cloud base mass flux and set iflag to 3 if
     1673! ccc    cloud base mass flux is exceedingly small and is decreasing (i.e. if
     1674! ccc    the final mass flux (cbmflast) is greater than the target mass flux
     1675! ccc    (cbmf) ??).
     1676! cc
     1677! c      do i = 1,ncum
     1678! c       cbmflast(i) = 0.
     1679! c      enddo
     1680! cc
     1681! c      do k= 1,nl
     1682! c       do i = 1,ncum
     1683! c        IF (k .ge. icb(i) .and. k .le. inb(i)) THEN
     1684! c         cbmflast(i) = cbmflast(i)+M(i,k)
     1685! c        ENDIF
     1686! c       enddo
     1687! c      enddo
     1688! cc
     1689! c      do i = 1,ncum
     1690! c       IF (cbmflast(i) .lt. 1.e-6) THEN
     1691! c         iflag(i) = 3
     1692! c       ENDIF
     1693! c      enddo
     1694! cc
     1695! c      do k= 1,nl
     1696! c       do i = 1,ncum
     1697! c        IF (iflag(i) .ge. 3) THEN
     1698! c         M(i,k) = 0.
     1699! c         sig(i,k) = 0.
     1700! c         w0(i,k) = 0.
     1701! c        ENDIF
     1702! c       enddo
     1703! c      enddo
     1704! cc
     1705!!      cape=0.0
     1706!!      do 98 i=icb+1,inb
     1707!!         deltap = min(pbase,ph(i-1))-min(pbase,ph(i))
     1708!!         cape=cape+rrd*buoy(i-1)*deltap/p(i-1)
     1709!!         dcape=rrd*buoy(i-1)*deltap/p(i-1)
     1710!!         dlnp=deltap/p(i-1)
     1711!!         cape=max(0.0,cape)
     1712!!         sigold=sig(i)
     1713
     1714!!         dtmin=100.0
     1715!!         do 97 j=icb,i-1
     1716!!            dtmin=amin1(dtmin,buoy(j))
     1717!!   97    continue
     1718
     1719!!         sig(i)=beta*sig(i)+alpha*dtmin*abs(dtmin)
     1720!!         sig(i)=max(sig(i),0.0)
     1721!!         sig(i)=amin1(sig(i),0.01)
     1722!!         fac=amin1(((dtcrit-dtmin)/dtcrit),1.0)
     1723!!         w=(1.-beta)*fac*sqrt(cape)+beta*w0(i)
     1724!!         amu=0.5*(sig(i)+sigold)*w
     1725!!         m(i)=amu*0.007*p(i)*(ph(i)-ph(i+1))/tv(i)
     1726!!         w0(i)=w
     1727!!   98 continue
     1728!!      w0(icb)=0.5*w0(icb+1)
     1729!!      m(icb)=0.5*m(icb+1)*(ph(icb)-ph(icb+1))/(ph(icb+1)-ph(icb+2))
     1730!!      sig(icb)=sig(icb+1)
     1731!!      sig(icb-1)=sig(icb)
    16601732
    16611733  RETURN
    16621734END SUBROUTINE cv3_closure
    16631735
    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)
     1736SUBROUTINE cv3_mixing(nloc, ncum, nd, na, ntra, icb, nk, inb, &
     1737                      ph, t, rr, rs, u, v, tra, h, lv, lf, frac, qnk, &
     1738                      unk, vnk, hp, tv, tvp, ep, clw, m, sig, &
     1739                      ment, qent, uent, vent, nent, sij, elij, ments, qents, traent)
    16671740  IMPLICIT NONE
    16681741
    1669   ! ---------------------------------------------------------------------
    1670   ! a faire:
    1671   ! - vectorisation de la partie normalisation des flux (do 789...)
    1672   ! ---------------------------------------------------------------------
     1742! ---------------------------------------------------------------------
     1743! a faire:
     1744! - vectorisation de la partie normalisation des flux (do 789...)
     1745! ---------------------------------------------------------------------
    16731746
    16741747  include "cvthermo.h"
    … …  
    16761749  include "cvflag.h"
    16771750
    1678   ! inputs:
     1751!inputs:
    16791752  INTEGER ncum, nd, na, ntra, nloc
    16801753  INTEGER icb(nloc), inb(nloc), nk(nloc)
    … …  
    16901763  REAL m(nloc, na) ! input of convect3
    16911764
    1692   ! outputs:
     1765!outputs:
    16931766  REAL ment(nloc, na, na), qent(nloc, na, na)
    16941767  REAL uent(nloc, na, na), vent(nloc, na, na)
    … …  
    16991772  INTEGER nent(nloc, nd)
    17001773
    1701   ! local variables:
     1774!local variables:
    17021775  INTEGER i, j, k, il, im, jm
    17031776  INTEGER num1, num2
    … …  
    17101783  LOGICAL lwork(nloc)
    17111784
    1712   ! =====================================================================
    1713   ! --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS
    1714   ! =====================================================================
    1715 
    1716   ! ori        do 360 i=1,ncum*nlp
     1785! =====================================================================
     1786! --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS
     1787! =====================================================================
     1788
     1789! ori        do 360 i=1,ncum*nlp
    17171790  DO j = 1, nl
    17181791    DO i = 1, ncum
    17191792      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
     1793! in convect3, m is computed in cv3_closure
     1794! ori          m(i,1)=0.0
     1795    END DO
     1796  END DO
     1797
     1798! ori      do 400 k=1,nlp
     1799! ori       do 390 j=1,nlp
    17271800  DO j = 1, nl
    17281801    DO k = 1, nl
    … …  
    17321805        vent(i, k, j) = v(i, j)
    17331806        elij(i, k, j) = 0.0
    1734         ! ym            ment(i,k,j)=0.0
    1735         ! ym            sij(i,k,j)=0.0
     1807!ym            ment(i,k,j)=0.0
     1808!ym            sij(i,k,j)=0.0
    17361809      END DO
    17371810    END DO
    17381811  END DO
    17391812
    1740   ! ym
     1813!ym
    17411814  ment(1:ncum, 1:nd, 1:nd) = 0.0
    17421815  sij(1:ncum, 1:nd, 1:nd) = 0.0
    17431816
    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
     1817!AC!      do k=1,ntra
     1818!AC!       do j=1,nd  ! instead nlp
     1819!AC!        do i=1,nd ! instead nlp
     1820!AC!         do il=1,ncum
     1821!AC!            traent(il,i,j,k)=tra(il,j,k)
     1822!AC!         enddo
     1823!AC!        enddo
     1824!AC!       enddo
     1825!AC!      enddo
    17531826  zm(:, :) = 0.
    17541827
    1755   ! =====================================================================
    1756   ! --- CALCULATE ENTRAINED AIR MASS FLUX (ment), TOTAL WATER MIXING
    1757   ! --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
    1758   ! --- FRACTION (sij)
    1759   ! =====================================================================
     1828! =====================================================================
     1829! --- CALCULATE ENTRAINED AIR MASS FLUX (ment), TOTAL WATER MIXING
     1830! --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
     1831! --- FRACTION (sij)
     1832! =====================================================================
    17601833
    17611834  DO i = minorig + 1, nl
    … …  
    17631836    DO j = minorig, nl
    17641837      DO il = 1, ncum
    1765         IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (j>=(icb(il)- &
    1766             1)) .AND. (j<=inb(il))) THEN
     1838        IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (j>=(icb(il)-1)) .AND. (j<=inb(il))) THEN
    17671839
    17681840          rti = qnk(il) - ep(il, i)*clw(il, i)
    … …  
    17711843
    17721844          IF (cvflag_ice) THEN
    1773             ! print*,cvflag_ice,'cvflag_ice dans do 700'
     1845! print*,cvflag_ice,'cvflag_ice dans do 700'
    17741846            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)
     1847              bf2 = 1. + (lf(il,j)+lv(il,j))*(lv(il,j)+frac(il,j)* &
     1848                   lf(il,j))*rs(il, j)/(rrv*t(il,j)*t(il,j)*cpd)
    17771849            END IF
    17781850          END IF
    … …  
    17911863
    17921864            IF (cvflag_ice) THEN
    1793               anum = anum - (lv(il,j)+frac(il,j)*lf(il,j))*(rti-rs(il,j)-cwat &
    1794                 *bf2)
     1865              anum = anum - (lv(il,j)+frac(il,j)*lf(il,j))*(rti-rs(il,j)-cwat*bf2)
    17951866              denom = denom + (lv(il,j)+frac(il,j)*lf(il,j))*(rr(il,i)-rti)
    17961867            ELSE
    … …  
    18011872            IF (abs(denom)<0.01) denom = 0.01
    18021873            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)
     1874            altem = sij(il, i, j)*rr(il, i) + (1.-sij(il,i,j))*rti - rs(il, j)
    18051875            altem = altem - (bf2-1.)*cwat
    18061876          END IF
    18071877          IF (sij(il,i,j)>0.0 .AND. sij(il,i,j)<0.95) THEN
    18081878            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
     1879            uent(il, i, j) = sij(il, i, j)*u(il, i) + (1.-sij(il,i,j))*unk(il)
     1880            vent(il, i, j) = sij(il, i, j)*v(il, i) + (1.-sij(il,i,j))*vnk(il)
     1881!!!!      do k=1,ntra
     1882!!!!      traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k)
     1883!!!!     :      +(1.-sij(il,i,j))*tra(il,nk(il),k)
     1884!!!!      end do
    18171885            elij(il, i, j) = altem
    18181886            elij(il, i, j) = max(0.0, elij(il,i,j))
    … …  
    18261894    END DO
    18271895
    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)
     1896!AC!       do k=1,ntra
     1897!AC!        do j=minorig,nl
     1898!AC!         do il=1,ncum
     1899!AC!          if( (i.ge.icb(il)).and.(i.le.inb(il)).and.
     1900!AC!     :       (j.ge.(icb(il)-1)).and.(j.le.inb(il)))then
     1901!AC!            traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k)
     1902!AC!     :            +(1.-sij(il,i,j))*tra(il,nk(il),k)
     1903!AC!          endif
     1904!AC!         enddo
     1905!AC!        enddo
     1906!AC!       enddo
     1907
     1908
     1909! ***   if no air can entrain at level i assume that updraft detrains  ***
     1910! ***   at that level and calculate detrained air flux and properties  ***
     1911
     1912
     1913! @      do 170 i=icb(il),inb(il)
    18481914
    18491915    DO il = 1, ncum
    18501916      IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (nent(il,i)==0)) THEN
    1851         ! @      if(nent(il,i).eq.0)then
     1917! @      if(nent(il,i).eq.0)then
    18521918        ment(il, i, i) = m(il, i)
    18531919        qent(il, i, i) = qnk(il) - ep(il, i)*clw(il, i)
    … …  
    18551921        vent(il, i, i) = vnk(il)
    18561922        elij(il, i, i) = clw(il, i)
    1857         ! MAF      sij(il,i,i)=1.0
     1923! MAF      sij(il,i,i)=1.0
    18581924        sij(il, i, i) = 0.0
    18591925      END IF
    … …  
    18611927  END DO
    18621928
    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
     1929!AC!      do j=1,ntra
     1930!AC!       do i=minorig+1,nl
     1931!AC!        do il=1,ncum
     1932!AC!         if (i.ge.icb(il) .and. i.le.inb(il) .and. nent(il,i).eq.0) then
     1933!AC!          traent(il,i,i,j)=tra(il,nk(il),j)
     1934!AC!         endif
     1935!AC!        enddo
     1936!AC!       enddo
     1937!AC!      enddo
    18731938
    18741939  DO j = minorig, nl
    18751940    DO i = minorig, nl
    18761941      DO il = 1, ncum
    1877         IF ((j>=(icb(il)-1)) .AND. (j<=inb(il)) .AND. (i>=icb(il)) .AND. (i<= &
    1878             inb(il))) THEN
     1942        IF ((j>=(icb(il)-1)) .AND. (j<=inb(il)) .AND. (i>=icb(il)) .AND. (i<=inb(il))) THEN
    18791943          sigij(il, i, j) = sij(il, i, j)
    18801944        END IF
    … …  
    18821946    END DO
    18831947  END DO
    1884   ! @      enddo
    1885 
    1886   ! @170   continue
    1887 
    1888   ! =====================================================================
    1889   ! ---  NORMALIZE ENTRAINED AIR MASS FLUXES
    1890   ! ---  TO REPRESENT EQUAL PROBABILITIES OF MIXING
    1891   ! =====================================================================
     1948! @      enddo
     1949
     1950! @170   continue
     1951
     1952! =====================================================================
     1953! ---  NORMALIZE ENTRAINED AIR MASS FLUXES
     1954! ---  TO REPRESENT EQUAL PROBABILITIES OF MIXING
     1955! =====================================================================
    18921956
    18931957  CALL zilch(asum, nloc*nd)
    … …  
    19151979        IF (cvflag_ice) THEN
    19161980
    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)
     1981          anum = h(il, i) - hp(il, i) - (lv(il,i)+frac(il,i)*lf(il,i))* &
     1982                       (qp-rs(il,i)) + (cpv-cpd)*t(il, i)*(qp-rr(il,i))
     1983          denom = h(il, i) - hp(il, i) + (lv(il,i)+frac(il,i)*lf(il,i))* &
     1984                       (rr(il,i)-qp) + (cpd-cpv)*t(il, i)*(rr(il,i)-qp)
    19211985        ELSE
    19221986
    19231987          anum = h(il, i) - hp(il, i) - lv(il, i)*(qp-rs(il,i)) + &
    1924             (cpv-cpd)*t(il, i)*(qp-rr(il,i))
     1988                       (cpv-cpd)*t(il, i)*(qp-rr(il,i))
    19251989          denom = h(il, i) - hp(il, i) + lv(il, i)*(rr(il,i)-qp) + &
    1926             (cpd-cpv)*t(il, i)*(rr(il,i)-qp)
     1990                       (cpd-cpv)*t(il, i)*(rr(il,i)-qp)
    19271991        END IF
    19281992
    … …  
    19402004      num2 = 0
    19412005      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
     2006        IF (i>=icb(il) .AND. i<=inb(il) .AND. &
     2007            j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
     2008            lwork(il)) num2 = num2 + 1
    19442009      END DO
    19452010      IF (num2<=0) GO TO 175
    19462011
    19472012      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
     2013        IF (i>=icb(il) .AND. i<=inb(il) .AND. &
     2014            j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
     2015            lwork(il)) THEN
    19502016
    19512017          IF (sij(il,i,j)>1.0E-16 .AND. sij(il,i,j)<0.95) THEN
    … …  
    19882054    DO j = minorig, nl
    19892055      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
     2056        IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
     2057            j>=(icb(il)-1) .AND. j<=inb(il)) THEN
    19922058          ment(il, i, j) = ment(il, i, j)*asij(il)
    19932059        END IF
    … …  
    19972063    DO j = minorig, nl
    19982064      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
     2065        IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
     2066            j>=(icb(il)-1) .AND. j<=inb(il)) THEN
    20012067          asum(il, i) = asum(il, i) + ment(il, i, j)
    20022068          ment(il, i, j) = ment(il, i, j)*sig(il, j)
    … …  
    20152081    DO j = minorig, nl
    20162082      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
     2083        IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
     2084            j>=(icb(il)-1) .AND. j<=inb(il)) THEN
    20192085          ment(il, i, j) = ment(il, i, j)*asum(il, i)*bsum(il, i)
    20202086        END IF
    … …  
    20242090    DO j = minorig, nl
    20252091      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
     2092        IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
     2093            j>=(icb(il)-1) .AND. j<=inb(il)) THEN
    20282094          csum(il, i) = csum(il, i) + ment(il, i, j)
    20292095        END IF
    … …  
    20402106        vent(il, i, i) = vnk(il)
    20412107        elij(il, i, i) = clw(il, i)
    2042         ! MAF        sij(il,i,i)=1.0
     2108! MAF        sij(il,i,i)=1.0
    20432109        sij(il, i, i) = 0.0
    20442110      END IF
    20452111    END DO ! il
    20462112
    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
     2113!AC!      do j=1,ntra
     2114!AC!       do il=1,ncum
     2115!AC!        if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il)
     2116!AC!     :     .and. csum(il,i).lt.m(il,i) ) then
     2117!AC!         traent(il,i,i,j)=tra(il,nk(il),j)
     2118!AC!        endif
     2119!AC!       enddo
     2120!AC!      enddo
    20552121789 END DO
    20562122
    2057   ! MAF: renormalisation de MENT
     2123! MAF: renormalisation de MENT
    20582124  CALL zilch(zm, nloc*na)
    20592125  DO jm = 1, nd
    … …  
    20872153END SUBROUTINE cv3_mixing
    20882154
    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
     2155SUBROUTINE cv3_unsat(nloc, ncum, nd, na, ntra, icb, inb, iflag, &
     2156                     t, rr, rs, gz, u, v, tra, p, ph, &
     2157                     th, tv, lv, lf, cpn, ep, sigp, clw, &
     2158                     m, ment, elij, delt, plcl, coef_clos, &
     2159                     mp, rp, up, vp, trap, wt, water, evap, fondue, ice, &
     2160                     faci, b, sigd, &
     2161                     wdtrainA, wdtrainM)                                      ! RomP
    20932162  IMPLICIT NONE
    20942163
    … …  
    20982167  include "cvflag.h"
    20992168
    2100   ! inputs:
     2169!inputs:
    21012170  INTEGER ncum, nd, na, ntra, nloc
    21022171  INTEGER icb(nloc), inb(nloc)
    … …  
    21122181  REAL coef_clos(nloc)
    21132182
    2114   ! input/output
     2183!input/output
    21152184  INTEGER iflag(nloc)
    21162185
    2117   ! outputs:
     2186!outputs:
    21182187  REAL mp(nloc, na), rp(nloc, na), up(nloc, na), vp(nloc, na)
    21192188  REAL water(nloc, na), evap(nloc, na), wt(nloc, na)
    … …  
    21212190  REAL trap(nloc, na, ntra)
    21222191  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
     2192! 25/08/10 - RomP---- ajout des masses precipitantes ejectees
     2193! de l ascendance adiabatique et des flux melanges Pa et Pm.
     2194! Distinction des wdtrain
     2195! Pa = wdtrainA     Pm = wdtrainM
     2196  REAL wdtrainA(nloc, na), wdtrainM(nloc, na)
     2197
     2198!local variables
    21302199  INTEGER i, j, k, il, num1, ndp1
    21312200  REAL tinv, delti, coef
    … …  
    21432212
    21442213
    2145   ! ------------------------------------------------------
     2214! ------------------------------------------------------
    21462215
    21472216  delti = 1./delt
    … …  
    21702239    END DO
    21712240  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 >>>
     2241!AC!        do k=1,ntra
     2242!AC!         do i=1,nd
     2243!AC!          do il=1,ncum
     2244!AC!           trap(il,i,k)=tra(il,i,k)
     2245!AC!          enddo
     2246!AC!         enddo
     2247!AC!        enddo
     2248!! RomP >>>
    21802249  DO i = 1, nd
    21812250    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                      ***
     2251      wdtrainA(il, i) = 0.0
     2252      wdtrainM(il, i) = 0.0
     2253    END DO
     2254  END DO
     2255!! RomP <<<
     2256
     2257! ***  check whether ep(inb)=0, if so, skip precipitating    ***
     2258! ***             downdraft calculation                      ***
    21902259
    21912260
    21922261  DO il = 1, ncum
    2193     ! !          lwork(il)=.TRUE.
    2194     ! !          if(ep(il,inb(il)).lt.0.0001)lwork(il)=.FALSE.
     2262!!          lwork(il)=.TRUE.
     2263!!          if(ep(il,inb(il)).lt.0.0001)lwork(il)=.FALSE.
    21952264    lwork(il) = ep(il, inb(il)) >= 0.0001
    21962265  END DO
    21972266
    2198   ! ***  Set the fractionnal area sigd of precipitating downdraughts
     2267! ***  Set the fractionnal area sigd of precipitating downdraughts
    21992268  DO il = 1, ncum
    22002269    sigd(il) = sigdz*coef_clos(il)
    … …  
    22022271
    22032272
    2204   ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    2205 
    2206   ! ***                    begin downdraft loop                    ***
    2207 
    2208   ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     2273! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     2274!
     2275! ***                    begin downdraft loop                    ***
     2276!
     2277! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    22092278
    22102279  DO i = nl + 1, 1, -1
    … …  
    22192288
    22202289
    2221     ! ***  integrate liquid water equation to find condensed water   ***
    2222     ! ***                and condensed water flux                    ***
    2223 
    2224 
    2225     ! ***              calculate detrained precipitation             ***
     2290! ***  integrate liquid water equation to find condensed water   ***
     2291! ***                and condensed water flux                    ***
     2292!
     2293!
     2294! ***              calculate detrained precipitation             ***
    22262295
    22272296    DO il = 1, ncum
    … …  
    22292298        IF (cvflag_grav) THEN
    22302299          wdtrain(il) = grav*ep(il, i)*m(il, i)*clw(il, i)
    2231           wdtraina(il, i) = wdtrain(il)/grav !   Pa   RomP
     2300          wdtrainA(il, i) = wdtrain(il)/grav                        !   Pa   RomP
    22322301        ELSE
    22332302          wdtrain(il) = 10.0*ep(il, i)*m(il, i)*clw(il, i)
    2234           wdtraina(il, i) = wdtrain(il)/10. !   Pa   RomP
     2303          wdtrainA(il, i) = wdtrain(il)/10.                        !   Pa   RomP
    22352304        END IF
    22362305      END IF
    … …  
    22452314            IF (cvflag_grav) THEN
    22462315              wdtrain(il) = wdtrain(il) + grav*awat*ment(il, j, i)
    2247               wdtrainm(il, i) = wdtrain(il)/grav - wdtraina(il, i) !   Pm  RomP
     2316              wdtrainM(il, i) = wdtrain(il)/grav - wdtrainA(il, i) !   Pm  RomP
    22482317            ELSE
    22492318              wdtrain(il) = wdtrain(il) + 10.0*awat*ment(il, j, i)
    2250               wdtrainm(il, i) = wdtrain(il)/10. - wdtraina(il, i) !   Pm  RomP
     2319              wdtrainM(il, i) = wdtrain(il)/10. - wdtrainA(il, i)  !   Pm  RomP
    22512320            END IF
    22522321          END IF
    … …  
    22562325
    22572326
    2258     ! ***    find rain water and evaporation using provisional   ***
    2259     ! ***              estimates of rp(i)and rp(i-1)             ***
     2327! ***    find rain water and evaporation using provisional   ***
     2328! ***              estimates of rp(i)and rp(i-1)             ***
    22602329
    22612330
    … …  
    22832352          END IF
    22842353
    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)
     2354          rp(il, i) = rp(il, i+1) + &
     2355                      (cpd*(t(il,i+1)-t(il,i))+gz(il,i+1)-gz(il,i))/lv(il, i)
    22872356          rp(il, i) = 0.5*(rp(il,i)+rr(il,i))
    22882357        END IF
    … …  
    22962365          afac = p(il, 1)*(rs(il,1)-rp(il,1))/(1.0E4+2000.0*p(il,1)*rs(il,1))
    22972366          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))
     2367            afac1 = p(il, i)*(rs(il,1)-rp(il,1))/(1.0E4+2000.0*p(il,1)*rs(il,1))
    23002368          END IF
    23012369        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)
     2370          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)
    23042371          rp(il, i-1) = 0.5*(rp(il,i-1)+rr(il,i-1))
    23052372          rp(il, i-1) = amin1(rp(il,i-1), rs(il,i-1))
    23062373          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))
     2374          afac1 = p(il, i)*(rs(il,i)-rp(il,i))/(1.0E4+2000.0*p(il,i)*rs(il,i))
     2375          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))
    23112376          afac = 0.5*(afac1+afac2)
    23122377        END IF
    … …  
    23152380        bfac = 1./(sigd(il)*wt(il,i))
    23162381
    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.
     2382!JYG1
     2383! cc        sigt=1.0
     2384! cc        if(i.ge.icb)sigt=sigp(i)
     2385! prise en compte de la variation progressive de sigt dans
     2386! les couches icb et icb-1:
     2387! pour plcl<ph(i+1), pr1=0 & pr2=1
     2388! pour plcl>ph(i),   pr1=1 & pr2=0
     2389! pour ph(i+1)<plcl<ph(i), pr1 est la proportion a cheval
     2390! sur le nuage, et pr2 est la proportion sous la base du
     2391! nuage.
    23272392        pr1 = (plcl(il)-ph(il,i+1))/(ph(il,i)-ph(il,i+1))
    23282393        pr1 = max(0., min(1.,pr1))
    … …  
    23302395        pr2 = max(0., min(1.,pr2))
    23312396        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.
     2397!JYG2
     2398
     2399!JYG----
     2400!    b6 = bfac*100.*sigd(il)*(ph(il,i)-ph(il,i+1))*sigt*afac
     2401!    c6 = water(il,i+1) + wdtrain(il)*bfac
     2402!    c6 = prec(il,i+1) + wdtrain(il)*bfac
     2403!    revap=0.5*(-b6+sqrt(b6*b6+4.*c6))
     2404!    evap(il,i)=sigt*afac*revap
     2405!    water(il,i)=revap*revap
     2406!    prec(il,i)=revap*revap
     2407!!        print *,' i,b6,c6,revap,evap(il,i),water(il,i),wdtrain(il) ', &
     2408!!                 i,b6,c6,revap,evap(il,i),water(il,i),wdtrain(il)
     2409!!---end jyg---
     2410
     2411! --------retour à la formulation originale d'Emanuel.
    23482412        IF (cvflag_ice) THEN
    23492413
    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)
     2414!  b6=bfac*50.*sigd(il)*(ph(il,i)-ph(il,i+1))*sigt*afac
     2415!   c6=prec(il,i+1)+bfac*wdtrain(il) &
     2416!       -50.*sigd(il)*bfac*(ph(il,i)-ph(il,i+1))*evap(il,i+1)
     2417!  if(c6.gt.0.0)then
     2418!  revap=0.5*(-b6+sqrt(b6*b6+4.*c6))
     2419
     2420!JAM  Attention: evap=sigt*E
     2421!    Modification: evap devient l'évaporation en milieu de couche
     2422!    car nécessaire dans cv3_yield
     2423!    Du coup, il faut modifier pas mal d'équations...
     2424!    et l'expression de afac qui devient afac1
     2425!    revap=sqrt((prec(i+1)+prec(i))/2)
    23622426
    23632427          b6 = bfac*50.*sigd(il)*(ph(il,i)-ph(il,i+1))*sigt*afac1
    23642428          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
     2429! print *,'bfac,sigd(il),sigt,afac1 ',bfac,sigd(il),sigt,afac1
     2430! print *,'prec(il,i+1),wdtrain(il) ',prec(il,i+1),wdtrain(il)
     2431! print *,'b6,c6,b6*b6+4.*c6 ',b6,c6,b6*b6+4.*c6
    23682432          IF (c6>b6*b6+1.E-20) THEN
    23692433            revap = 2.*c6/(b6+sqrt(b6*b6+4.*c6))
    … …  
    23722436          END IF
    23732437          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)
     2438! print*,prec(il,i),'neige'
     2439
     2440!JYG    Dans sa formulation originale, Emanuel calcule l'evaporation par:
     2441! c             evap(il,i)=sigt*afac*revap
     2442! ce qui n'est pas correct. Dans cv_routines, la formulation a été modifiee.
     2443! Ici,l'evaporation evap est simplement calculee par l'equation de
     2444! conservation.
     2445! prec(il,i)=revap*revap
     2446! else
     2447!JYG----   Correction : si c6 <= 0, water(il,i)=0.
     2448! prec(il,i)=0.
     2449! endif
     2450
     2451!JYG---   Dans tous les cas, evaporation = [tt ce qui entre dans la couche i]
     2452! moins [tt ce qui sort de la couche i]
     2453! print *, 'evap avec ice'
     2454          evap(il, i) = (wdtrain(il)+sigd(il)*wt(il,i)*(prec(il,i+1)-prec(il,i))) / &
     2455                        (sigd(il)*(ph(il,i)-ph(il,i+1))*100.)
     2456
     2457          d6 = bfac*wdtrain(il) - 100.*sigd(il)*bfac*(ph(il,i)-ph(il,i+1))*evap(il, i)
    23982458          e6 = bfac*wdtrain(il)
    23992459          f6 = -100.*sigd(il)*bfac*(ph(il,i)-ph(il,i+1))*evap(il, i)
    … …  
    24152475          END IF
    24162476
    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
     2477!          water(il,i)=water(il,i+1)+(1.-fraci(il,i))*e6+(1.-faci(il,i))*f6
     2478!          water(il,i)=max(water(il,i),0.)
     2479!          ice(il,i)=ice(il,i+1)+fraci(il,i)*e6+faci(il,i)*f6
     2480!          ice(il,i)=max(ice(il,i),0.)
     2481!          fondue(il,i)=ice(il,i)*thaw
     2482!          water(il,i)=water(il,i)+fondue(il,i)
     2483!          ice(il,i)=ice(il,i)-fondue(il,i)
     2484           
     2485!          if((water(il,i)+ice(il,i)).lt.1.e-30)then
     2486!            faci(il,i)=0.
     2487!          else
     2488!            faci(il,i)=ice(il,i)/(water(il,i)+ice(il,i))
     2489!          endif
    24302490
    24312491        ELSE
    24322492          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)
     2493          c6 = water(il, i+1) + bfac*wdtrain(il) - &
     2494               50.*sigd(il)*bfac*(ph(il,i)-ph(il,i+1))*evap(il, i+1)
    24352495          IF (c6>0.0) THEN
    24362496            revap = 0.5*(-b6+sqrt(b6*b6+4.*c6))
    … …  
    24392499            water(il, i) = 0.
    24402500          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.)
     2501! print *, 'evap sans ice'
     2502          evap(il, i) = (wdtrain(il)+sigd(il)*wt(il,i)*(water(il,i+1)-water(il,i)))/ &
     2503                        (sigd(il)*(ph(il,i)-ph(il,i+1))*100.)
    24442504
    24452505        END IF
    24462506      END IF !(i.le.inb(il) .and. lwork(il))
    24472507    END DO
    2448     ! ----------------------------------------------------------------
    2449 
    2450     ! cc
    2451     ! ***  calculate precipitating downdraft mass flux under     ***
    2452     ! ***              hydrostatic approximation                 ***
     2508! ----------------------------------------------------------------
     2509
     2510! cc
     2511! ***  calculate precipitating downdraft mass flux under     ***
     2512! ***              hydrostatic approximation                 ***
    24532513
    24542514    DO il = 1, ncum
    … …  
    24592519        IF (cvflag_ice) THEN
    24602520          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)))
     2521            mp(il, i) = 100.*ginv*(lvcp(il,i)*sigd(il)*tevap(il)* &
     2522                                               (p(il,i-1)-p(il,i))/delth + &
     2523                                   lfcp(il,i)*sigd(il)*faci(il,i)*tevap(il)* &
     2524                                               (p(il,i-1)-p(il,i))/delth + &
     2525                                   lfcp(il,i)*sigd(il)*wt(il,i)/100.*fondue(il,i)* &
     2526                                               (p(il,i-1)-p(il,i))/delth/(ph(il,i)-ph(il,i+1)))
    24652527          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)))
     2528            mp(il, i) = 10.*(lvcp(il,i)*sigd(il)*tevap(il)* &
     2529                                                (p(il,i-1)-p(il,i))/delth + &
     2530                             lfcp(il,i)*sigd(il)*faci(il,i)*tevap(il)* &
     2531                                                (p(il,i-1)-p(il,i))/delth + &
     2532                             lfcp(il,i)*sigd(il)*wt(il,i)/100.*fondue(il,i)* &
     2533                                                (p(il,i-1)-p(il,i))/delth/(ph(il,i)-ph(il,i+1)))
    24702534
    24712535          END IF
    … …  
    24732537          IF (cvflag_grav) THEN
    24742538            mp(il, i) = 100.*ginv*lvcp(il, i)*sigd(il)*tevap(il)* &
    2475               (p(il,i-1)-p(il,i))/delth
     2539                                                (p(il,i-1)-p(il,i))/delth
    24762540          ELSE
    24772541            mp(il, i) = 10.*lvcp(il, i)*sigd(il)*tevap(il)* &
    2478               (p(il,i-1)-p(il,i))/delth
     2542                                                (p(il,i-1)-p(il,i))/delth
    24792543          END IF
    24802544
    … …  
    24832547      END IF !(i.le.inb(il) .and. lwork(il) .and. i.ne.1)
    24842548    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 ***
     2549! ----------------------------------------------------------------
     2550
     2551! ***           if hydrostatic assumption fails,             ***
     2552! ***   solve cubic difference equation for downdraft theta  ***
     2553! ***  and mass flux from two simultaneous differential eqns ***
    24902554
    24912555    DO il = 1, ncum
    … …  
    24932557
    24942558        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))
     2559                         (th(il,i)-th(il,i-1))/(tv(il,i)*th(il,i))
    24962560        amp2 = abs(mp(il,i+1)*mp(il,i+1)-mp(il,i)*mp(il,i))
    24972561
    24982562        IF (amp2>(0.1*amfac)) THEN
    24992563          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))
     2564          tf = b(il, i) - 5.0*(th(il,i)-th(il,i-1))*t(il, i) / &
     2565                              (lvcp(il,i)*sigd(il)*th(il,i))
    25022566          af = xf*tf + mp(il, i+1)*mp(il, i+1)*tinv
    25032567
    25042568          IF (cvflag_ice) THEN
    25052569            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)))
     2570                 50.*(p(il,i-1)-p(il,i))*xf*(tevap(il)*(1.+(lf(il,i)/lv(il,i))*faci(il,i)) + &
     2571                (lf(il,i)/lv(il,i))*wt(il,i)/100.*fondue(il,i)/(ph(il,i)-ph(il,i+1)))
    25092572          ELSE
    25102573
    25112574            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)
     2575                                           50.*(p(il,i-1)-p(il,i))*xf*tevap(il)
    25132576          END IF
    25142577
    … …  
    25222585            IF ((0.5*bf-sru)<0.0) fac = -1.0
    25232586            mp(il, i) = mp(il, i+1)*tinv + (0.5*bf+sru)**tinv + &
    2524               fac*(abs(0.5*bf-sru))**tinv
     2587                                           fac*(abs(0.5*bf-sru))**tinv
    25252588          ELSE
    25262589            d = atan(2.*sqrt(-ur)/(bf+1.0E-28))
    … …  
    25322595          IF (cvflag_ice) THEN
    25332596            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))
     2597!JYG : il y a vraisemblablement une erreur dans la ligne 2 suivante:
     2598! il faut diviser par (mp(il,i)*sigd(il)*grav) et non par (mp(il,i)+sigd(il)*0.1).
     2599! Et il faut bien revoir les facteurs 100.
     2600              b(il, i-1) = b(il, i) + 100.0*(p(il,i-1)-p(il,i))* &
     2601                           (tevap(il)*(1.+(lf(il,i)/lv(il,i))*faci(il,i)) + &
     2602                           (lf(il,i)/lv(il,i))*wt(il,i)/100.*fondue(il,i) / &
     2603                           (ph(il,i)-ph(il,i+1))) / &
     2604                           (mp(il,i)+sigd(il)*0.1) - &
     2605                           10.0*(th(il,i)-th(il,i-1))*t(il, i) / &
     2606                           (lvcp(il,i)*sigd(il)*th(il,i))
    25442607            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))
     2608              b(il, i-1) = b(il, i) + 100.0*(p(il,i-1)-p(il,i))*&
     2609                           (tevap(il)*(1.+(lf(il,i)/lv(il,i))*faci(il,i)) + &
     2610                           (lf(il,i)/lv(il,i))*wt(il,i)/100.*fondue(il,i) / &
     2611                           (ph(il,i)-ph(il,i+1))) / &
     2612                           (mp(il,i)+sigd(il)*0.1) - &
     2613                           10.0*(th(il,i)-th(il,i-1))*t(il, i) / &
     2614                           (lvcp(il,i)*sigd(il)*th(il,i))
    25502615            END IF
    25512616          ELSE
    25522617            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))
     2618              b(il, i-1) = b(il, i) + 100.0*(p(il,i-1)-p(il,i))*tevap(il) / &
     2619                           (mp(il,i)+sigd(il)*0.1) - &
     2620                           10.0*(th(il,i)-th(il,i-1))*t(il, i) / &
     2621                           (lvcp(il,i)*sigd(il)*th(il,i))
    25562622            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))
     2623              b(il, i-1) = b(il, i) + 100.0*(p(il,i-1)-p(il,i))*tevap(il) / &
     2624                           (mp(il,i)+sigd(il)*0.1) - &
     2625                           10.0*(th(il,i)-th(il,i-1))*t(il, i) / &
     2626                           (lvcp(il,i)*sigd(il)*th(il,i))
    25602627            END IF
    25612628          END IF
    … …  
    25642631        END IF !(amp2.gt.(0.1*amfac))
    25652632
    2566         ! ***         limit magnitude of mp(i) to meet cfl condition      ***
     2633! ***         limit magnitude of mp(i) to meet cfl condition      ***
    25672634
    25682635        ampmax = 2.0*(ph(il,i)-ph(il,i+1))*delti
    … …  
    25712638        mp(il, i) = min(mp(il,i), ampmax)
    25722639
    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
     2640! ***      force mp to decrease linearly to zero                 ***
     2641! ***       between cloud base and the surface                   ***
     2642
     2643
     2644! c      if(p(il,i).gt.p(il,icb(il)))then
     2645! c       mp(il,i)=mp(il,icb(il))*(p(il,1)-p(il,i))/(p(il,1)-p(il,icb(il)))
     2646! c      endif
    25812647        IF (ph(il,i)>0.9*plcl(il)) THEN
    25822648          mp(il, i) = mp(il, i)*(ph(il,1)-ph(il,i))/(ph(il,1)-0.9*plcl(il))
    … …  
    25852651      END IF ! (i.le.inb(il) .and. lwork(il) .and. i.ne.1)
    25862652    END DO
    2587     ! ----------------------------------------------------------------
    2588 
    2589     ! ***       find mixing ratio of precipitating downdraft     ***
     2653! ----------------------------------------------------------------
     2654
     2655! ***       find mixing ratio of precipitating downdraft     ***
    25902656
    25912657    DO il = 1, ncum
    … …  
    26032669
    26042670          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))
     2671            rp(il, i) = rp(il, i+1)*mp(il, i+1) + rr(il, i)*(mp(il,i)-mp(il,i+1)) + &
     2672              100.*ginv*0.5*sigd(il)*(ph(il,i)-ph(il,i+1))*(evap(il,i+1)+evap(il,i))
    26082673          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))
     2674            rp(il, i) = rp(il, i+1)*mp(il, i+1) + rr(il, i)*(mp(il,i)-mp(il,i+1)) + &
     2675              5.*sigd(il)*(ph(il,i)-ph(il,i+1))*(evap(il,i+1)+evap(il,i))
    26122676          END IF
    26132677          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             )
     2678          up(il, i) = up(il, i+1)*mp(il, i+1) + u(il, i)*(mp(il,i)-mp(il,i+1))
    26162679          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             )
     2680          vp(il, i) = vp(il, i+1)*mp(il, i+1) + v(il, i)*(mp(il,i)-mp(il,i+1))
    26192681          vp(il, i) = vp(il, i)/mp(il, i)
    26202682
    … …  
    26232685          IF (mp(il,i+1)>1.0E-16) THEN
    26242686            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)
     2687              rp(il, i) = rp(il,i+1) + 100.*ginv*0.5*sigd(il)*(ph(il,i)-ph(il,i+1)) * &
     2688                                       (evap(il,i+1)+evap(il,i))/mp(il,i+1)
    26272689            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)
     2690              rp(il, i) = rp(il,i+1) + 5.*sigd(il)*(ph(il,i)-ph(il,i+1)) * &
     2691                                       (evap(il,i+1)+evap(il,i))/mp(il, i+1)
    26302692            END IF
    26312693            up(il, i) = up(il, i+1)
    … …  
    26392701      END IF ! (i.lt.inb(il) .and. lwork(il))
    26402702    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
     2703! ----------------------------------------------------------------
     2704
     2705! ***       find tracer concentrations in precipitating downdraft     ***
     2706
     2707!AC!      do j=1,ntra
     2708!AC!       do il = 1,ncum
     2709!AC!       if (i.lt.inb(il) .and. lwork(il)) then
     2710!AC!c
     2711!AC!         if(mplus(il))then
     2712!AC!          trap(il,i,j)=trap(il,i+1,j)*mp(il,i+1)
     2713!AC!     :              +trap(il,i,j)*(mp(il,i)-mp(il,i+1))
     2714!AC!          trap(il,i,j)=trap(il,i,j)/mp(il,i)
     2715!AC!         else ! if (mplus(il))
     2716!AC!          if(mp(il,i+1).gt.1.0e-16)then
     2717!AC!           trap(il,i,j)=trap(il,i+1,j)
     2718!AC!          endif
     2719!AC!         endif ! (mplus(il)) else if (.not.mplus(il))
     2720!AC!c
     2721!AC!        endif ! (i.lt.inb(il) .and. lwork(il))
     2722!AC!       enddo
     2723!AC!      end do
    26622724
    26632725400 END DO
    2664   ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    2665 
    2666   ! ***                    end of downdraft loop                    ***
    2667 
    2668   ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     2726! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     2727
     2728! ***                    end of downdraft loop                    ***
     2729
     2730! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    26692731
    26702732
    … …  
    26722734END SUBROUTINE cv3_unsat
    26732735
    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)
     2736SUBROUTINE cv3_yield(nloc, ncum, nd, na, ntra, ok_conserv_q, &
     2737                     icb, inb, delt, &
     2738                     t, rr, t_wake, rr_wake, s_wake, u, v, tra, &
     2739                     gz, p, ph, h, hp, lv, lf, cpn, th, th_wake, &
     2740                     ep, clw, m, tp, mp, rp, up, vp, trap, &
     2741                     wt, water, ice, evap, fondue, faci, b, sigd, &
     2742                     ment, qent, hent, iflag_mix, uent, vent, &
     2743                     nent, elij, traent, sig, &
     2744                     tv, tvp, wghti, &
     2745                     iflag, precip, Vprecip, ft, fr, fu, fv, ftra, &
     2746                     cbmf, upwd, dnwd, dnwd0, ma, mip, &
     2747                     tls, tps, qcondc, wd, &
     2748                     ftd, fqd)
    26802749
    26812750  IMPLICIT NONE
    … …  
    26862755  include "conema3.h"
    26872756
    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:
     2757!inputs:
     2758      INTEGER iflag_mix
     2759      INTEGER ncum, nd, na, ntra, nloc
     2760      LOGICAL ok_conserv_q
     2761      INTEGER icb(nloc), inb(nloc)
     2762      REAL delt
     2763      REAL t(nloc, nd), rr(nloc, nd), u(nloc, nd), v(nloc, nd)
     2764      REAL t_wake(nloc, nd), rr_wake(nloc, nd)
     2765      REAL s_wake(nloc)
     2766      REAL tra(nloc, nd, ntra), sig(nloc, nd)
     2767      REAL gz(nloc, na), ph(nloc, nd+1), h(nloc, na), hp(nloc, na)
     2768      REAL th(nloc, na), p(nloc, nd), tp(nloc, na)
     2769      REAL lv(nloc, na), cpn(nloc, na), ep(nloc, na), clw(nloc, na)
     2770      REAL lf(nloc, na)
     2771      REAL m(nloc, na), mp(nloc, na), rp(nloc, na), up(nloc, na)
     2772      REAL vp(nloc, na), wt(nloc, nd), trap(nloc, nd, ntra)
     2773      REAL water(nloc, na), evap(nloc, na), b(nloc, na), sigd(nloc)
     2774      REAL fondue(nloc, na), faci(nloc, na), ice(nloc, na)
     2775      REAL ment(nloc, na, na), qent(nloc, na, na), uent(nloc, na, na)
     2776      REAL hent(nloc, na, na)
     2777!IM bug   real vent(nloc,na,na), nent(nloc,na), elij(nloc,na,na)
     2778      REAL vent(nloc, na, na), elij(nloc, na, na)
     2779      INTEGER nent(nloc, nd)
     2780      REAL traent(nloc, na, na, ntra)
     2781      REAL tv(nloc, nd), tvp(nloc, nd), wghti(nloc, nd)
     2782!
     2783!input/output:
     2784      INTEGER iflag(nloc)
     2785!
     2786!outputs:
     2787      REAL precip(nloc)
     2788      REAL ft(nloc, nd), fr(nloc, nd), fu(nloc, nd), fv(nloc, nd)
     2789      REAL ftd(nloc, nd), fqd(nloc, nd)
     2790      REAL ftra(nloc, nd, ntra)
     2791      REAL upwd(nloc, nd), dnwd(nloc, nd), ma(nloc, nd)
     2792      REAL dnwd0(nloc, nd), mip(nloc, nd)
     2793      REAL Vprecip(nloc, nd+1)
     2794      REAL tls(nloc, nd), tps(nloc, nd)
     2795      REAL qcondc(nloc, nd) ! cld
     2796      REAL wd(nloc) ! gust
     2797      REAL cbmf(nloc)
     2798!
     2799!local variables:
     2800      INTEGER i, k, il, n, j, num1
     2801      REAL rat, delti
     2802      REAL ax, bx, cx, dx, ex
     2803      REAL cpinv, rdcp, dpinv
     2804      REAL awat(nloc)
     2805      REAL lvcp(nloc, na), lfcp(nloc, na), mke(nloc, na)
     2806      REAL am(nloc), work(nloc), ad(nloc), amp1(nloc)
     2807!!      real up1(nloc), dn1(nloc)
     2808      REAL up1(nloc, nd, nd), dn1(nloc, nd, nd)
     2809      REAL asum(nloc), bsum(nloc), csum(nloc), dsum(nloc)
     2810      REAL esum(nloc), fsum(nloc), gsum(nloc), hsum(nloc)
     2811      REAL th_wake(nloc, nd)
     2812      REAL alpha_qpos(nloc), alpha_qpos1(nloc)
     2813      REAL qcond(nloc, nd), nqcond(nloc, nd), wa(nloc, nd) ! cld
     2814      REAL siga(nloc, nd), sax(nloc, nd), mac(nloc, nd) ! cld
     2815
     2816      REAL sumdq !jyg
     2817!
     2818! -------------------------------------------------------------
     2819
     2820! initialization:
    27512821
    27522822  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  )
     2823! print*,'cv3_yield initialisation delt', delt
     2824!
    27582825  DO il = 1, ncum
    27592826    precip(il) = 0.0
    2760     vprecip(il, nd+1) = 0.0
     2827    Vprecip(il, nd+1) = 0.0
    27612828    wd(il) = 0.0 ! gust
    27622829  END DO
    … …  
    27642831  DO i = 1, nd
    27652832    DO il = 1, ncum
    2766       vprecip(il, i) = 0.0
     2833      Vprecip(il, i) = 0.0
    27672834      ft(il, i) = 0.0
    27682835      fr(il, i) = 0.0
    … …  
    27802847    END DO
    27812848  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'
     2849! print*,'cv3_yield initialisation 2'
     2850!AC!      do j=1,ntra
     2851!AC!       do i=1,nd
     2852!AC!        do il=1,ncum
     2853!AC!          ftra(il,i,j)=0.0
     2854!AC!        enddo
     2855!AC!       enddo
     2856!AC!      enddo
     2857! print*,'cv3_yield initialisation 3'
    27912858  DO i = 1, nl
    27922859    DO il = 1, ncum
    … …  
    27982865
    27992866
    2800   ! ***  calculate surface precipitation in mm/day     ***
     2867! ***  calculate surface precipitation in mm/day     ***
    28012868
    28022869  DO il = 1, ncum
    28032870    IF (ep(il,inb(il))>=0.0001 .AND. iflag(il)<=1) THEN
    28042871      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
     2872        precip(il) = wt(il, 1)*sigd(il)*(water(il,1)+ice(il,1)) &
     2873                              *86400.*1000./(rowl*grav)
    28112874      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
     2875        precip(il) = wt(il, 1)*sigd(il)*water(il, 1) &
     2876                              *86400.*1000./(rowl*grav)
    28182877      END IF
    28192878    END IF
    28202879  END DO
    2821   ! print*,'cv3_yield apres calcul precip'
    2822 
    2823 
    2824   ! ===  calculate vertical profile of  precipitation in kg/m2/s  ===
     2880! print*,'cv3_yield apres calcul precip'
     2881
     2882
     2883! ===  calculate vertical profile of  precipitation in kg/m2/s  ===
    28252884
    28262885  DO i = 1, nl
    … …  
    28282887      IF (ep(il,inb(il))>=0.0001 .AND. i<=inb(il) .AND. iflag(il)<=1) THEN
    28292888        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
     2889          Vprecip(il, i) = wt(il, i)*sigd(il)*(water(il,i)+ice(il,i))/grav
    28352890        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
     2891          Vprecip(il, i) = wt(il, i)*sigd(il)*water(il, i)/grav
    28412892        END IF
    28422893      END IF
    … …  
    28452896
    28462897
    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                        ***
     2898! ***  Calculate downdraft velocity scale    ***
     2899! ***  NE PAS UTILISER POUR L'INSTANT ***
     2900
     2901!!      do il=1,ncum
     2902!!        wd(il)=betad*abs(mp(il,icb(il)))*0.01*rrd*t(il,icb(il)) &
     2903!!                                       /(sigd(il)*p(il,icb(il)))
     2904!!      enddo
     2905
     2906
     2907! ***  calculate tendencies of lowest level potential temperature  ***
     2908! ***                      and mixing ratio                        ***
    28582909
    28592910  DO il = 1, ncum
    … …  
    28702921  END DO
    28712922
    2872   ! print*,'cv3_yield avant ft'
    2873   ! AM is the part of cbmf taken from the first level
     2923!    print*,'cv3_yield avant ft'
     2924! am is the part of cbmf taken from the first level
    28742925  DO il = 1, ncum
    28752926    am(il) = cbmf(il)*wghti(il, 1)
    … …  
    28782929  DO il = 1, ncum
    28792930    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
     2931! convect3      if((0.1*dpinv*am).ge.delti)iflag(il)=4
     2932!JYG  Correction pour conserver l'eau
     2933! cc       ft(il,1)=-0.5*lvcp(il,1)*sigd(il)*(evap(il,1)+evap(il,2))          !precip
    28842934      IF (cvflag_ice) THEN
    28852935        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
     2936                     lfcp(il, 1)*sigd(il)*evap(il, 1)*faci(il, 1) - &
     2937                     lfcp(il, 1)*sigd(il)*(fondue(il,1)*wt(il,1)) / &
     2938                       (100.*(ph(il,1)-ph(il,2)))                            !precip
    28892939      ELSE
    28902940        ft(il, 1) = -lvcp(il, 1)*sigd(il)*evap(il, 1)
    28912941      END IF
    28922942
    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)
     2943      ft(il, 1) = ft(il, 1) - 0.009*grav*sigd(il)*mp(il, 2)*t_wake(il, 1)*b(il, 1)*work(il)
     2944
     2945      IF (cvflag_ice) THEN
     2946        ft(il, 1) = ft(il, 1) + 0.01*sigd(il)*wt(il, 1)*(cl-cpd)*water(il, 2) * &
     2947                                     (t_wake(il,2)-t_wake(il,1))*work(il)/cpn(il, 1) + &
     2948                                0.01*sigd(il)*wt(il, 1)*(ci-cpd)*ice(il, 2) * &
     2949                                     (t_wake(il,2)-t_wake(il,1))*work(il)/cpn(il, 1)
    28962950      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)
     2951        ft(il, 1) = ft(il, 1) + 0.01*sigd(il)*wt(il, 1)*(cl-cpd)*water(il, 2) * &
     2952                                     (t_wake(il,2)-t_wake(il,1))*work(il)/cpn(il, 1)
    28992953      END IF
    29002954
    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
     2955      ftd(il, 1) = ft(il, 1)                                                  ! fin precip
     2956
     2957      IF ((0.01*grav*work(il)*am(il))>=delti) iflag(il) = 1 !consist vect
     2958      ft(il, 1) = ft(il, 1) + 0.01*grav*work(il)*am(il) * &
     2959                                   (t(il,2)-t(il,1)+(gz(il,2)-gz(il,1))/cpn(il,1))
    29222960    END IF ! iflag
    29232961  END DO
    … …  
    29272965    IF (iflag_mix>0) THEN
    29282966      DO il = 1, ncum
    2929         ! FH WARNING a modifier :
     2967! FH WARNING a modifier :
    29302968        cpinv = 0.
    2931         ! cpinv=1.0/cpn(il,1)
     2969! cpinv=1.0/cpn(il,1)
    29322970        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
     2971          ft(il, 1) = ft(il, 1) + 0.01*grav*work(il)*ment(il, j, 1) * &
     2972                     (hent(il,j,1)-h(il,1)+t(il,1)*(cpv-cpd)*(rr(il,1)-qent(il,j,1)))*cpinv
    29412973        END IF ! j
    29422974      END DO
    29432975    END IF
    29442976  END DO
    2945     ! fin sature
     2977! fin sature
    29462978
    29472979
    29482980  DO il = 1, ncum
    29492981    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
     2982!JYG1  Correction pour mieux conserver l'eau (conformite avec CONVECT4.3)
     2983      fr(il, 1) = 0.01*grav*mp(il, 2)*(rp(il,2)-rr_wake(il,1))*work(il) + &
     2984                  sigd(il)*evap(il, 1)
     2985!!!                  sigd(il)*0.5*(evap(il,1)+evap(il,2))
     2986
     2987      fqd(il, 1) = fr(il, 1) !precip
     2988
     2989      fr(il, 1) = fr(il, 1) + 0.01*grav*am(il)*(rr(il,2)-rr(il,1))*work(il)        !sature
     2990
     2991      fu(il, 1) = fu(il, 1) + 0.01*grav*work(il)*(mp(il,2)*(up(il,2)-u(il,1)) + &
     2992                                                  am(il)*(u(il,2)-u(il,1)))
     2993      fv(il, 1) = fv(il, 1) + 0.01*grav*work(il)*(mp(il,2)*(vp(il,2)-v(il,1)) + &
     2994                                                  am(il)*(v(il,2)-v(il,1)))
    29762995    END IF ! iflag
    29772996  END DO ! il
    29782997
    29792998
    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
     2999!AC!     do j=1,ntra
     3000!AC!      do il=1,ncum
     3001!AC!       if (iflag(il) .le. 1) then
     3002!AC!       if (cvflag_grav) then
     3003!AC!        ftra(il,1,j)=ftra(il,1,j)+0.01*grav*work(il)
     3004!AC!    :                     *(mp(il,2)*(trap(il,2,j)-tra(il,1,j))
     3005!AC!    :             +am(il)*(tra(il,2,j)-tra(il,1,j)))
     3006!AC!       else
     3007!AC!        ftra(il,1,j)=ftra(il,1,j)+0.1*work(il)
     3008!AC!    :                     *(mp(il,2)*(trap(il,2,j)-tra(il,1,j))
     3009!AC!    :             +am(il)*(tra(il,2,j)-tra(il,1,j)))
     3010!AC!       endif
     3011!AC!       endif  ! iflag
     3012!AC!      enddo
     3013!AC!     enddo
    29953014
    29963015  DO j = 2, nl
    29973016    DO il = 1, ncum
    29983017      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
     3018        fr(il, 1) = fr(il, 1) + 0.01*grav*work(il)*ment(il, j, 1)*(qent(il,j,1)-rr(il,1))
     3019        fu(il, 1) = fu(il, 1) + 0.01*grav*work(il)*ment(il, j, 1)*(uent(il,j,1)-u(il,1))
     3020        fv(il, 1) = fv(il, 1) + 0.01*grav*work(il)*ment(il, j, 1)*(vent(il,j,1)-v(il,1))
    30143021      END IF ! j
    30153022    END DO
    30163023  END DO
    30173024
    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                          ***
     3025!AC!      do k=1,ntra
     3026!AC!       do j=2,nl
     3027!AC!        do il=1,ncum
     3028!AC!         if (j.le.inb(il) .and. iflag(il) .le. 1) then
     3029!AC!
     3030!AC!          if (cvflag_grav) then
     3031!AC!           ftra(il,1,k)=ftra(il,1,k)+0.01*grav*work(il)*ment(il,j,1)
     3032!AC!     :                *(traent(il,j,1,k)-tra(il,1,k))
     3033!AC!          else
     3034!AC!           ftra(il,1,k)=ftra(il,1,k)+0.1*work(il)*ment(il,j,1)
     3035!AC!     :                *(traent(il,j,1,k)-tra(il,1,k))
     3036!AC!          endif
     3037!AC!
     3038!AC!         endif
     3039!AC!        enddo
     3040!AC!       enddo
     3041!AC!      enddo
     3042! print*,'cv3_yield apres ft'
     3043
     3044! ***  calculate tendencies of potential temperature and mixing ratio  ***
     3045! ***               at levels above the lowest level                   ***
     3046
     3047! ***  first find the net saturated updraft and downdraft mass fluxes  ***
     3048! ***                      through each level                          ***
    30423049
    30433050
    … …  
    30603067          END IF
    30613068        ELSE
    3062           ! AMP1 is the part of cbmf taken from layers I and lower
     3069! AMP1 is the part of cbmf taken from layers I and lower
    30633070          IF (k<=i) THEN
    30643071            amp1(il) = amp1(il) + cbmf(il)*wghti(il, k)
    … …  
    30933100        cpinv = 1.0/cpn(il, i)
    30943101
    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))
     3102! convect3      if((0.1*dpinv*amp1).ge.delti)iflag(il)=4
     3103        IF ((0.01*grav*dpinv*amp1(il))>=delti) iflag(il) = 1 ! vecto
     3104
     3105! precip
     3106! cc       ft(il,i)= -0.5*sigd(il)*lvcp(il,i)*(evap(il,i)+evap(il,i+1))
    31053107        IF (cvflag_ice) THEN
    31063108          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)))
     3109                       sigd(il)*lfcp(il, i)*evap(il, i)*faci(il, i) - &
     3110                       sigd(il)*lfcp(il, i)*fondue(il, i)*wt(il, i)/(100.*(p(il,i-1)-p(il,i)))
    31103111        ELSE
    31113112          ft(il, i) = -sigd(il)*lvcp(il, i)*evap(il, i)
    … …  
    31143115        rat = cpn(il, i-1)*cpinv
    31153116
    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
     3117        ft(il, i) = ft(il, i) - 0.009*grav*sigd(il) * &
     3118                     (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
     3119        IF (cvflag_ice) THEN
     3120          ft(il, i) = ft(il, i) + 0.01*sigd(il)*wt(il, i)*(cl-cpd)*water(il, i+1) * &
     3121                                       (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, i+1) * &
     3126                                       (t_wake(il,i+1)-t_wake(il,i))*dpinv* &
     3127            cpinv
     3128        END IF
     3129
     3130        ftd(il, i) = ft(il, i)
     3131! fin precip
     3132
     3133! sature
     3134        ft(il, i) = ft(il, i) + 0.01*grav*dpinv * &
     3135                     (amp1(il)*(t(il,i+1)-t(il,i) + (gz(il,i+1)-gz(il,i))*cpinv) - &
     3136                      ad(il)*(t(il,i)-t(il,i-1)+(gz(il,i)-gz(il,i-1))*cpinv))
     3137
     3138
     3139        IF (iflag_mix==0) THEN
     3140          ft(il, i) = ft(il, i) + 0.01*grav*dpinv*ment(il, i, i)*(hp(il,i)-h(il,i) + &
     3141                                    t(il,i)*(cpv-cpd)*(rr(il,i)-qent(il,i,i)))*cpinv
     3142        END IF
     3143
     3144
     3145
     3146! sb: on ne fait pas encore la correction permettant de mieux
     3147! conserver l'eau:
     3148!JYG: correction permettant de mieux conserver l'eau:
     3149! cc         fr(il,i)=0.5*sigd(il)*(evap(il,i)+evap(il,i+1))
     3150        fr(il, i) = sigd(il)*evap(il, i) + 0.01*grav*(mp(il,i+1)*(rp(il,i+1)-rr_wake(il,i)) - &
     3151                                                      mp(il,i)*(rp(il,i)-rr_wake(il,i-1)))*dpinv
     3152        fqd(il, i) = fr(il, i)                                                                     ! precip
     3153
     3154        fu(il, i) = 0.01*grav*(mp(il,i+1)*(up(il,i+1)-u(il,i)) - &
     3155                               mp(il,i)*(up(il,i)-u(il,i-1)))*dpinv
     3156        fv(il, i) = 0.01*grav*(mp(il,i+1)*(vp(il,i+1)-v(il,i)) - &
     3157                               mp(il,i)*(vp(il,i)-v(il,i-1)))*dpinv
     3158
     3159
     3160        fr(il, i) = fr(il, i) + 0.01*grav*dpinv*(amp1(il)*(rr(il,i+1)-rr(il,i)) - &
     3161                                                 ad(il)*(rr(il,i)-rr(il,i-1)))
     3162        fu(il, i) = fu(il, i) + 0.01*grav*dpinv*(amp1(il)*(u(il,i+1)-u(il,i)) - &
     3163                                                 ad(il)*(u(il,i)-u(il,i-1)))
     3164        fv(il, i) = fv(il, i) + 0.01*grav*dpinv*(amp1(il)*(v(il,i+1)-v(il,i)) - &
     3165                                                 ad(il)*(v(il,i)-v(il,i-1)))
    32143166
    32153167      END IF ! i
    32163168    END DO
    32173169
    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
     3170!AC!      do k=1,ntra
     3171!AC!       do il=1,ncum
     3172!AC!        if (i.le.inb(il) .and. iflag(il) .le. 1) then
     3173!AC!         dpinv=1.0/(ph(il,i)-ph(il,i+1))
     3174!AC!         cpinv=1.0/cpn(il,i)
     3175!AC!         if (cvflag_grav) then
     3176!AC!           ftra(il,i,k)=ftra(il,i,k)+0.01*grav*dpinv
     3177!AC!     :         *(amp1(il)*(tra(il,i+1,k)-tra(il,i,k))
     3178!AC!     :           -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
     3179!AC!         else
     3180!AC!           ftra(il,i,k)=ftra(il,i,k)+0.1*dpinv
     3181!AC!     :         *(amp1(il)*(tra(il,i+1,k)-tra(il,i,k))
     3182!AC!     :           -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
     3183!AC!         endif
     3184!AC!        endif
     3185!AC!       enddo
     3186!AC!      enddo
    32353187
    32363188    DO k = 1, i - 1
    … …  
    32463198            dpinv = 1.0/(ph(il,i)-ph(il,i+1))
    32473199            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
     3200            ft(il, i) = ft(il, i) + 0.01*grav*dpinv*ment(il, k, i) * &
     3201                 (hent(il,k,i)-h(il,i)+t(il,i)*(cpv-cpd)*(rr(il,i)+awat(il)-qent(il,k,i)))*cpinv
     3202!
     3203!
    32603204          END IF ! i
    32613205        END DO
    … …  
    32663210          dpinv = 1.0/(ph(il,i)-ph(il,i+1))
    32673211          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
     3212          fr(il, i) = fr(il, i) + 0.01*grav*dpinv*ment(il, k, i) * &
     3213                                                       (qent(il,k,i)-awat(il)-rr(il,i))
     3214          fu(il, i) = fu(il, i) + 0.01*grav*dpinv*ment(il, k, i)*(uent(il,k,i)-u(il,i))
     3215          fv(il, i) = fv(il, i) + 0.01*grav*dpinv*ment(il, k, i)*(vent(il,k,i)-v(il,i))
     3216
     3217! (saturated updrafts resulting from mixing)                                   ! cld
     3218          qcond(il, i) = qcond(il, i) + (elij(il,k,i)-awat(il))                ! cld
    32863219          nqcond(il, i) = nqcond(il, i) + 1. ! cld
    32873220        END IF ! i
    … …  
    32893222    END DO
    32903223
    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
     3224!AC!      do j=1,ntra
     3225!AC!       do k=1,i-1
     3226!AC!        do il=1,ncum
     3227!AC!         if (i.le.inb(il) .and. iflag(il) .le. 1) then
     3228!AC!          dpinv=1.0/(ph(il,i)-ph(il,i+1))
     3229!AC!          cpinv=1.0/cpn(il,i)
     3230!AC!          if (cvflag_grav) then
     3231!AC!           ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv*ment(il,k,i)
     3232!AC!     :        *(traent(il,k,i,j)-tra(il,i,j))
     3233!AC!          else
     3234!AC!           ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv*ment(il,k,i)
     3235!AC!     :        *(traent(il,k,i,j)-tra(il,i,j))
     3236!AC!          endif
     3237!AC!         endif
     3238!AC!        enddo
     3239!AC!       enddo
     3240!AC!      enddo
    33083241
    33093242    DO k = i, nl + 1
    … …  
    33143247            dpinv = 1.0/(ph(il,i)-ph(il,i+1))
    33153248            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
     3249            ft(il, i) = ft(il, i) + 0.01*grav*dpinv*ment(il, k, i) * &
     3250                  (hent(il,k,i)-h(il,i)+t(il,i)*(cpv-cpd)*(rr(il,i)-qent(il,k,i)))*cpinv
     3251
     3252
    33263253          END IF ! i
    33273254        END DO
    … …  
    33333260          cpinv = 1.0/cpn(il, i)
    33343261
    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
     3262          fr(il, i) = fr(il, i) + 0.01*grav*dpinv*ment(il, k, i)*(qent(il,k,i)-rr(il,i))
     3263          fu(il, i) = fu(il, i) + 0.01*grav*dpinv*ment(il, k, i)*(uent(il,k,i)-u(il,i))
     3264          fv(il, i) = fv(il, i) + 0.01*grav*dpinv*ment(il, k, i)*(vent(il,k,i)-v(il,i))
    33503265        END IF ! i and k
    33513266      END DO
    33523267    END DO
    33533268
    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
     3269!AC!      do j=1,ntra
     3270!AC!       do k=i,nl+1
     3271!AC!        do il=1,ncum
     3272!AC!         if (i.le.inb(il) .and. k.le.inb(il)
     3273!AC!     $                .and. iflag(il) .le. 1) then
     3274!AC!          dpinv=1.0/(ph(il,i)-ph(il,i+1))
     3275!AC!          cpinv=1.0/cpn(il,i)
     3276!AC!          if (cvflag_grav) then
     3277!AC!           ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv*ment(il,k,i)
     3278!AC!     :         *(traent(il,k,i,j)-tra(il,i,j))
     3279!AC!          else
     3280!AC!           ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv*ment(il,k,i)
     3281!AC!     :             *(traent(il,k,i,j)-tra(il,i,j))
     3282!AC!          endif
     3283!AC!         endif ! i and k
     3284!AC!        enddo
     3285!AC!       enddo
     3286!AC!      enddo
     3287
     3288! sb: interface with the cloud parameterization:                               ! cld
    33743289
    33753290    DO k = i + 1, nl
    33763291      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
     3292        IF (k<=inb(il) .AND. i<=inb(il) .AND. iflag(il)<=1) THEN               ! cld
     3293! (saturated downdrafts resulting from mixing)                                 ! cld
     3294          qcond(il, i) = qcond(il, i) + elij(il, k, i)                         ! cld
    33803295          nqcond(il, i) = nqcond(il, i) + 1. ! cld
    33813296        END IF ! cld
    … …  
    33833298    END DO ! cld
    33843299
    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
     3300! (particular case: no detraining level is found)                              ! cld
     3301    DO il = 1, ncum                                                            ! cld
     3302      IF (i<=inb(il) .AND. nent(il,i)==0 .AND. iflag(il)<=1) THEN              ! cld
     3303        qcond(il, i) = qcond(il, i) + (1.-ep(il,i))*clw(il, i)                 ! cld
     3304        nqcond(il, i) = nqcond(il, i) + 1.                                     ! cld
     3305      END IF                                                                   ! cld
     3306    END DO                                                                     ! cld
     3307
     3308    DO il = 1, ncum                                                            ! cld
     3309      IF (i<=inb(il) .AND. nqcond(il,i)/=0 .AND. iflag(il)<=1) THEN            ! cld
     3310        qcond(il, i) = qcond(il, i)/nqcond(il, i)                              ! cld
     3311      END IF                                                                   ! cld
     3312    END DO
     3313
     3314!AC!      do j=1,ntra
     3315!AC!       do il=1,ncum
     3316!AC!        if (i.le.inb(il) .and. iflag(il) .le. 1) then
     3317!AC!         dpinv=1.0/(ph(il,i)-ph(il,i+1))
     3318!AC!         cpinv=1.0/cpn(il,i)
     3319!AC!
     3320!AC!         if (cvflag_grav) then
     3321!AC!          ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv
     3322!AC!     :     *(mp(il,i+1)*(trap(il,i+1,j)-tra(il,i,j))
     3323!AC!     :     -mp(il,i)*(trap(il,i,j)-trap(il,i-1,j)))
     3324!AC!         else
     3325!AC!          ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv
     3326!AC!     :     *(mp(il,i+1)*(trap(il,i+1,j)-tra(il,i,j))
     3327!AC!     :     -mp(il,i)*(trap(il,i,j)-trap(il,i-1,j)))
     3328!AC!         endif
     3329!AC!        endif ! i
     3330!AC!       enddo
     3331!AC!      enddo
    34173332
    34183333
    34193334500 END DO
    34203335
    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
     3336!JYG<
     3337!Conservation de l'eau
     3338!   sumdq = 0.
     3339!   DO k = 1, nl
     3340!     sumdq = sumdq + fr(1, k)*100.*(ph(1,k)-ph(1,k+1))/grav
     3341!   END DO
     3342!   PRINT *, 'cv3_yield, apres 500, sum(dq), precip, somme ', sumdq, Vprecip(1, 1), sumdq + vprecip(1, 1)
     3343!JYG>
     3344! ***   move the detrainment at level inb down to level inb-1   ***
     3345! ***        in such a way as to preserve the vertically        ***
     3346! ***          integrated enthalpy and water tendencies         ***
     3347
     3348! Correction bug le 18-03-09
    34273349  DO il = 1, ncum
    34283350    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
     3351      ax = 0.01*grav*ment(il, inb(il), inb(il))* &
     3352           (hp(il,inb(il))-h(il,inb(il))+t(il,inb(il))*(cpv-cpd)*(rr(il,inb(il))-qent(il,inb(il),inb(il))))/ &
     3353                                (cpn(il,inb(il))*(ph(il,inb(il))-ph(il,inb(il)+1)))
     3354      ft(il, inb(il)) = ft(il, inb(il)) - ax
     3355      ft(il, inb(il)-1) = ft(il, inb(il)-1) + ax*cpn(il, inb(il))*(ph(il,inb(il))-ph(il,inb(il)+1))/ &
     3356                              (cpn(il,inb(il)-1)*(ph(il,inb(il)-1)-ph(il,inb(il))))
     3357
     3358      bx = 0.01*grav*ment(il, inb(il), inb(il))*(qent(il,inb(il),inb(il))-rr(il,inb(il)))/ &
     3359                                                 (ph(il,inb(il))-ph(il,inb(il)+1))
     3360      fr(il, inb(il)) = fr(il, inb(il)) - bx
     3361      fr(il, inb(il)-1) = fr(il, inb(il)-1) + bx*(ph(il,inb(il))-ph(il,inb(il)+1))/ &
     3362                                                 (ph(il,inb(il)-1)-ph(il,inb(il)))
     3363
     3364      cx = 0.01*grav*ment(il, inb(il), inb(il))*(uent(il,inb(il),inb(il))-u(il,inb(il)))/ &
     3365                                                 (ph(il,inb(il))-ph(il,inb(il)+1))
     3366      fu(il, inb(il)) = fu(il, inb(il)) - cx
     3367      fu(il, inb(il)-1) = fu(il, inb(il)-1) + cx*(ph(il,inb(il))-ph(il,inb(il)+1))/ &
     3368                                                 (ph(il,inb(il)-1)-ph(il,inb(il)))
     3369
     3370      dx = 0.01*grav*ment(il, inb(il), inb(il))*(vent(il,inb(il),inb(il))-v(il,inb(il)))/ &
     3371                                                 (ph(il,inb(il))-ph(il,inb(il)+1))
     3372      fv(il, inb(il)) = fv(il, inb(il)) - dx
     3373      fv(il, inb(il)-1) = fv(il, inb(il)-1) + dx*(ph(il,inb(il))-ph(il,inb(il)+1))/ &
     3374                                                 (ph(il,inb(il)-1)-ph(il,inb(il)))
    34823375    END IF !iflag
    34833376  END DO
    34843377
    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    ***
     3378!JYG<
     3379!Conservation de l'eau
     3380!   sumdq = 0.
     3381!   DO k = 1, nl
     3382!     sumdq = sumdq + fr(1, k)*100.*(ph(1,k)-ph(1,k+1))/grav
     3383!   END DO
     3384!   PRINT *, 'cv3_yield, apres 503, sum(dq), precip, somme ', sumdq, Vprecip(1, 1), sumdq + vprecip(1, 1)
     3385!JYG>
     3386
     3387!AC!      do j=1,ntra
     3388!AC!       do il=1,ncum
     3389!AC!        IF (iflag(il) .le. 1) THEN
     3390!AC!    IF (cvflag_grav) then
     3391!AC!        ex=0.01*grav*ment(il,inb(il),inb(il))
     3392!AC!     :      *(traent(il,inb(il),inb(il),j)-tra(il,inb(il),j))
     3393!AC!     :      /(ph(i l,inb(il))-ph(il,inb(il)+1))
     3394!AC!        ftra(il,inb(il),j)=ftra(il,inb(il),j)-ex
     3395!AC!        ftra(il,inb(il)-1,j)=ftra(il,inb(il)-1,j)
     3396!AC!     :       +ex*(ph(il,inb(il))-ph(il,inb(il)+1))
     3397!AC!     :          /(ph(il,inb(il)-1)-ph(il,inb(il)))
     3398!AC!    else
     3399!AC!        ex=0.1*ment(il,inb(il),inb(il))
     3400!AC!     :      *(traent(il,inb(il),inb(il),j)-tra(il,inb(il),j))
     3401!AC!     :      /(ph(i l,inb(il))-ph(il,inb(il)+1))
     3402!AC!        ftra(il,inb(il),j)=ftra(il,inb(il),j)-ex
     3403!AC!        ftra(il,inb(il)-1,j)=ftra(il,inb(il)-1,j)
     3404!AC!     :       +ex*(ph(il,inb(il))-ph(il,inb(il)+1))
     3405!AC!     :          /(ph(il,inb(il)-1)-ph(il,inb(il)))
     3406!AC!        ENDIF   !cvflag grav
     3407!AC!        ENDIF    !iflag
     3408!AC!       enddo
     3409!AC!      enddo
     3410
     3411
     3412! ***    homogenize tendencies below cloud base    ***
    35113413
    35123414
    … …  
    35223424  END DO
    35233425
    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
     3426!do i=1,nl
     3427!do il=1,ncum
     3428!th_wake(il,i)=t_wake(il,i)*(1000.0/p(il,i))**rdcp
     3429!enddo
     3430!enddo
    35293431
    35303432  DO i = 1, nl
    35313433    DO il = 1, ncum
    35323434      IF (i<=(icb(il)-1) .AND. iflag(il)<=1) THEN
    3533         ! jyg  Saturated part : use T profile
     3435!jyg  Saturated part : use T profile
    35343436        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))
     3437!jyg<20140311
     3438!Correction pour conserver l eau
     3439        IF (ok_conserv_q) THEN
     3440          bsum(il) = bsum(il) + (fr(il,i)-fqd(il,i))*(ph(il,i)-ph(il,i+1))
     3441          csum(il) = csum(il) + (ph(il,i)-ph(il,i+1))
     3442
     3443        ELSE
     3444          bsum(il)=bsum(il)+(fr(il,i)-fqd(il,i))*(lv(il,i)+(cl-cpd)*(t(il,i)-t(il,1)))* &
     3445                            (ph(il,i)-ph(il,i+1))
     3446          csum(il)=csum(il)+(lv(il,i)+(cl-cpd)*(t(il,i)-t(il,1)))* &
     3447                            (ph(il,i)-ph(il,i+1))
     3448        ENDIF ! (ok_conserv_q)
     3449!jyg>
    35393450        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
     3451!jyg  Unsaturated part : use T_wake profile
    35413452        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)
     3453!jyg<20140311
     3454!Correction pour conserver l eau
     3455        IF (ok_conserv_q) THEN
     3456          fsum(il) = fsum(il) + fqd(il, i)*(ph(il,i)-ph(il,i+1))
     3457          gsum(il) = gsum(il) + (ph(il,i)-ph(il,i+1))
     3458        ELSE
     3459          fsum(il)=fsum(il)+fqd(il,i)*(lv(il,i)+(cl-cpd)*(t_wake(il,i)-t_wake(il,1)))* &
     3460                            (ph(il,i)-ph(il,i+1))
     3461          gsum(il)=gsum(il)+(lv(il,i)+(cl-cpd)*(t_wake(il,i)-t_wake(il,1)))* &
     3462                            (ph(il,i)-ph(il,i+1))
     3463        ENDIF ! (ok_conserv_q)
     3464!jyg>
     3465        hsum(il) = hsum(il) + t_wake(il, i)*(ph(il,i)-ph(il,i+1))/th_wake(il, i)
    35483466      END IF
    35493467    END DO
    35503468  END DO
    35513469
    3552   ! !!!      do 700 i=1,icb(il)-1
     3470!!!!      do 700 i=1,icb(il)-1
    35533471  DO i = 1, nl
    35543472    DO il = 1, ncum
    … …  
    35623480  END DO
    35633481
    3564 
    3565   ! ***   Check that moisture stays positive. If not, scale tendencies
    3566   ! in order to ensure moisture positivity
     3482!jyg<
     3483!Conservation de l'eau
     3484!!  sumdq = 0.
     3485!!  DO k = 1, nl
     3486!!    sumdq = sumdq + fr(1, k)*100.*(ph(1,k)-ph(1,k+1))/grav
     3487!!  END DO
     3488!!  PRINT *, 'cv3_yield, apres hom, sum(dq), precip, somme ', sumdq, Vprecip(1, 1), sumdq + vprecip(1, 1)
     3489!jyg>
     3490
     3491
     3492! ***   Check that moisture stays positive. If not, scale tendencies
     3493! in order to ensure moisture positivity
    35673494  DO il = 1, ncum
    35683495    alpha_qpos(il) = 1.
    35693496    IF (iflag(il)<=1) THEN
    35703497      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)))
     3498        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)))
    35733499      END IF
    35743500    END IF
    … …  
    35783504      IF (iflag(il)<=1) THEN
    35793505        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)
     3506          alpha_qpos1(il) = max(1., (-delt*fr(il,i))/(s_wake(il)*rr_wake(il,i)+(1.-s_wake(il))*rr(il,i)))
     3507          IF (alpha_qpos1(il)>=alpha_qpos(il)) alpha_qpos(il) = alpha_qpos1(il)
    35843508        END IF
    35853509      END IF
    … …  
    36083532        m(il, i) = m(il, i)/alpha_qpos(il)
    36093533        mp(il, i) = mp(il, i)/alpha_qpos(il)
    3610         vprecip(il, i) = vprecip(il, i)/alpha_qpos(il)
     3534        Vprecip(il, i) = vprecip(il, i)/alpha_qpos(il)
    36113535      END IF
    36123536    END DO
    … …  
    36223546  END DO
    36233547
    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           ***
     3548!AC!      DO j = 1,ntra
     3549!AC!      DO i = 1,nl
     3550!AC!       DO il = 1,ncum
     3551!AC!        IF (iflag(il) .le. 1) THEN
     3552!AC!         ftra(il,i,j) = ftra(il,i,j)/alpha_qpos(il)
     3553!AC!        ENDIF
     3554!AC!       ENDDO
     3555!AC!      ENDDO
     3556!AC!      ENDDO
     3557
     3558
     3559! ***           reset counter and return           ***
    36363560
    36373561  DO il = 1, ncum
    … …  
    37023626          END IF
    37033627        END IF
    3704         ! c        print *,'cbmf',il,i,k,cbmf(il),wghti(il,k)
     3628! c        print *,'cbmf',il,i,k,cbmf(il),wghti(il,k)
    37053629      END DO
    37063630    END DO
    … …  
    37103634    DO k = i, nl
    37113635      DO il = 1, ncum
    3712         ! test         if (i.ge.icb(il).and.i.le.inb(il).and.k.le.inb(il))
    3713         ! then
     3636! test         if (i.ge.icb(il).and.i.le.inb(il).and.k.le.inb(il)) then
    37143637        IF (i<=inb(il) .AND. k<=inb(il)) THEN
    37153638          upwd(il, i) = upwd(il, i) + up1(il, k, i)
    37163639          dnwd(il, i) = dnwd(il, i) + dn1(il, k, i)
    37173640        END IF
    3718         ! c         print
    3719         ! *,'upwd',il,i,k,inb(il),upwd(il,i),m(il,k),up1(il,k,i)
     3641! c         print *,'upwd',il,i,k,inb(il),upwd(il,i),m(il,k),up1(il,k,i)
    37203642      END DO
    37213643    END DO
    … …  
    37233645
    37243646
    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
     3647!!!!      DO il=1,ncum
     3648!!!!      do i=icb(il),inb(il)
     3649!!!!
     3650!!!!      upwd(il,i)=0.0
     3651!!!!      dnwd(il,i)=0.0
     3652!!!!      do k=i,inb(il)
     3653!!!!      up1=0.0
     3654!!!!      dn1=0.0
     3655!!!!      do n=1,i-1
     3656!!!!      up1=up1+ment(il,n,k)
     3657!!!!      dn1=dn1-ment(il,k,n)
     3658!!!!      enddo
     3659!!!!      upwd(il,i)=upwd(il,i)+m(il,k)+up1
     3660!!!!      dnwd(il,i)=dnwd(il,i)+dn1
     3661!!!!      enddo
     3662!!!!      enddo
     3663!!!!
     3664!!!!      ENDDO
     3665
     3666! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
     3667! determination de la variation de flux ascendant entre
     3668! deux niveau non dilue mip
     3669! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    37483670
    37493671  DO i = 1, nl
    … …  
    37873709  END DO
    37883710
    3789   ! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    3790   ! icb represente de niveau ou se trouve la
    3791   ! base du nuage , et inb le top du nuage
    3792   ! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
     3711! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
     3712! icb represente de niveau ou se trouve la
     3713! base du nuage , et inb le top du nuage
     3714! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    37933715
    37943716  DO i = 1, nd
    … …  
    38003722  DO i = 1, nd
    38013723    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)
     3724      rdcp = (rrd*(1.-rr(il,i))-rr(il,i)*rrv)/(cpd*(1.-rr(il,i))+rr(il,i)*cpv)
    38043725      tls(il, i) = t(il, i)*(1000.0/p(il,i))**rdcp
    38053726      tps(il, i) = tp(il, i)
    … …  
    38083729
    38093730
    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
     3731! *** diagnose the in-cloud mixing ratio   ***                       ! cld
     3732! ***           of condensed water         ***                       ! cld
     3733!! cld                                                               
     3734                                                                     
     3735  DO i = 1, nd                                                       ! cld
     3736    DO il = 1, ncum                                                  ! cld
     3737      mac(il, i) = 0.0                                               ! cld
     3738      wa(il, i) = 0.0                                                ! cld
     3739      siga(il, i) = 0.0                                              ! cld
     3740      sax(il, i) = 0.0                                               ! cld
     3741    END DO                                                           ! cld
     3742  END DO                                                             ! cld
     3743                                                                     
     3744  DO i = minorig, nl                                                 ! cld
     3745    DO k = i + 1, nl + 1                                             ! cld
     3746      DO il = 1, ncum                                                ! cld
    38263747        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
     3748          mac(il, i) = mac(il, i) + m(il, k)                         ! cld
     3749        END IF                                                       ! cld
     3750      END DO                                                         ! cld
     3751    END DO                                                           ! cld
     3752  END DO                                                             ! cld
     3753
     3754  DO i = 1, nl                                                       ! cld
     3755    DO j = 1, i                                                      ! cld
     3756      DO il = 1, ncum                                                ! cld
     3757        IF (i>=icb(il) .AND. i<=(inb(il)-1) &                        ! cld
     3758            .AND. j>=icb(il) .AND. iflag(il)<=1) THEN                ! cld
     3759          sax(il, i) = sax(il, i) + rrd*(tvp(il,j)-tv(il,j)) &       ! cld
     3760            *(ph(il,j)-ph(il,j+1))/p(il, j)                          ! cld
     3761        END IF                                                       ! cld
     3762      END DO                                                         ! cld
     3763    END DO                                                           ! cld
     3764  END DO                                                             ! cld
     3765
     3766  DO i = 1, nl                                                       ! cld
     3767    DO il = 1, ncum                                                  ! cld
     3768      IF (i>=icb(il) .AND. i<=(inb(il)-1) &                          ! cld
     3769          .AND. sax(il,i)>0.0 .AND. iflag(il)<=1) THEN               ! cld
     3770        wa(il, i) = sqrt(2.*sax(il,i))                               ! cld
     3771      END IF                                                         ! cld
     3772    END DO                                                           ! cld
     3773  END DO                                                             ! cld
     3774
     3775  DO i = 1, nl                                                       ! cld
     3776    DO il = 1, ncum                                                  ! cld
     3777      IF (wa(il,i)>0.0 .AND. iflag(il)<=1) &                         ! cld
     3778        siga(il, i) = mac(il, i)/wa(il, i) &                         ! cld
     3779        *rrd*tvp(il, i)/p(il, i)/100./delta                          ! cld
     3780      siga(il, i) = min(siga(il,i), 1.0)                             ! cld
     3781! IM cf. FH                                                         
     3782      IF (iflag_clw==0) THEN                                         ! cld
     3783        qcondc(il, i) = siga(il, i)*clw(il, i)*(1.-ep(il,i)) &       ! cld
     3784          +(1.-siga(il,i))*qcond(il, i)                              ! cld
     3785      ELSE IF (iflag_clw==1) THEN                                    ! cld
     3786        qcondc(il, i) = qcond(il, i)                                 ! cld
     3787      END IF                                                         ! cld
     3788
     3789    END DO                                                           ! cld
     3790  END DO
     3791! print*,'cv3_yield fin'
     3792
    38723793  RETURN
    38733794END SUBROUTINE cv3_yield
    38743795
    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)
     3796!AC! et !RomP >>>
     3797SUBROUTINE cv3_tracer(nloc, len, ncum, nd, na, &
     3798                      ment, sigij, da, phi, phi2, d1a, dam, &
     3799                      ep, Vprecip, elij, clw, epmlmMm, eplaMm, &
     3800                      icb, inb)
    38783801  IMPLICIT NONE
    38793802
    38803803  include "cv3param.h"
    38813804
    3882   ! inputs:
     3805!inputs:
    38833806  INTEGER ncum, nd, na, nloc, len
    38843807  REAL ment(nloc, na, na), sigij(nloc, na, na)
    … …  
    38863809  REAL ep(nloc, na)
    38873810  INTEGER icb(nloc), inb(nloc)
    3888   REAL vprecip(nloc, nd+1)
    3889   ! ouputs:
     3811  REAL Vprecip(nloc, nd+1)
     3812!ouputs:
    38903813  REAL da(nloc, na), phi(nloc, na, na)
    38913814  REAL phi2(nloc, na, na)
    38923815  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:
     3816  REAL epmlmMm(nloc, na, na), eplaMm(nloc, na)
     3817! variables pour tracer dans precip de l'AA et des mel
     3818!local variables:
    38963819  INTEGER i, j, k
    38973820  REAL epm(nloc, na, na)
    38983821
    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
     3822! variables d'Emanuel : du second indice au troisieme
     3823! --->    tab(i,k,j) -> de l origine k a l arrivee j
     3824! ment, sigij, elij
     3825! variables personnelles : du troisieme au second indice
     3826! --->    tab(i,j,k) -> de k a j
     3827! phi, phi2
     3828
     3829! initialisations
    39073830
    39083831  da(:, :) = 0.
    … …  
    39103833  dam(:, :) = 0.
    39113834  epm(:, :, :) = 0.
    3912   eplamm(:, :) = 0.
    3913   epmlmmm(:, :, :) = 0.
     3835  eplaMm(:, :) = 0.
     3836  epmlmMm(:, :, :) = 0.
    39143837  phi(:, :, :) = 0.
    39153838  phi2(:, :, :) = 0.
    39163839
    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
     3840! fraction deau condensee dans les melanges convertie en precip : epm
     3841! et eau condensée précipitée dans masse d'air saturé : l_m*dM_m/dzdz.dzdz
    39193842  DO j = 1, na
    39203843    DO k = 1, na
    39213844      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)
     3845        IF (k>=icb(i) .AND. k<=inb(i) .AND. &
     3846!!jyg              j.ge.k.and.j.le.inb(i)) then
     3847!!jyg             epm(i,j,k)=1.-(1.-ep(i,j))*clw(i,j)/elij(i,k,j)
    39253848            j>k .AND. j<=inb(i)) THEN
    39263849          epm(i, j, k) = 1. - (1.-ep(i,j))*clw(i, j)/max(elij(i,k,j), 1.E-16)
    3927           ! !
     3850!!
    39283851          epm(i, j, k) = max(epm(i,j,k), 0.0)
    39293852        END IF
    … …  
    39373860      DO i = 1, ncum
    39383861        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))
     3862          eplaMm(i, j) = eplamm(i, j) + &
     3863                         ep(i, j)*clw(i, j)*ment(i, j, k)*(1.-sigij(i,j,k))
    39413864        END IF
    39423865      END DO
    … …  
    39483871      DO i = 1, ncum
    39493872        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)
     3873          epmlmMm(i, j, k) = epm(i, j, k)*elij(i, k, j)*ment(i, k, j)
    39513874        END IF
    39523875      END DO
    … …  
    39543877  END DO
    39553878
    3956   ! matrices pour calculer la tendance des concentrations dans cvltr.F90
     3879! matrices pour calculer la tendance des concentrations dans cvltr.F90
    39573880  DO j = 1, na
    39583881    DO k = 1, na
    … …  
    39623885        d1a(i, j) = d1a(i, j) + ment(i, k, j)*ep(i, k)*(1.-sigij(i,k,j))
    39633886        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 
     3887          dam(i, j) = dam(i, j) + ment(i, k, j)*epm(i, k, j)*(1.-ep(i,k))*(1.-sigij(i,k,j))
    39673888          phi2(i, j, k) = phi(i, j, k)*epm(i, j, k)
    39683889        END IF
    … …  
    39733894  RETURN
    39743895END 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)
     3896!AC! et !RomP <<<
     3897
     3898SUBROUTINE cv3_uncompress(nloc, len, ncum, nd, ntra, idcum, &
     3899                          iflag, &
     3900                          precip, sig, w0, &
     3901                          ft, fq, fu, fv, ftra, &
     3902                          Ma, upwd, dnwd, dnwd0, qcondc, wd, cape, &
     3903                          iflag1, &
     3904                          precip1, sig1, w01, &
     3905                          ft1, fq1, fu1, fv1, ftra1, &
     3906                          Ma1, upwd1, dnwd1, dnwd01, qcondc1, wd1, cape1)
    39813907  IMPLICIT NONE
    39823908
    39833909  include "cv3param.h"
    39843910
    3985   ! inputs:
     3911!inputs:
    39863912  INTEGER len, ncum, nd, ntra, nloc
    39873913  INTEGER idcum(nloc)
    … …  
    39963922  REAL wd(nloc), cape(nloc)
    39973923
    3998   ! outputs:
     3924!outputs:
    39993925  INTEGER iflag1(len)
    40003926  REAL precip1(len)
    … …  
    40073933  REAL wd1(nloc), cape1(nloc)
    40083934
    4009   ! local variables:
     3935!local variables:
    40103936  INTEGER i, k, j
    40113937
    … …  
    40383964
    40393965
    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
     3966!AC!        do 2100 j=1,ntra
     3967!AC!c oct3         do 2110 k=1,nl
     3968!AC!         do 2110 k=1,nd ! oct3
     3969!AC!          do 2120 i=1,ncum
     3970!AC!            ftra1(idcum(i),k,j)=ftra(i,k,j)
     3971!AC! 2120     continue
     3972!AC! 2110    continue
     3973!AC! 2100   continue
     3974!
    40483975  RETURN
    40493976END SUBROUTINE cv3_uncompress
Note: See TracChangeset for help on using the changeset viewer.