Ignore:
Timestamp:
Nov 21, 2019, 4:43:45 PM (5 years ago)
Author:
lguez
Message:

Merge revisions 3427:3600 of trunk into branch Ocean_skin

Location:
LMDZ6/branches/Ocean_skin
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/branches/Ocean_skin

  • LMDZ6/branches/Ocean_skin/libf/phylmd/cv3_routines.F90

    r3345 r3605  
    3535
    3636  include "cv3param.h"
     37  include "cvflag.h"
    3738  include "conema3.h"
    3839
    … …  
    125126     tlcrit=-55.0
    126127     CALL getin_p('tlcrit',tlcrit)
     128     ejectliq=0.
     129     CALL getin_p('ejectliq',ejectliq)
     130     ejectice=0.
     131     CALL getin_p('ejectice',ejectice)
     132     cvflag_prec_eject = .FALSE.
     133     CALL getin_p('cvflag_prec_eject',cvflag_prec_eject)
     134     qsat_depends_on_qt = .FALSE.
     135     CALL getin_p('qsat_depends_on_qt',qsat_depends_on_qt)
     136     adiab_ascent_mass_flux_depends_on_ejectliq = .FALSE.
     137     CALL getin_p('adiab_ascent_mass_flux_depends_on_ejectliq',adiab_ascent_mass_flux_depends_on_ejectliq)
     138     keepbug_ice_frac = .TRUE.
     139     CALL getin_p('keepbug_ice_frac', keepbug_ice_frac)
    127140
    128141    WRITE (*, *) 't_top_max=', t_top_max
    … …  
    144157    WRITE (*, *) 'elcrit=', elcrit
    145158    WRITE (*, *) 'tlcrit=', tlcrit
     159    WRITE (*, *) 'ejectliq=', ejectliq
     160    WRITE (*, *) 'ejectice=', ejectice
     161    WRITE (*, *) 'cvflag_prec_eject =', cvflag_prec_eject
     162    WRITE (*, *) 'qsat_depends_on_qt =', qsat_depends_on_qt
     163    WRITE (*, *) 'adiab_ascent_mass_flux_depends_on_ejectliq =', adiab_ascent_mass_flux_depends_on_ejectliq
     164    WRITE (*, *) 'keepbug_ice_frac =', keepbug_ice_frac
     165
    146166    first = .FALSE.
    147167  END IF ! (first)
    … …  
    170190
    171191  include "cv3param.h"
     192  include "cvflag.h"
    172193
    173194!inputs:
    … …  
    236257! debug          lv(i,k)= lv0-clmcpv*(t(i,k)-t0)
    237258      lv(i, k) = lv0 - clmcpv*(t(i,k)-273.15)
    238       lf(i, k) = lf0 - clmci*(t(i,k)-273.15)
     259!!      lf(i, k) = lf0 - clmci*(t(i,k)-273.15)   ! erreur de signe !!
     260      lf(i, k) = lf0 + clmci*(t(i,k)-273.15)
    239261      cpn(i, k) = cpd*(1.0-q(i,k)) + cpv*q(i, k)
    240262      cpx(i, k) = cpd*(1.0-q(i,k)) + cl*q(i, k)
    … …  
    289311  USE mod_phys_lmdz_transfert_para, ONLY : bcast
    290312  USE add_phys_tend_mod, ONLY: fl_cor_ebil
     313  USE print_control_mod, ONLY: prt_level
    291314  IMPLICIT NONE
    292315
    … …  
    516539    END DO
    517540  ENDIF
     541  IF (prt_level .GE. 10) THEN
     542    print *,'cv3_feed : iflag(1), pfeed(1), plcl(1), wghti(1,k) ', &
     543                        iflag(1), pfeed(1), plcl(1), (wghti(1,k),k=1,10)
     544  ENDIF
    518545
    519546! -------------------------------------------------------------------
    … …  
    11051132                         tnk, qnk, gznk, hnk, t, q, qs, gz, &
    11061133                         p, ph, h, tv, lv, lf, pbase, buoybase, plcl, &
    1107                          inb, tp, tvp, clw, hp, ep, sigp, buoy, frac)
     1134                         inb, tp, tvp, clw, hp, ep, sigp, buoy, &
     1135                         frac_a, frac_s, qpreca, qta)
    11081136  USE print_control_mod, ONLY: prt_level
    11091137  IMPLICIT NONE
    … …  
    11531181  REAL, DIMENSION (nloc, nd), INTENT (OUT)           :: ep, sigp, hp
    11541182  REAL, DIMENSION (nloc, nd), INTENT (OUT)           :: buoy
    1155   REAL, DIMENSION (nloc, nd), INTENT (OUT)           :: frac
     1183  REAL, DIMENSION (nloc, nd), INTENT (OUT)           :: frac_a, frac_s
     1184  REAL, DIMENSION (nloc, nd), INTENT (OUT)           :: qpreca
     1185  REAL, DIMENSION (nloc, nd), INTENT (OUT)           :: qta
    11561186
    11571187!local variables:
    11581188  INTEGER i, j, k
    1159   REAL tg, qg, ahg, alv, alf, s, tc, es, esi, denom, rg, tca, elacrit
     1189  REAL smallestreal
     1190  REAL tg, qg, dqgdT, ahg, alv, alf, s, tc, es, esi, denom, rg, tca, elacrit
     1191  REAL                                               :: phinu2p
    11601192  REAL als
    1161   REAL qsat_new, snew, qi(nloc, nd)
    1162   REAL by, defrac, pden, tbis
    1163   REAL ah0(nloc), cape(nloc), capem(nloc), byp(nloc)
    1164   LOGICAL lcape(nloc)
    1165   INTEGER iposit(nloc)
    1166   REAL fracg
    1167   REAL deltap
     1193  REAL                                               :: qsat_new, snew
     1194  REAL, DIMENSION (nloc,nd)                          :: qi
     1195  REAL, DIMENSION (nloc,nd)                          :: ha    ! moist static energy of adiabatic ascents
     1196                                                              ! taking into account precip ejection
     1197  REAL, DIMENSION (nloc,nd)                          :: hla   ! liquid water static energy of adiabatic ascents
     1198                                                              ! taking into account precip ejection
     1199  REAL, DIMENSION (nloc,nd)                          :: qcld  ! specific cloud water
     1200  REAL, DIMENSION (nloc,nd)                          :: qhsat    ! specific humidity at saturation
     1201  REAL, DIMENSION (nloc,nd)                          :: dqhsatdT ! dqhsat/dT
     1202  REAL, DIMENSION (nloc,nd)                          :: frac  ! ice fraction function of envt temperature
     1203  REAL, DIMENSION (nloc,nd)                          :: qps   ! specific solid precipitation
     1204  REAL, DIMENSION (nloc,nd)                          :: qpl   ! specific liquid precipitation
     1205  REAL, DIMENSION (nloc)                             :: ah0, cape, capem, byp
     1206  LOGICAL, DIMENSION (nloc)                          :: lcape
     1207  INTEGER, DIMENSION (nloc)                          :: iposit
     1208  REAL                                               :: denomm1
     1209  REAL                                               :: by, defrac, pden, tbis
     1210  REAL                                               :: fracg
     1211  REAL                                               :: deltap
     1212  REAL, SAVE                                         :: Tx, Tm
     1213  DATA Tx/263.15/, Tm/243.15/
     1214!$OMP THREADPRIVATE(Tx, Tm)
     1215  REAL                                               :: aa, bb, dd, ddelta, discr
     1216  REAL                                               :: ff, fp
     1217  REAL                                               :: coefx, coefm, Zx, Zm, Ux, U, Um
    11681218
    11691219  IF (prt_level >= 10) THEN
    1170     print *,'cv3_undilute2.0. t(1,k), q(1,k), qs(1,k) ', &
    1171                         (k, t(1,k), q(1,k), qs(1,k), k = 1,nl)
     1220    print *,'cv3_undilute2.0. icvflag_Tpa, t(1,k), q(1,k), qs(1,k) ', &
     1221                        icvflag_Tpa, (k, t(1,k), q(1,k), qs(1,k), k = 1,nl)
    11721222  ENDIF
     1223  smallestreal=tiny(smallestreal)
    11731224
    11741225! =====================================================================
    … …  
    11811232    END DO
    11821233  END DO
     1234
    11831235
    11841236! =====================================================================
    … …  
    11971249             qnk(i)*(lv0-clmcpv*(tnk(i)-273.15)) + gznk(i)
    11981250  END DO
    1199 
     1251!
     1252!  Ice fraction
     1253!
     1254  IF (cvflag_ice) THEN
     1255    DO k = minorig, nl
     1256      DO i = 1, ncum
     1257          frac(i, k) = (Tx - t(i,k))/(Tx - Tm)
     1258          frac(i, k) = min(max(frac(i,k),0.0), 1.0)
     1259      END DO
     1260    END DO
     1261! Below cloud base, set ice fraction to cloud base value
     1262    DO k = 1, nl
     1263      DO i = 1, ncum
     1264        IF (k<icb(i)) THEN
     1265          frac(i,k) = frac(i,icb(i))
     1266        END IF
     1267      END DO
     1268    END DO
     1269  ELSE
     1270    DO k = 1, nl
     1271      DO i = 1, ncum
     1272          frac(i,k) = 0.
     1273      END DO
     1274    END DO
     1275  ENDIF ! (cvflag_ice)
     1276
     1277
     1278  DO k = minorig, nl
     1279    DO i = 1,ncum
     1280      ha(i,k) = ah0(i)
     1281      hla(i,k) = hnk(i)
     1282      qta(i,k) = qnk(i)
     1283      qpreca(i,k) = 0.
     1284      frac_a(i,k) = 0.
     1285      frac_s(i,k) = frac(i,k)
     1286      qpl(i,k) = 0.
     1287      qps(i,k) = 0.
     1288      qhsat(i,k) = qs(i,k)
     1289      qcld(i,k) = max(qta(i,k)-qhsat(i,k),0.)
     1290      IF (k <= icb(i)+1) THEN
     1291        qhsat(i,k) = qnk(i)-clw(i,k)
     1292        qcld(i,k) = clw(i,k)
     1293      ENDIF
     1294    ENDDO
     1295  ENDDO
     1296
     1297!jyg<
     1298! =====================================================================
     1299! --- SET THE THE FRACTION OF PRECIPITATION FALLING OUTSIDE OF CLOUD
     1300! =====================================================================
     1301  DO k = 1, nl
     1302    DO i = 1, ncum
     1303      ep(i, k) = 0.0
     1304      sigp(i, k) = spfac
     1305    END DO
     1306  END DO
     1307!>jyg
     1308!
    12001309
    12011310! ***  Find lifted parcel quantities above cloud base    ***
    12021311
    1203 
     1312!----------------------------------------------------------------------------
     1313!
     1314  IF (icvflag_Tpa == 2) THEN
     1315!
     1316!----------------------------------------------------------------------------
     1317!
     1318    DO k = minorig + 1, nl
     1319      DO i = 1,ncum
     1320        tp(i,k) = t(i,k)
     1321      ENDDO
     1322!!      alv = lv0 - clmcpv*(t(i,k)-273.15)
     1323!!      alf = lf0 + clmci*(t(i,k)-273.15)
     1324!!      als = alf + alv
     1325      DO j = 1,4
     1326        DO i = 1, ncum
     1327! ori       if(k.ge.(icb(i)+1))then
     1328          IF (k>=(icbs(i)+1)) THEN                                ! convect3
     1329            tg = tp(i, k)
     1330            IF (tg .gt. Tx) THEN
     1331              es = 6.112*exp(17.67*(tg - 273.15)/(tg + 243.5 - 273.15))
     1332              qg = eps*es/(p(i,k)-es*(1.-eps))
     1333            ELSE
     1334              esi = exp(23.33086-(6111.72784/tg)+0.15215*log(tg))
     1335              qg = eps*esi/(p(i,k)-esi*(1.-eps))
     1336            ENDIF
     1337! Ice fraction
     1338            ff = 0.
     1339            fp = 1./(Tx - Tm)
     1340            IF (tg < Tx) THEN
     1341              IF (tg > Tm) THEN
     1342                ff = (Tx - tg)*fp
     1343              ELSE
     1344                ff = 1.
     1345              ENDIF ! (tg > Tm)
     1346            ENDIF ! (tg < Tx)
     1347! Intermediate variables
     1348            aa = cpd + (cl-cpd)*qnk(i) + lv(i,k)*lv(i,k)*qg/(rrv*tg*tg)
     1349            ahg = (cpd + (cl-cpd)*qnk(i))*tg + lv(i,k)*qg - &
     1350                  lf(i,k)*ff*(qnk(i) - qg) + gz(i,k)
     1351            dd = lf(i,k)*lv(i,k)*qg/(rrv*tg*tg)
     1352            ddelta = lf(i,k)*(qnk(i) - qg)
     1353            bb = aa + ddelta*fp + dd*fp*(Tx-tg)
     1354! Compute Zx and Zm
     1355            coefx = aa
     1356            coefm = aa + dd
     1357            IF (tg .gt. Tx) THEN
     1358              Zx = ahg            + coefx*(Tx - tg)
     1359              Zm = ahg - ddelta   + coefm*(Tm - tg)
     1360            ELSE
     1361              IF (tg .gt. Tm) THEN
     1362                Zx = ahg          + (coefx +fp*ddelta)*(Tx - Tg)
     1363                Zm = ahg          + (coefm +fp*ddelta)*(Tm - Tg)
     1364              ELSE
     1365                Zx = ahg + ddelta + coefx*(Tx - tg)
     1366                Zm = ahg          + coefm*(Tm - tg)
     1367              ENDIF ! (tg .gt. Tm)
     1368            ENDIF ! (tg .gt. Tx)
     1369! Compute the masks Um, U, Ux
     1370            Um = (sign(1., Zm-ah0(i))+1.)/2.
     1371            Ux = (sign(1., ah0(i)-Zx)+1.)/2.
     1372            U = (1. - Um)*(1. - Ux)
     1373! Compute the updated parcell temperature Tp : 3 cases depending on tg value
     1374            IF (tg .gt. Tx) THEN
     1375              discr = bb*bb - 4*dd*fp*(ah0(i) - ahg + ddelta*fp*(Tx-tg))
     1376              Tp(i,k) = tg + &
     1377                  Um*  (ah0(i) - ahg + ddelta)           /(aa + dd) + &
     1378                  U *2*(ah0(i) - ahg + ddelta*fp*(Tx-tg))/(bb + sqrt(discr)) + &
     1379                  Ux*  (ah0(i) - ahg)                    /aa
     1380            ELSEIF (tg .gt. Tm) THEN
     1381              discr = bb*bb - 4*dd*fp*(ah0(i) - ahg)
     1382              Tp(i,k) = tg + &
     1383                  Um*  (ah0(i) - ahg + ddelta*fp*(tg-Tm))/(aa + dd) + &
     1384                  U *2*(ah0(i) - ahg)                    /(bb + sqrt(discr)) + &
     1385                  Ux*  (ah0(i) - ahg + ddelta*fp*(tg-Tx))/aa
     1386            ELSE
     1387              discr = bb*bb - 4*dd*fp*(ah0(i) - ahg + ddelta*fp*(Tm-tg))
     1388              Tp(i,k) = tg + &
     1389                  Um*  (ah0(i) - ahg)                    /(aa + dd) + &
     1390                  U *2*(ah0(i) - ahg + ddelta*fp*(Tm-tg))/(bb + sqrt(discr)) + &
     1391                  Ux*  (ah0(i) - ahg - ddelta)           /aa
     1392            ENDIF ! (tg .gt. Tx)
     1393!
     1394!!     print *,' j, k, Um, U, Ux, aa, bb, discr, dd, ddelta ', j, k, Um, U, Ux, aa, bb, discr, dd, ddelta
     1395!!     print *,' j, k, ah0(i), ahg, tg, qg, tp(i,k), ff ', j, k, ah0(i), ahg, tg, qg, tp(i,k), ff
     1396          END IF ! (k>=(icbs(i)+1))
     1397        END DO ! i = 1, ncum
     1398      END DO ! j = 1,4
     1399      DO i = 1, ncum
     1400        IF (k>=(icbs(i)+1)) THEN                                ! convect3
     1401          tg = tp(i, k)
     1402          IF (tg .gt. Tx) THEN
     1403            es = 6.112*exp(17.67*(tg - 273.15)/(tg + 243.5 - 273.15))
     1404            qg = eps*es/(p(i,k)-es*(1.-eps))
     1405          ELSE
     1406            esi = exp(23.33086-(6111.72784/tg)+0.15215*log(tg))
     1407            qg = eps*esi/(p(i,k)-esi*(1.-eps))
     1408          ENDIF
     1409          clw(i, k) = qnk(i) - qg
     1410          clw(i, k) = max(0.0, clw(i,k))
     1411          tvp(i, k) = max(0., tp(i,k)*(1.+qg/eps-qnk(i)))
     1412! print*,tvp(i,k),'tvp'
     1413          IF (clw(i,k)<1.E-11) THEN
     1414            tp(i, k) = tv(i, k)
     1415            tvp(i, k) = tv(i, k)
     1416          END IF ! (clw(i,k)<1.E-11)
     1417        END IF ! (k>=(icbs(i)+1))
     1418      END DO ! i = 1, ncum
     1419    END DO ! k = minorig + 1, nl
     1420!----------------------------------------------------------------------------
     1421!
     1422  ELSE IF (icvflag_Tpa == 1) THEN  ! (icvflag_Tpa == 2)
     1423!
     1424!----------------------------------------------------------------------------
     1425!
     1426    DO k = minorig + 1, nl
     1427      DO i = 1,ncum
     1428        tp(i,k) = t(i,k)
     1429      ENDDO
     1430!!      alv = lv0 - clmcpv*(t(i,k)-273.15)
     1431!!      alf = lf0 + clmci*(t(i,k)-273.15)
     1432!!      als = alf + alv
     1433      DO j = 1,4
     1434        DO i = 1, ncum
     1435! ori       if(k.ge.(icb(i)+1))then
     1436          IF (k>=(icbs(i)+1)) THEN                                ! convect3
     1437            tg = tp(i, k)
     1438            IF (tg .gt. Tx .OR. .NOT.cvflag_ice) THEN
     1439              es = 6.112*exp(17.67*(tg - 273.15)/(tg + 243.5 - 273.15))
     1440              qg = eps*es/(p(i,k)-es*(1.-eps))
     1441              dqgdT = lv(i,k)*qg/(rrv*tg*tg)
     1442            ELSE
     1443              esi = exp(23.33086-(6111.72784/tg)+0.15215*log(tg))
     1444              qg = eps*esi/(p(i,k)-esi*(1.-eps))
     1445              dqgdT = (lv(i,k)+lf(i,k))*qg/(rrv*tg*tg)
     1446            ENDIF
     1447            IF (qsat_depends_on_qt) THEN
     1448              dqgdT = dqgdT*(1.-qta(i,k-1))/(1.-qg)**2
     1449              qg = qg*(1.-qta(i,k-1))/(1.-qg)           
     1450            ENDIF
     1451            ahg = (cpd + (cl-cpd)*qta(i,k-1))*tg + lv(i,k)*qg - &
     1452                  lf(i,k)*frac(i,k)*(qta(i,k-1) - qg) + gz(i,k)
     1453            Tp(i,k) = tg + (ah0(i) - ahg)/ &
     1454                    (cpd + (cl-cpd)*qta(i,k-1) + (lv(i,k)+frac(i,k)*lf(i,k))*dqgdT)
     1455!!   print *,'undilute2 iterations k, Tp(i,k), ah0(i), ahg ', &
     1456!!                                 k, Tp(i,k), ah0(i), ahg
     1457          END IF ! (k>=(icbs(i)+1))
     1458        END DO ! i = 1, ncum
     1459      END DO ! j = 1,4
     1460      DO i = 1, ncum
     1461        IF (k>=(icbs(i)+1)) THEN                                ! convect3
     1462          tg = tp(i, k)
     1463          IF (tg .gt. Tx .OR. .NOT.cvflag_ice) THEN
     1464            es = 6.112*exp(17.67*(tg - 273.15)/(tg + 243.5 - 273.15))
     1465            qg = eps*es/(p(i,k)-es*(1.-eps))
     1466          ELSE
     1467            esi = exp(23.33086-(6111.72784/tg)+0.15215*log(tg))
     1468            qg = eps*esi/(p(i,k)-esi*(1.-eps))
     1469          ENDIF
     1470          IF (qsat_depends_on_qt) THEN
     1471            qg = qg*(1.-qta(i,k-1))/(1.-qg)           
     1472          ENDIF
     1473          qhsat(i,k) = qg
     1474        END IF ! (k>=(icbs(i)+1))
     1475      END DO ! i = 1, ncum
     1476      DO i = 1, ncum
     1477        IF (k>=(icbs(i)+1)) THEN                                ! convect3
     1478          clw(i, k) = qta(i,k-1) - qhsat(i,k)
     1479          clw(i, k) = max(0.0, clw(i,k))
     1480          tvp(i, k) = max(0., tp(i,k)*(1.+qhsat(i,k)/eps-qta(i,k-1)))
     1481! print*,tvp(i,k),'tvp'
     1482          IF (clw(i,k)<1.E-11) THEN
     1483            tp(i, k) = tv(i, k)
     1484            tvp(i, k) = tv(i, k)
     1485          END IF ! (clw(i,k)<1.E-11)
     1486        END IF ! (k>=(icbs(i)+1))
     1487      END DO ! i = 1, ncum
     1488!
     1489      IF (cvflag_prec_eject) THEN
     1490        DO i = 1, ncum
     1491          IF (k>=(icbs(i)+1)) THEN                                ! convect3
     1492!  Specific precipitation (liquid and solid) and ice content
     1493!  before ejection of precipitation                                                     !!jygprl
     1494            elacrit = elcrit*min(max(1.-(tp(i,k)-T0)/Tlcrit, 0.), 1.)                   !!jygprl
     1495!!!!            qcld(i,k) = min(clw(i,k), elacrit)                                          !!jygprl
     1496            qcld(i,k) = min(clw(i,k), elacrit*(1.-qta(i,k-1))/(1.-elacrit))             !!jygprl
     1497            phinu2p = qhsat(i,k-1) + qcld(i,k-1) - (qhsat(i,k) + qcld(i,k))                     !!jygprl
     1498            qpl(i,k) = qpl(i,k-1) + (1.-frac(i,k))*phinu2p                            !!jygprl
     1499            qps(i,k) = qps(i,k-1) + frac(i,k)     *phinu2p                            !!jygprl
     1500            qi(i,k) = (1.-ejectliq)*clw(i,k)*frac(i,k) + &                            !!jygprl
     1501                     ejectliq*(qps(i,k-1) + frac(i,k)*(phinu2p+qcld(i,k)))            !!jygprl
     1502!!
     1503!  =====================================================================================
     1504!  Ejection of precipitation from adiabatic ascents if requested (cvflag_prec_eject=True):
     1505!  Compute the steps of total water (qta), of moist static energy (ha), of specific
     1506!  precipitation (qpl and qps) and of specific cloud water (qcld) associated with precipitation
     1507!   ejection.
     1508!  =====================================================================================
     1509! 
     1510!   Verif
     1511            qpreca(i,k) = ejectliq*qpl(i,k) + ejectice*qps(i,k)                                   !!jygprl
     1512            frac_a(i,k) = ejectice*qps(i,k)/max(qpreca(i,k),smallestreal)                         !!jygprl
     1513            frac_s(i,k) = (1.-ejectliq)*frac(i,k) + &                                             !!jygprl
     1514               ejectliq*(1. - (qpl(i,k)+(1.-frac(i,k))*qcld(i,k))/max(clw(i,k),smallestreal))     !!jygprl
     1515!         
     1516            denomm1 = 1./(1. - qpreca(i,k))
     1517!         
     1518            qta(i,k) = qta(i,k-1) - &
     1519                      qpreca(i,k)*(1.-qta(i,k-1))*denomm1
     1520            ha(i,k)  = ha(i,k-1) + &
     1521                      ( qpreca(i,k)*(-(1.-qta(i,k-1))*(cl-cpd)*tp(i,k) + &
     1522                                  lv(i,k)*qhsat(i,k) - lf(i,k)*(frac_s(i,k)*qcld(i,k)+qps(i,k))) + &
     1523                        lf(i,k)*ejectice*qps(i,k))*denomm1
     1524            hla(i,k) = hla(i,k-1) + &
     1525                      ( qpreca(i,k)*(-(1.-qta(i,k-1))*(cpv-cpd)*tp(i,k) - &
     1526                                  lv(i,k)*((1.-frac_s(i,k))*qcld(i,k)+qpl(i,k)) - &
     1527                                  (lv(i,k)+lf(i,k))*(frac_s(i,k)*qcld(i,k)+qps(i,k))) + &
     1528                        lv(i,k)*ejectliq*qpl(i,k) + (lv(i,k)+lf(i,k))*ejectice*qps(i,k))*denomm1
     1529            qpl(i,k) = qpl(i,k)*(1.-ejectliq)*denomm1
     1530            qps(i,k) = qps(i,k)*(1.-ejectice)*denomm1
     1531            qcld(i,k) = qcld(i,k)*denomm1
     1532            qhsat(i,k) = qhsat(i,k)*(1.-qta(i,k))/(1.-qta(i,k-1))
     1533         END IF ! (k>=(icbs(i)+1))
     1534        END DO ! i = 1, ncum
     1535      ENDIF  ! (cvflag_prec_eject)
     1536!
     1537    END DO ! k = minorig + 1, nl
     1538!
     1539!----------------------------------------------------------------------------
     1540!
     1541  ELSE IF (icvflag_Tpa == 0) THEN! (icvflag_Tpa == 2) ELSE IF(icvflag_Tpa == 1)
     1542!
     1543!----------------------------------------------------------------------------
     1544!
    12041545  DO k = minorig + 1, nl
    12051546    DO i = 1, ncum
    … …  
    13581699  END DO
    13591700
    1360   IF (prt_level >= 10) THEN
    1361     print *,'cv3_undilute2.1. tp(1,k), tvp(1,k) ', &
    1362                           (k, tp(1,k), tvp(1,k), k = 1,nl)
    1363   ENDIF
    1364 
     1701!----------------------------------------------------------------------------
     1702!
     1703  ENDIF ! (icvflag_Tpa == 2) ELSEIF (icvflag_Tpa == 1) ELSE (icvflag_Tpa == 0)
     1704!
     1705!----------------------------------------------------------------------------
     1706!
    13651707! =====================================================================
    1366 ! --- SET THE PRECIPITATION EFFICIENCIES AND THE FRACTION OF
    1367 ! --- PRECIPITATION FALLING OUTSIDE OF CLOUD
     1708! --- SET THE PRECIPITATION EFFICIENCIES
    13681709! --- THESE MAY BE FUNCTIONS OF TP(I), P(I) AND CLW(I)
    13691710! =====================================================================
    1370 !
    1371 !jyg<
    1372   DO k = 1, nl
    1373     DO i = 1, ncum
    1374       ep(i, k) = 0.0
    1375       sigp(i, k) = spfac
    1376     END DO
    1377   END DO
    1378 !>jyg
    13791711!
    13801712  IF (flag_epkeorig/=1) THEN
    … …  
    14131745    END DO
    14141746  END IF
     1747!
     1748!   =========================================================================
     1749  IF (prt_level >= 10) THEN
     1750    print *,'cv3_undilute2.1. tp(1,k), tvp(1,k) ', &
     1751                          (k, tp(1,k), tvp(1,k), k = 1,nl)
     1752  ENDIF
    14151753!
    14161754! =====================================================================
    … …  
    16481986  IF (cvflag_ice) THEN
    16491987!
     1988  IF (cvflag_prec_eject) THEN
     1989!!    DO k = minorig + 1, nl
     1990!!      DO i = 1, ncum
     1991!!        IF ((k>=icb(i)) .AND. (k<=inb(i))) THEN
     1992!!          frac_s(i,k) = qi(i,k)/max(clw(i,k),smallestreal)   
     1993!!          frac_s(i,k) = 1. - (qpl(i,k)+(1.-frac_s(i,k))*qcld(i,k))/max(clw(i,k),smallestreal)   
     1994!!        END IF
     1995!!      END DO
     1996!!    END DO
     1997  ELSE    ! (cvflag_prec_eject)
    16501998    DO k = minorig + 1, nl
    16511999      DO i = 1, ncum
    16522000        IF ((k>=icb(i)) .AND. (k<=inb(i))) THEN
    1653           frac(i, k) = 1. - (t(i,k)-243.15)/(263.15-243.15)
    1654           frac(i, k) = min(max(frac(i,k),0.0), 1.0)
    1655           hp(i, k) = hnk(i) + (lv(i,k)+(cpd-cpv)*t(i,k)+frac(i,k)*lf(i,k))* &
    1656                               ep(i, k)*clw(i, k)
     2001!jyg< frac computation moved to beginning of cv3_undilute2.
     2002!     kept here for compatibility test with CMip6 version
     2003          frac_s(i, k) = 1. - (t(i,k)-243.15)/(263.15-243.15)
     2004          frac_s(i, k) = min(max(frac_s(i,k),0.0), 1.0)
    16572005        END IF
    16582006      END DO
    16592007    END DO
    1660 ! Below cloud base, set ice fraction to cloud base value
    1661     DO k = 1, nl
     2008  ENDIF  ! (cvflag_prec_eject) ELSE
     2009    DO k = minorig + 1, nl
    16622010      DO i = 1, ncum
    1663         IF (k<icb(i)) THEN
    1664           frac(i,k) = frac(i,icb(i))
     2011        IF ((k>=icb(i)) .AND. (k<=inb(i))) THEN
     2012!!          hp(i, k) = hnk(i) + (lv(i,k)+(cpd-cpv)*t(i,k)+frac_s(i,k)*lf(i,k))* &     !!jygprl
     2013!!                              ep(i, k)*clw(i, k)                                    !!jygprl
     2014          hp(i, k) = hla(i,k-1) + (lv(i,k)+(cpd-cpv)*t(i,k)+frac_s(i,k)*lf(i,k))* &   !!jygprl
     2015                              ep(i, k)*clw(i, k)                                      !!jygprl
    16652016        END IF
    16662017      END DO
    16672018    END DO
    16682019!
    1669   ELSE
     2020  ELSE   ! (cvflag_ice)
    16702021!
    16712022    DO k = minorig + 1, nl
    … …  
    23502701SUBROUTINE cv3_unsat(nloc, ncum, nd, na, ntra, icb, inb, iflag, &
    23512702                     t, rr, rs, gz, u, v, tra, p, ph, &
    2352                      th, tv, lv, lf, cpn, ep, sigp, clw, &
     2703                     th, tv, lv, lf, cpn, ep, sigp, clw, frac_s, qpreca, frac_a, qta , &                       !!jygprl
    23532704                     m, ment, elij, delt, plcl, coef_clos, &
    23542705                     mp, rp, up, vp, trap, wt, water, evap, fondue, ice, &
    23552706                     faci, b, sigd, &
    2356                      wdtrainA, wdtrainM)                                      ! RomP
     2707                     wdtrainA, wdtrainS, wdtrainM)                                      ! RomP
    23572708  USE print_control_mod, ONLY: prt_level, lunout
    23582709  IMPLICIT NONE
    … …  
    23722723  REAL, DIMENSION (nloc, na), INTENT (IN)            :: gz
    23732724  REAL, DIMENSION (nloc, nd), INTENT (IN)            :: u, v
    2374   REAL tra(nloc, nd, ntra)
    2375   REAL p(nloc, nd), ph(nloc, nd+1)
    2376   REAL, DIMENSION (nloc, na), INTENT (IN)            :: ep, sigp, clw
     2725  REAL, DIMENSION (nloc, nd, ntra), INTENT(IN)       :: tra
     2726  REAL, DIMENSION (nloc, nd), INTENT (IN)            :: p
     2727  REAL, DIMENSION (nloc, nd+1), INTENT (IN)          :: ph
     2728  REAL, DIMENSION (nloc, na), INTENT (IN)            :: ep, sigp, clw   !adiab ascent shedding
     2729  REAL, DIMENSION (nloc, na), INTENT (IN)            :: frac_s          !ice fraction in adiab ascent shedding !!jygprl
     2730  REAL, DIMENSION (nloc, na), INTENT (IN)            :: qpreca          !adiab ascent precip                   !!jygprl
     2731  REAL, DIMENSION (nloc, na), INTENT (IN)            :: frac_a          !ice fraction in adiab ascent precip   !!jygprl
     2732  REAL, DIMENSION (nloc, na), INTENT (IN)            :: qta             !adiab ascent specific total water     !!jygprl
    23772733  REAL, DIMENSION (nloc, na), INTENT (IN)            :: th, tv, lv, cpn
    23782734  REAL, DIMENSION (nloc, na), INTENT (IN)            :: lf
    … …  
    23872743  REAL, DIMENSION (nloc, na), INTENT (OUT)           :: mp, rp, up, vp
    23882744  REAL, DIMENSION (nloc, na), INTENT (OUT)           :: water, evap, wt
    2389   REAL, DIMENSION (nloc, na), INTENT (OUT)           :: ice, fondue, faci
     2745  REAL, DIMENSION (nloc, na), INTENT (OUT)           :: ice, fondue
     2746  REAL, DIMENSION (nloc, na), INTENT (OUT)           :: faci            ! ice fraction in precipitation
    23902747  REAL, DIMENSION (nloc, na, ntra), INTENT (OUT)     :: trap
    23912748  REAL, DIMENSION (nloc, na), INTENT (OUT)           :: b
    … …  
    23952752! Distinction des wdtrain
    23962753! Pa = wdtrainA     Pm = wdtrainM
    2397   REAL, DIMENSION (nloc, na), INTENT (OUT)           :: wdtrainA, wdtrainM
     2754  REAL, DIMENSION (nloc, na), INTENT (OUT)           :: wdtrainA, wdtrainS, wdtrainM
    23982755
    23992756!local variables
    24002757  INTEGER i, j, k, il, num1, ndp1
     2758  REAL smallestreal
    24012759  REAL tinv, delti, coef
    24022760  REAL awat, afac, afac1, afac2, bfac
    … …  
    24052763  REAL ampmax, thaw
    24062764  REAL tevap(nloc)
    2407   REAL lvcp(nloc, na), lfcp(nloc, na)
    2408   REAL h(nloc, na), hm(nloc, na)
    2409   REAL frac(nloc, na)
    2410   REAL fraci(nloc, na), prec(nloc, na)
     2765  REAL, DIMENSION (nloc, na)      :: lvcp, lfcp
     2766  REAL, DIMENSION (nloc, na)      :: h, hm
     2767  REAL, DIMENSION (nloc, na)      :: ma
     2768  REAL, DIMENSION (nloc, na)      :: frac          ! ice fraction in precipitation source
     2769  REAL, DIMENSION (nloc, na)      :: fraci         ! provisionnal ice fraction in precipitation
     2770  REAL, DIMENSION (nloc, na)      :: prec
    24112771  REAL wdtrain(nloc)
    24122772  LOGICAL lwork(nloc), mplus(nloc)
    … …  
    24152775! ------------------------------------------------------
    24162776IF (prt_level .GE. 10) print *,' ->cv3_unsat, iflag(1) ', iflag(1)
     2777
     2778smallestreal=tiny(smallestreal)
    24172779
    24182780! =============================
    … …  
    24342796!! RomP >>>
    24352797wdtrainA(:,:) = 0.
     2798wdtrainS(:,:) = 0.
    24362799wdtrainM(:,:) = 0.
    24372800!! RomP <<<
    … …  
    24892852  END DO
    24902853
     2854!
     2855! Get adiabatic ascent mass flux
     2856!
     2857!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     2858  IF (adiab_ascent_mass_flux_depends_on_ejectliq) THEN
     2859!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     2860!!! Warning : this option leads to water conservation violation
     2861!!!           Expert only
     2862!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     2863    DO il = 1, ncum
     2864      ma(il, nlp) = 0.
     2865      ma(il, 1)   = 0.
     2866    END DO
     2867
     2868  DO i = nl, 2, -1
     2869      DO il = 1, ncum
     2870        ma(il, i) = ma(il, i+1)*(1.-qta(il,i))/(1.-qta(il,i-1)) + m(il, i)
     2871      END DO
     2872  END DO
     2873!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     2874  ELSE ! (adiab_ascent_mass_flux_depends_on_ejectliq)
     2875!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     2876    DO il = 1, ncum
     2877      ma(il, nlp) = 0.
     2878      ma(il, 1)   = 0.
     2879    END DO
     2880
     2881  DO i = nl, 2, -1
     2882      DO il = 1, ncum
     2883        ma(il, i) = ma(il, i+1) + m(il, i)
     2884      END DO
     2885  END DO
     2886
     2887  ENDIF ! (adiab_ascent_mass_flux_depends_on_ejectliq) ELSE
     2888!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    24912889
    24922890! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    … …  
    25132911! ***              calculate detrained precipitation             ***
    25142912
    2515     DO il = 1, ncum
    2516       IF (i<=inb(il) .AND. lwork(il)) THEN
    2517         IF (cvflag_grav) THEN
    2518           wdtrain(il) = grav*ep(il, i)*m(il, i)*clw(il, i)
    2519           wdtrainA(il, i) = wdtrain(il)/grav                        !   Pa   RomP
    2520         ELSE
    2521           wdtrain(il) = 10.0*ep(il, i)*m(il, i)*clw(il, i)
    2522           wdtrainA(il, i) = wdtrain(il)/10.                         !   Pa   RomP
    2523         END IF
    2524       END IF
    2525     END DO
     2913
     2914    DO il = 1, ncum                                                   
     2915      IF (i<=inb(il) .AND. lwork(il)) THEN                           
     2916        wdtrain(il) = grav*ep(il, i)*m(il, i)*clw(il, i)           
     2917        wdtrainS(il, i) = wdtrain(il)/grav                                            !   Ps   jyg
     2918!!        wdtrainA(il, i) = wdtrain(il)/grav                                          !   Ps   RomP
     2919      END IF                                                         
     2920    END DO                                                           
    25262921
    25272922    IF (i>1) THEN
    … …  
    25312926            awat = elij(il, j, i) - (1.-ep(il,i))*clw(il, i)
    25322927            awat = max(awat, 0.0)
    2533             IF (cvflag_grav) THEN
    2534               wdtrain(il) = wdtrain(il) + grav*awat*ment(il, j, i)
    2535               wdtrainM(il, i) = wdtrain(il)/grav - wdtrainA(il, i)  !   Pm  RomP
    2536             ELSE
    2537               wdtrain(il) = wdtrain(il) + 10.0*awat*ment(il, j, i)
    2538               wdtrainM(il, i) = wdtrain(il)/10. - wdtrainA(il, i)   !   Pm  RomP
    2539             END IF
     2928            wdtrain(il) = wdtrain(il) + grav*awat*ment(il, j, i)
     2929            wdtrainM(il, i) = wdtrain(il)/grav - wdtrainS(il, i)    !   Pm  jyg
     2930!!            wdtrainM(il, i) = wdtrain(il)/grav - wdtrainA(il, i)  !   Pm  RomP
    25402931          END IF
    25412932        END DO
    … …  
    25432934    END IF
    25442935
     2936    IF (cvflag_prec_eject) THEN
     2937!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     2938      IF (adiab_ascent_mass_flux_depends_on_ejectliq) THEN
     2939!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     2940!!! Warning : this option leads to water conservation violation
     2941!!!           Expert only
     2942!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     2943          IF ( i > 1) THEN
     2944            DO il = 1, ncum
     2945              IF (i<=inb(il) .AND. lwork(il)) THEN
     2946                wdtrainA(il,i) = ma(il, i+1)*(qta(il, i-1)-qta(il,i))/(1. - qta(il, i-1))    !   Pa   jygprl
     2947                wdtrain(il) = wdtrain(il) + grav*wdtrainA(il,i)
     2948              END IF
     2949            END DO
     2950          ENDIF  ! ( i > 1)
     2951!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     2952      ELSE ! (adiab_ascent_mass_flux_depends_on_ejectliq)
     2953!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     2954          IF ( i > 1) THEN
     2955            DO il = 1, ncum
     2956              IF (i<=inb(il) .AND. lwork(il)) THEN
     2957                wdtrainA(il,i) = ma(il, i+1)*(qta(il, i-1)-qta(il,i))                        !   Pa   jygprl
     2958                wdtrain(il) = wdtrain(il) + grav*wdtrainA(il,i)
     2959              END IF
     2960            END DO
     2961          ENDIF  ! ( i > 1)
     2962
     2963      ENDIF ! (adiab_ascent_mass_flux_depends_on_ejectliq) ELSE
     2964!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     2965    ENDIF  ! (cvflag_prec_eject)
     2966
    25452967
    25462968! ***    find rain water and evaporation using provisional   ***
    … …  
    25482970
    25492971
     2972    IF (cvflag_ice) THEN                                                                                !!jygprl
     2973      IF (cvflag_prec_eject) THEN
     2974        DO il = 1, ncum                                                                                   !!jygprl
     2975          IF (i<=inb(il) .AND. lwork(il)) THEN                                                            !!jygprl
     2976            frac(il, i) = (frac_a(il,i)*wdtrainA(il,i)+frac_s(il,i)*(wdtrainS(il,i)+wdtrainM(il,i))) / &  !!jygprl
     2977                          max(wdtrainA(il,i)+wdtrainS(il,i)+wdtrainM(il,i),smallestreal)                  !!jygprl
     2978            fraci(il, i) = frac(il, i)                                                                    !!jygprl
     2979          END IF                                                                                          !!jygprl
     2980        END DO                                                                                            !!jygprl
     2981      ELSE  ! (cvflag_prec_eject)
     2982        DO il = 1, ncum                                                                                   !!jygprl
     2983          IF (i<=inb(il) .AND. lwork(il)) THEN                                                            !!jygprl
     2984!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     2985            IF (keepbug_ice_frac) THEN
     2986              frac(il, i) = frac_s(il, i)
     2987!       Ice fraction computed again here as a function of the temperature seen by unsaturated downdraughts
     2988!       (i.e. the cold pool temperature) for compatibility with earlier versions.
     2989              fraci(il, i) = 1. - (t(il,i)-243.15)/(263.15-243.15)
     2990              fraci(il, i) = min(max(fraci(il,i),0.0), 1.0)
     2991!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     2992            ELSE  ! (keepbug_ice_frac)
     2993!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     2994              frac(il, i) = frac_s(il, i)
     2995              fraci(il, i) = frac(il, i)                                                                    !!jygprl
     2996            ENDIF  ! (keepbug_ice_frac)
     2997!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     2998          END IF                                                                                          !!jygprl
     2999        END DO                                                                                            !!jygprl
     3000      ENDIF  ! (cvflag_prec_eject)
     3001    END IF                                                                                              !!jygprl
     3002
     3003
    25503004    DO il = 1, ncum
    25513005      IF (i<=inb(il) .AND. lwork(il)) THEN
    … …  
    25533007        wt(il, i) = 45.0
    25543008
    2555         IF (cvflag_ice) THEN
    2556           frac(il, inb(il)) = 1. - (t(il,inb(il))-243.15)/(263.15-243.15)
    2557           frac(il, inb(il)) = min(max(frac(il,inb(il)),0.), 1.)
    2558           fraci(il, inb(il)) = frac(il, inb(il))
    2559         ELSE
    2560           CONTINUE
    2561         END IF
    2562 
    25633009        IF (i<inb(il)) THEN
    2564 
    2565           IF (cvflag_ice) THEN
    2566 !CR:tmax_fonte_cv: T for which ice is totally melted (used to be 275.15)
    2567             thaw = (t(il,i)-273.15)/(tmax_fonte_cv-273.15)
    2568             thaw = min(max(thaw,0.0), 1.0)
    2569             frac(il, i) = frac(il, i)*(1.-thaw)
    2570           ELSE
    2571             CONTINUE
    2572           END IF
    2573 
    25743010          rp(il, i) = rp(il, i+1) + &
    25753011                      (cpd*(t(il,i+1)-t(il,i))+gz(il,i+1)-gz(il,i))/lv(il, i)
    25763012          rp(il, i) = 0.5*(rp(il,i)+rr(il,i))
    25773013        END IF
    2578         fraci(il, i) = 1. - (t(il,i)-243.15)/(263.15-243.15)
    2579         fraci(il, i) = min(max(fraci(il,i),0.0), 1.0)
    25803014        rp(il, i) = max(rp(il,i), 0.0)
    25813015        rp(il, i) = amin1(rp(il,i), rs(il,i))
    … …  
    29983432
    29993433  RETURN
     3434
    30003435END SUBROUTINE cv3_unsat
    30013436
    … …  
    30043439                     t, rr, t_wake, rr_wake, s_wake, u, v, tra, &
    30053440                     gz, p, ph, h, hp, lv, lf, cpn, th, th_wake, &
    3006                      ep, clw, m, tp, mp, rp, up, vp, trap, &
     3441                     ep, clw, qpreca, m, tp, mp, rp, up, vp, trap, &
    30073442                     wt, water, ice, evap, fondue, faci, b, sigd, &
    30083443                     ment, qent, hent, iflag_mix, uent, vent, &
    … …  
    30143449!!                     tls, tps,                             ! useless . jyg
    30153450                     qcondc, wd, &
    3016                      ftd, fqd, qnk, qtc, sigt, tau_cld_cv, coefw_cld_cv)
     3451                     ftd, fqd, qta, qtc, sigt, tau_cld_cv, coefw_cld_cv)
    30173452
    30183453    USE print_control_mod, ONLY: lunout, prt_level
    … …  
    30543489      REAL, DIMENSION (nloc, na, na, ntra), INTENT (IN)  :: traent
    30553490      REAL, DIMENSION (nloc, nd), INTENT (IN)            :: tv, tvp, wghti
    3056       REAL,INTENT(IN)                                    :: tau_cld_cv, coefw_cld_cv
     3491      REAL, DIMENSION (nloc, nd), INTENT (IN)            :: qta
     3492      REAL, DIMENSION (nloc, na),INTENT(IN)              :: qpreca
     3493      REAL, INTENT(IN)                                   :: tau_cld_cv, coefw_cld_cv
    30573494!
    30583495!input/output:
    … …  
    30833520      REAL                                               :: ax, bx, cx, dx, ex
    30843521      REAL                                               :: cpinv, rdcp, dpinv
     3522      REAL                                               :: sigaq
    30853523      REAL, DIMENSION (nloc)                             ::  awat
    30863524      REAL, DIMENSION (nloc, nd)                         :: lvcp, lfcp              ! , mke ! unused . jyg
    … …  
    31003538      REAL, DIMENSION (nloc)                             :: sument
    31013539      REAL, DIMENSION (nloc, nd)                         :: sigment, qtment             ! cld
    3102       REAL, DIMENSION (nloc)                             :: qnk
    31033540      REAL sumdq !jyg
    31043541!
    … …  
    32113648  END DO
    32123649
     3650! - Adiabatic ascent mass flux "ma" and cloud base mass flux "cbmf"
     3651!-----------------------------------------------------------------
     3652!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     3653  IF (adiab_ascent_mass_flux_depends_on_ejectliq) THEN
     3654!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     3655!!! Warning : this option leads to water conservation violation
     3656!!!           Expert only
     3657!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     3658  DO il = 1, ncum
     3659    ma(il, nlp) = 0.
     3660    ma(il, 1)   = 0.
     3661  END DO
     3662  DO k = nl, 2, -1
     3663    DO il = 1, ncum
     3664      ma(il, k) = ma(il, k+1)*(1.-qta(il, k))/(1.-qta(il, k-1)) + m(il, k)
     3665      cbmf(il) = max(cbmf(il), ma(il,k))
     3666    END DO
     3667  END DO
     3668  DO k = 2,nl
     3669    DO il = 1, ncum
     3670      IF (k <icb(il)) THEN
     3671        ma(il, k) = ma(il, k-1) + wghti(il,k-1)*cbmf(il)
     3672      ENDIF
     3673    END DO
     3674  END DO
     3675!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     3676  ELSE ! (adiab_ascent_mass_flux_depends_on_ejectliq)
     3677!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     3678!! Line kept for compatibility with earlier versions
    32133679  DO k = 2, nl
    32143680    DO il = 1, ncum
    … …  
    32193685  END DO
    32203686
     3687  DO il = 1, ncum
     3688    ma(il, nlp) = 0.
     3689    ma(il, 1)   = 0.
     3690  END DO
     3691  DO k = nl, 2, -1
     3692    DO il = 1, ncum
     3693      ma(il, k) = ma(il, k+1) + m(il, k)
     3694    END DO
     3695  END DO
     3696  DO k = 2,nl
     3697    DO il = 1, ncum
     3698      IF (k <icb(il)) THEN
     3699        ma(il, k) = ma(il, k-1) + wghti(il,k-1)*cbmf(il)
     3700      ENDIF
     3701    END DO
     3702  END DO
     3703
     3704  ENDIF ! (adiab_ascent_mass_flux_depends_on_ejectliq) ELSE
     3705!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     3706!
    32213707!    print*,'cv3_yield avant ft'
    32223708! am is the part of cbmf taken from the first level
    … …  
    33553841!***    Compute convective mass fluxes upwd and dnwd      ***
    33563842
     3843!
     3844! =================================================
     3845!              upward fluxes                      |
     3846! ------------------------------------------------
     3847!
    33573848upwd(:,:) = 0.
    33583849up_to(:,:) = 0.
    33593850up_from(:,:) = 0.
    3360 dnwd(:,:) = 0.
    3361 dn_to(:,:) = 0.
    3362 dn_from(:,:) = 0.
    3363 !
    3364 ! =================================================
    3365 !              upward fluxes                      |
    3366 ! ------------------------------------------------
     3851!
     3852!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     3853  IF (adiab_ascent_mass_flux_depends_on_ejectliq) THEN
     3854!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     3855!! The decrease of the adiabatic ascent mass flux due to ejection of precipitation
     3856!! is taken into account.
     3857!! WARNING : in the present version, taking into account the mass-flux decrease due to
     3858!! precipitation ejection leads to water conservation violation.
     3859!
     3860! - Upward mass flux of mixed draughts
     3861!---------------------------------------
     3862DO i = 2, nl
     3863  DO j = 1, i-1
     3864    DO il = 1, ncum
     3865      IF (i<=inb(il)) THEN
     3866        up_to(il,i) = up_to(il,i) + ment(il,j,i)
     3867      ENDIF
     3868    ENDDO
     3869  ENDDO
     3870ENDDO
     3871!
     3872DO j = 3, nl
     3873  DO i = 2, j-1
     3874    DO il = 1, ncum
     3875      IF (j<=inb(il)) THEN
     3876        up_from(il,i) = up_from(il,i) + ment(il,i,j)
     3877      ENDIF
     3878    ENDDO
     3879  ENDDO
     3880ENDDO
     3881!
     3882! The difference between upwd(il,i) and upwd(il,i-1) is due to updrafts ending in layer
     3883!(i-1) (theses drafts cross interface (i-1) but not interface(i)) and to updrafts starting
     3884!from layer (i-1) (theses drafts cross interface (i) but not interface(i-1)):
     3885!
     3886DO i = 2, nlp
     3887  DO il = 1, ncum
     3888    IF (i<=inb(il)+1) THEN
     3889      upwd(il,i) = max(0., upwd(il,i-1) - up_to(il,i-1) + up_from(il,i-1))
     3890    ENDIF
     3891  ENDDO
     3892ENDDO
     3893!
     3894! - Total upward mass flux
     3895!---------------------------
     3896DO i = 2, nlp
     3897  DO il = 1, ncum
     3898    IF (i<=inb(il)+1) THEN
     3899      upwd(il,i) = upwd(il,i) + ma(il,i)
     3900    ENDIF
     3901  ENDDO
     3902ENDDO
     3903!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     3904  ELSE ! (adiab_ascent_mass_flux_depends_on_ejectliq)
     3905!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     3906!! The decrease of the adiabatic ascent mass flux due to ejection of precipitation
     3907!! is not taken into account.
     3908!
     3909! - Upward mass flux
     3910!-------------------
    33673911DO i = 2, nl
    33683912  DO il = 1, ncum
    … …  
    33873931  ENDDO
    33883932ENDDO
    3389 !!DO i = 2, nl
    3390 !!  DO j = i+1, nl          !! Permuter les boucles i et j
     3933!
    33913934DO j = 3, nl
    33923935  DO i = 2, j-1
    … …  
    34103953  ENDDO
    34113954ENDDO
     3955
     3956
     3957  ENDIF ! (adiab_ascent_mass_flux_depends_on_ejectliq) ELSE
     3958!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     3959
    34123960!
    34133961! =================================================
    34143962!              downward fluxes                    |
    34153963! ------------------------------------------------
     3964dnwd(:,:) = 0.
     3965dn_to(:,:) = 0.
     3966dn_from(:,:) = 0.
    34163967DO i = 1, nl
    34173968  DO j = i+1, nl
    … …  
    34243975ENDDO
    34253976!
    3426 !!DO i = 2, nl
    3427 !!  DO j = 1, i-1          !! Permuter les boucles i et j
    34283977DO j = 1, nl
    34293978  DO i = j+1, nl
    … …  
    37494298    END DO ! cld
    37504299
     4300!ym BIG Warning : it seems that the k loop is missing !!!
     4301!ym Strong advice to check this
     4302!ym add a k loop temporary
     4303
    37514304! (particular case: no detraining level is found)                              ! cld
     4305! Verif merge Dynamico<<<<<<< .working
    37524306    DO il = 1, ncum                                                            ! cld
    37534307      IF (i<=inb(il) .AND. nent(il,i)==0 .AND. iflag(il)<=1) THEN              ! cld
    … …  
    37614315      END IF                                                                   ! cld
    37624316    END DO                                                                     ! cld
     4317! Verif merge Dynamico =======
     4318! Verif merge Dynamico     DO k = i + 1, nl
     4319! Verif merge Dynamico       DO il = 1, ncum        !ym k loop added                                    ! cld
     4320! Verif merge Dynamico         IF (i<=inb(il) .AND. nent(il,i)==0 .AND. iflag(il)<=1) THEN              ! cld
     4321! Verif merge Dynamico           qcond(il, i) = qcond(il, i) + (1.-ep(il,i))*clw(il, i)                 ! cld
     4322! Verif merge Dynamico           qtment(il, i) = qent(il,k,i) + qtment(il,i)                          ! cld
     4323! Verif merge Dynamico           nqcond(il, i) = nqcond(il, i) + 1.                                     ! cld
     4324! Verif merge Dynamico         END IF                                                                   ! cld
     4325! Verif merge Dynamico       END DO
     4326! Verif merge Dynamico     ENDDO                                                                     ! cld
     4327! Verif merge Dynamico >>>>>>> .merge-right.r3413
    37634328
    37644329    DO il = 1, ncum                                                            ! cld
    … …  
    41814746!!!!
    41824747!!!!      ENDDO
     4748
     4749!!  DO i = 1, nlp
     4750!!    DO il = 1, ncum
     4751!!      ma(il, i) = 0
     4752!!    END DO
     4753!!  END DO
     4754!!
     4755!!  DO i = 1, nl
     4756!!    DO j = i, nl
     4757!!      DO il = 1, ncum
     4758!!        ma(il, i) = ma(il, i) + m(il, j)
     4759!!      END DO
     4760!!    END DO
     4761!!  END DO
     4762
     4763!jyg<  (loops stop at nl)
     4764!!  DO i = nl + 1, nd
     4765!!    DO il = 1, ncum
     4766!!      ma(il, i) = 0.
     4767!!    END DO
     4768!!  END DO
     4769!>jyg
     4770
     4771!!  DO i = 1, nl
     4772!!    DO il = 1, ncum
     4773!!      IF (i<=(icb(il)-1)) THEN
     4774!!        ma(il, i) = 0
     4775!!      END IF
     4776!!    END DO
     4777!!  END DO
     4778
    41834779!-----------------------------------------------------------
    41844780        ENDIF !(.NOT.ok_optim_yield)                      !|
    … …  
    42054801!>jyg
    42064802
    4207   DO i = 1, nlp
    4208     DO il = 1, ncum
    4209       ma(il, i) = 0
    4210     END DO
    4211   END DO
    4212 
    4213   DO i = 1, nl
    4214     DO j = i, nl
    4215       DO il = 1, ncum
    4216         ma(il, i) = ma(il, i) + m(il, j)
    4217       END DO
    4218     END DO
    4219   END DO
    4220 
    4221 !jyg<  (loops stop at nl)
    4222 !!  DO i = nl + 1, nd
    4223 !!    DO il = 1, ncum
    4224 !!      ma(il, i) = 0.
    4225 !!    END DO
    4226 !!  END DO
    4227 !>jyg
    4228 
    4229   DO i = 1, nl
    4230     DO il = 1, ncum
    4231       IF (i<=(icb(il)-1)) THEN
    4232         ma(il, i) = 0
    4233       END IF
    4234     END DO
    4235   END DO
    42364803
    42374804! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    … …  
    43204887! 14/01/15 AJ delta n'a rien à faire là...                                                 
    43214888    DO il = 1, ncum                                                  ! cld
    4322       IF (wa(il,i)>0.0 .AND. iflag(il)<=1) &                         ! cld
     4889!!      IF (wa(il,i)>0.0 .AND. iflag(il)<=1) &                         ! cld
     4890!!        siga(il, i) = mac(il, i)/(coefw_cld_cv*wa(il, i)) &          ! cld
     4891!!        *rrd*tvp(il, i)/p(il, i)/100.                                ! cld
     4892!!
     4893!!      siga(il, i) = min(siga(il,i), 1.0)                             ! cld
     4894      sigaq = 0.
     4895      IF (wa(il,i)>0.0 .AND. iflag(il)<=1)  THEN                     ! cld
    43234896        siga(il, i) = mac(il, i)/(coefw_cld_cv*wa(il, i)) &          ! cld
    4324         *rrd*tvp(il, i)/p(il, i)/100.                                ! cld
    4325 
    4326       siga(il, i) = min(siga(il,i), 1.0)                             ! cld
     4897                     *rrd*tvp(il, i)/p(il, i)/100.                   ! cld
     4898        siga(il, i) = min(siga(il,i), 1.0)                           ! cld
     4899        sigaq = siga(il,i)*qta(il,i-1)                               ! cld
     4900      ENDIF
    43274901
    43284902! IM cf. FH
    … …  
    43364910        sigment(il,i)=sument(il)*tau_cld_cv/(ph(il,i)-ph(il,i+1))    ! cld
    43374911        sigment(il, i) = min(1.e-4+sigment(il,i), 1.0 - siga(il,i))  ! cld
    4338         qtc(il, i) = (siga(il,i)*qnk(il)+sigment(il,i)*qtment(il,i)) & ! cld
     4912!!        qtc(il, i) = (siga(il,i)*qta(il,i-1)+sigment(il,i)*qtment(il,i)) & ! cld
     4913        qtc(il, i) = (sigaq+sigment(il,i)*qtment(il,i)) & ! cld
    43394914                     /(siga(il,i)+sigment(il,i))                     ! cld
    43404915        sigt(il,i) = sigment(il, i) + siga(il, i)
    43414916
    4342 !        qtc(il, i) = siga(il,i)*qnk(il)+(1.-siga(il,i))*qtment(il,i) ! cld
     4917!        qtc(il, i) = siga(il,i)*qta(il,i-1)+(1.-siga(il,i))*qtment(il,i) ! cld
    43434918!     print*,'BIGAUSSIAN CONV',siga(il,i),sigment(il,i),qtc(il,i) 
    43444919               
    … …  
    46295204      do k=1,nl
    46305205        do i=1,ncum
    4631           hp(i,k)=h(i,k)
    4632         enddo
     5206          hp(i,k)=h(i,k)
     5207        enddo
    46335208      enddo
    46345209
Note: See TracChangeset for help on using the changeset viewer.