Changeset 4910


Ignore:
Timestamp:
Apr 18, 2024, 6:16:25 PM (2 weeks ago)
Author:
evignon
Message:

debut du chantier de refonte de cloudth_mpc

Location:
LMDZ6/trunk/libf/phylmd
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/trunk/libf/phylmd/lmdz_cloudth.F90

    r4674 r4910  
    14971497
    14981498!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    1499 SUBROUTINE cloudth_mpc(klon,klev,ind2,iflag_mpc_bl,mpc_bl_points,           &
    1500 &           temp,ztv,po,zqta,fraca,zpspsk,paprs,pplay,ztla,zthl,            &
    1501 &           ratqs,zqs,snowflux,qcloud,qincloud,icefrac,ctot,ctot_vol,       &
     1499SUBROUTINE cloudth_mpc(klon,klev,ind2,mpc_bl_points,                        &
     1500&           temp,qt,qt_th,frac_th,zpspsk,paprsup, paprsdn,play,thetal_th,   &
     1501&           ratqs,qcloud,qincloud,icefrac,ctot,ctot_vol,                    &
    15021502&           cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
    15031503!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    1504 ! Author : Arnaud Octavio Jam (LMD/CNRS), Etienne Vignon (LMDZ/CNRS)
    1505 ! Date: Adapted from cloudth_vert_v3 in 2021
    1506 ! Aim : computes qc and rneb in thermals with cold microphysical considerations
    1507 !       + for mixed phase boundary layer clouds, calculate ql and qi from
    1508 !       a stationary MPC model
    1509 ! IMPORTANT NOTE: we assume iflag_clouth_vert=3
     1504! Author : Etienne Vignon (LMDZ/CNRS)
     1505! Date: April 2024
     1506! Date: Adapted from cloudth_vert_v3 in 2023 by Arnaud Otavio Jam
     1507! IMPORTANT NOTE: we assume iflag_cloudth_vert=7
    15101508!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    15111509
    15121510      use lmdz_cloudth_ini, only: iflag_cloudth_vert,iflag_ratqs
    1513       use lmdz_cloudth_ini, only:  C_mpc ,Ni,C_cap,Ei,d_top ,vert_alpha, vert_alpha_th ,sigma1s_factor,sigma1s_power,sigma2s_factor,sigma2s_power,cloudth_ratqsmin ,iflag_cloudth_vert_noratqs
    1514       use lmdz_lscp_tools, ONLY: CALC_QSAT_ECMWF, FALLICE_VELOCITY
    1515 
    1516       use phys_local_var_mod, ONLY : qlth, qith, qsith, wiceth
     1511      use lmdz_cloudth_ini, only: C_mpc ,Ni,C_cap,Ei,d_top ,vert_alpha, vert_alpha_th ,sigma1s_factor,sigma1s_power,sigma2s_factor,sigma2s_power,cloudth_ratqsmin ,iflag_cloudth_vert_noratqs
     1512      use lmdz_lscp_tools,  only: CALC_QSAT_ECMWF
     1513      use lmdz_lscp_ini,    only: RTT, RG, RPI, RD, RCPD, RLVTT, RLSTT, temp_nowater, min_frac_th_cld, min_neb_th
    15171514
    15181515      IMPLICIT NONE
    15191516
    1520       INCLUDE "YOMCST.h"
    1521       INCLUDE "YOETHF.h"
    1522       INCLUDE "FCTTRE.h"
    1523      
    15241517
    15251518!------------------------------------------------------------------------------
    … …  
    15321525     
    15331526
    1534       REAL, DIMENSION(klon,klev), INTENT(IN)      ::  temp          ! Temperature [K]
    1535       REAL, DIMENSION(klon,klev), INTENT(IN)      ::  ztv           ! Virtual potential temp [K]
    1536       REAL, DIMENSION(klon),      INTENT(IN)      ::  po            ! specific humidity [kg/kg]
    1537       REAL, DIMENSION(klon,klev), INTENT(IN)      ::  zqta          ! specific humidity within thermals [kg/kg]
    1538       REAL, DIMENSION(klon,klev+1), INTENT(IN)    ::  fraca         ! Fraction of the mesh covered by thermals [0-1]
    1539       REAL, DIMENSION(klon,klev), INTENT(IN)      ::  zpspsk
    1540       REAL, DIMENSION(klon,klev+1), INTENT(IN)    ::  paprs         ! Pressure at layer interfaces [Pa]
    1541       REAL, DIMENSION(klon,klev), INTENT(IN)      ::  pplay         ! Pressure at the center of layers [Pa]
    1542       REAL, DIMENSION(klon,klev), INTENT(IN)      ::  ztla          ! Liquid temp [K]
    1543       REAL, DIMENSION(klon,klev), INTENT(INOUT)      ::  zthl       ! Liquid potential temp [K]
    1544       REAL, DIMENSION(klon,klev), INTENT(IN)      ::  ratqs         ! Parameter that determines the width of the total water distrib.
    1545       REAL, DIMENSION(klon),      INTENT(IN)      ::  zqs           ! Saturation specific humidity in the mesh [kg/kg]
    1546       REAL, DIMENSION(klon,klev+1), INTENT(IN)    ::  snowflux      ! snow flux at the interface of the layer [kg/m2/s]
    1547       INTEGER,                      INTENT(IN)    ::  iflag_mpc_bl  ! option flag for mpc boundary layer clouds param.
     1527      REAL, DIMENSION(klon),      INTENT(IN)      ::  temp          ! Temperature (liquid temperature) in the mesh [K] : has seen evap of precip
     1528      REAL, DIMENSION(klon),      INTENT(IN)      ::  qt            ! total water specific humidity in the mesh [kg/kg]: has seen evap of precip
     1529      REAL, DIMENSION(klon),      INTENT(IN)      ::  qt_th         ! total water specific humidity in thermals [kg/kg]: has not seen evap of precip
     1530      REAL, DIMENSION(klon),      INTENT(IN)      ::  thetal_th     ! Liquid potential temperature in thermals [K]: has not seen the evap of precip
     1531      REAL, DIMENSION(klon),      INTENT(IN)      ::  frac_th       ! Fraction of the mesh covered by thermals [0-1]
     1532      REAL, DIMENSION(klon),      INTENT(IN)      ::  zpspsk        ! Exner potential
     1533      REAL, DIMENSION(klon),      INTENT(IN)      ::  paprsup       ! Pressure at top layer interface [Pa]
     1534      REAL, DIMENSION(klon),      INTENT(IN)      ::  paprsdn       ! Pressure at bottom layer interface [Pa]
     1535      REAL, DIMENSION(klon),      INTENT(IN)      ::  play          ! Pressure of layers [Pa]
     1536      REAL, DIMENSION(klon),      INTENT(IN)      ::  ratqs         ! Parameter that determines the width of the total water distrib.
    15481537
    15491538
    15501539      INTEGER, DIMENSION(klon,klev), INTENT(INOUT) :: mpc_bl_points  ! grid points with convective (thermals) mixed phase clouds
    15511540     
    1552       REAL, DIMENSION(klon,klev), INTENT(OUT)      ::  ctot         ! Cloud fraction [0-1]
    1553       REAL, DIMENSION(klon,klev), INTENT(OUT)      ::  ctot_vol     ! Volume cloud fraction [0-1]
     1541      REAL, DIMENSION(klon),      INTENT(OUT)      ::  ctot         ! Cloud fraction [0-1]
     1542      REAL, DIMENSION(klon),      INTENT(OUT)      ::  ctot_vol     ! Volume cloud fraction [0-1]
    15541543      REAL, DIMENSION(klon),      INTENT(OUT)      ::  qcloud       ! In cloud total water content [kg/kg]
    1555       REAL, DIMENSION(klon),      INTENT(OUT)      ::  qincloud       ! In cloud condensed water content [kg/kg]
    1556       REAL, DIMENSION(klon,klev), INTENT(OUT)      ::  icefrac      ! Fraction of ice in clouds [0-1]
    1557       real, dimension(klon,klev), intent(out) :: cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv
     1544      REAL, DIMENSION(klon),      INTENT(OUT)      ::  qincloud     ! In cloud condensed water content [kg/kg]
     1545      REAL, DIMENSION(klon),      INTENT(OUT)      ::  icefrac      ! Fraction of ice in clouds [0-1]
     1546      REAL, DIMENSION(klon),      INTENT(OUT)      ::  cloudth_sth  ! mean saturation deficit in thermals
     1547      REAL, DIMENSION(klon),      INTENT(OUT)      ::  cloudth_senv ! mean saturation deficit in environment
     1548      REAL, DIMENSION(klon),      INTENT(OUT)      ::  cloudth_sigmath ! std of saturation deficit in thermals
     1549      REAL, DIMENSION(klon),      INTENT(OUT)      ::  cloudth_sigmaenv ! std of saturation deficit in environment
    15581550
    15591551
    … …  
    15621554      INTEGER itap,ind1,l,ig,iter,k
    15631555      INTEGER iflag_topthermals, niter
    1564       LOGICAL falseklon(klon)
    1565 
    1566       REAL zqsatth(klon), zqsatenv(klon), zqsatenvonly(klon)
    1567       REAL sigma1(klon,klev)                                                         
    1568       REAL sigma2(klon,klev)
    1569       REAL qcth(klon,klev)
    1570       REAL qcenv(klon,klev)
    1571       REAL qctot(klon,klev)
    1572       REAL cth(klon,klev)
    1573       REAL cenv(klon,klev)   
    1574       REAL cth_vol(klon,klev)
    1575       REAL cenv_vol(klon,klev)
    1576       REAL rneb(klon,klev)
    1577       REAL zqenv(klon), zqth(klon), zqenvonly(klon)
    1578       REAL qsatmmussig1,qsatmmussig2,sqrtpi,sqrt2,sqrt2pi,pi
    1579       REAL rdd,cppd,Lv
     1556
     1557      REAL qcth(klon)
     1558      REAL qcenv(klon)
     1559      REAL qctot(klon)
     1560      REAL cth(klon)
     1561      REAL cenv(klon)   
     1562      REAL cth_vol(klon)
     1563      REAL cenv_vol(klon)
     1564      REAL qt_env(klon), thetal_env(klon)
     1565      REAL icefracenv(klon), icefracth(klon)
     1566      REAL sqrtpi,sqrt2,sqrt2pi
    15801567      REAL alth,alenv,ath,aenv
    15811568      REAL sth,senv,sigma1s,sigma2s,sigma1s_fraca,sigma1s_ratqs
    … …  
    15851572      REAL IntJ,IntI1,IntI2,IntI3,IntJ_CF,IntI1_CF,IntI3_CF,coeffqlenv,coeffqlth
    15861573      REAL zdelta,qsatbef,zcor
    1587       REAL Tbefth(klon), Tbefenv(klon), Tbefenvonly(klon), rhoth(klon)
     1574      REAL Tbefth(klon), Tbefenv(klon), zrho(klon), rhoth(klon)
    15881575      REAL qlbef
    1589       REAL dqsatenv(klon), dqsatth(klon), dqsatenvonly(klon)
     1576      REAL dqsatenv(klon), dqsatth(klon)
    15901577      REAL erf
    15911578      REAL zpdf_sig(klon),zpdf_k(klon),zpdf_delta(klon)
    15921579      REAL zpdf_a(klon),zpdf_b(klon),zpdf_e1(klon),zpdf_e2(klon)
    1593       REAL rhodz(klon,klev)
    1594       REAL zrho(klon,klev)
    1595       REAL dz(klon,klev)
    1596       REAL qslenv(klon), qslth(klon)
     1580      REAL rhodz(klon)
     1581      REAL rho(klon)
     1582      REAL dz(klon)
     1583      REAL qslenv(klon), qslth(klon), qsienv(klon), qsith(klon)
    15971584      REAL alenvl, aenvl
    15981585      REAL sthi, sthl, sthil, althl, athl, althi, athi, sthlc, deltasthc, sigma2sc
    15991586      REAL senvi, senvl, qbase, sbase, qliqth, qiceth
    16001587      REAL qmax, ttarget, stmp, cout, coutref, fraci
    1601       REAL maxi, mini, pas, temp_lim
    1602       REAL deltazlev_mpc(klev), temp_mpc(klev), pres_mpc(klev), fraca_mpc(klev+1), snowf_mpc(klev+1)
     1588      REAL maxi, mini, pas
    16031589
    16041590      ! Modifty the saturation deficit PDF in thermals
    … …  
    16151601! Few initial checksS
    16161602
    1617       IF (iflag_cloudth_vert.NE.3) THEN
    1618          abort_message = 'clouth_mpc cannot be used if iflag_cloudth_vert .NE. 3'
    1619          CALL abort_physic(routname,abort_message,1)
    1620       ENDIF
    1621 
    16221603      DO k = 1,klev
    16231604      DO ind1 = 1, klon
    1624         rhodz(ind1,k) = (paprs(ind1,k)-paprs(ind1,k+1))/rg !kg/m2
    1625         zrho(ind1,k) = pplay(ind1,k)/temp(ind1,k)/rd !kg/m3
    1626         dz(ind1,k) = rhodz(ind1,k)/zrho(ind1,k) !m : epaisseur de la couche en metre
     1605        rhodz(ind1) = (paprsdn(ind1)-paprsup(ind1))/rg !kg/m2
     1606        zrho(ind1)  = play(ind1)/temp(ind1)/rd !kg/m3
     1607        dz(ind1)    = rhodz(ind1)/zrho(ind1) !m : epaisseur de la couche en metre
    16271608      END DO
    16281609      END DO
    16291610
    16301611
    1631       sigma1(:,:)=0.
    1632       sigma2(:,:)=0.
    1633       qcth(:,:)=0.
    1634       qcenv(:,:)=0. 
    1635       qctot(:,:)=0.
    1636       qlth(:,ind2)=0.
    1637       qith(:,ind2)=0.
    1638       wiceth(:,ind2)=0.
    1639       rneb(:,:)=0.
    1640       qcloud(:)=0.
    1641       cth(:,:)=0.
    1642       cenv(:,:)=0.
    1643       ctot(:,:)=0.
    1644       cth_vol(:,:)=0.
    1645       cenv_vol(:,:)=0.
    1646       ctot_vol(:,:)=0.
    1647       falseklon(:)=.false.
    1648       qsatmmussig1=0.
    1649       qsatmmussig2=0.
    1650       rdd=287.04
    1651       cppd=1005.7
    1652       pi=3.14159
    1653       sqrt2pi=sqrt(2.*pi)
     1612      icefracth(:)=0.
     1613      icefracenv(:)=0.
     1614      sqrt2pi=sqrt(2.*rpi)
    16541615      sqrt2=sqrt(2.)
    1655       sqrtpi=sqrt(pi)
    1656       icefrac(:,ind2)=0.
    1657       althl=0.
    1658       althi=0.
    1659       athl=0.
    1660       athi=0.
    1661       senvl=0.
    1662       senvi=0.
    1663       sthi=0.
    1664       sthl=0.
    1665       sthil=0.
     1616      sqrtpi=sqrt(rpi)
     1617      icefrac(:)=0.
    16661618
    16671619!-------------------------------------------------------------------------------
    … …  
    16691621!-------------------------------------------------------------------------------
    16701622
    1671       temp_lim=RTT-40.0
    16721623
    16731624      DO ind1=1,klon
    1674             IF ((temp(ind1,ind2) .LT. RTT) .AND. (temp(ind1,ind2) .GT. temp_lim) &
    1675             .AND. (iflag_mpc_bl .GE. 1) .AND. (ind2<=klev-2)  &
    1676             .AND. (ztv(ind1,1).GT.ztv(ind1,2)) .AND.(fraca(ind1,ind2).GT.1.e-10)) THEN
     1625            IF ((temp(ind1) .LT. RTT) .AND. (temp(ind1) .GT. temp_nowater) &
     1626            .AND. (ind2<=klev-2)  &
     1627            .AND. (frac_th(ind1).GT.min_frac_th_cld)) THEN
    16771628                mpc_bl_points(ind1,ind2)=1
    16781629            ELSE
    … …  
    16861637!------------------------------------------------------------------------------- 
    16871638
    1688 ! calculation of temperature, humidity and saturation specific humidity
    1689 
    1690 Tbefenv(:)=zthl(:,ind2)*zpspsk(:,ind2)
    1691 Tbefth(:)=ztla(:,ind2)*zpspsk(:,ind2)
    1692 Tbefenvonly(:)=temp(:,ind2)
    1693 rhoth(:)=paprs(:,ind2)/Tbefth(:)/RD
    1694 
    1695 zqenv(:)=(po(:)-fraca(:,ind2)*zqta(:,ind2))/(1.-fraca(:,ind2)) !qt = a*qtth + (1-a)*qtenv
    1696 zqth(:)=zqta(:,ind2)
    1697 zqenvonly(:)=po(:)
    1698 
    1699 
    1700 CALL CALC_QSAT_ECMWF(klon,Tbefenvonly,zqenvonly,paprs(:,ind2),RTT,0,.false.,zqsatenvonly,dqsatenv)
    1701 
    1702 CALL CALC_QSAT_ECMWF(klon,Tbefenv,zqenv,paprs(:,ind2),RTT,0,.false.,zqsatenv,dqsatenv)
    1703 CALL CALC_QSAT_ECMWF(klon,Tbefenv,zqenv,paprs(:,ind2),RTT,1,.false.,qslenv,dqsatenv)
    1704 
    1705 CALL CALC_QSAT_ECMWF(klon,Tbefth,zqth,paprs(:,ind2),RTT,1,.false.,qslth,dqsatth)
    1706 CALL CALC_QSAT_ECMWF(klon,Tbefth,zqth,paprs(:,ind2),RTT,2,.false.,qsith(:,ind2),dqsatth)
    1707 CALL CALC_QSAT_ECMWF(klon,Tbefth,zqth,paprs(:,ind2),RTT,0,.false.,zqsatth,dqsatth)
    1708 
    1709 
    1710   DO ind1=1,klon
    1711 
    1712 
    1713     IF ((ztv(ind1,1).GT.ztv(ind1,2)).AND.(fraca(ind1,ind2).GT.1.e-10)) THEN !Thermal and environnement
    1714 
     1639! initialisations and calculation of temperature, humidity and saturation specific humidity
     1640
     1641DO ind1=1,klon
     1642 
     1643   rhodz(ind1)   = (paprsdn(ind1)-paprsup(ind1))/rg  ! kg/m2
     1644   rho(ind1)     = play(ind1)/temp(ind1)/rd          ! kg/m3
     1645   dz(ind1)      = rhodz(ind1)/rho(ind1)             ! m : epaisseur de la couche en metre
     1646   Tbefenv(ind1) = temp(ind1)
     1647   thetal_env(ind1) = Tbefenv(ind1)/zpspsk(ind1)
     1648   Tbefth(ind1)  = thetal_th(ind1)*zpspsk(ind1)
     1649   rhoth(ind1)   = play(ind1)/Tbefth(ind1)/RD
     1650   qt_env(ind1)  = (qt(ind1)-frac_th(ind1)*qt_th(ind1))/(1.-frac_th(ind1)) !qt = a*qtth + (1-a)*qtenv
     1651
     1652ENDDO
     1653
     1654! Calculation of saturation specific humidity
     1655CALL CALC_QSAT_ECMWF(klon,Tbefenv,qt_env,play,RTT,1,.false.,qslenv,dqsatenv)
     1656CALL CALC_QSAT_ECMWF(klon,Tbefenv,qt_env,play,RTT,2,.false.,qsienv,dqsatenv)
     1657CALL CALC_QSAT_ECMWF(klon,Tbefth,qt_th,play,RTT,1,.false.,qslth,dqsatth)
     1658
     1659
     1660DO ind1=1,klon
     1661
     1662
     1663    IF (frac_th(ind1).GT.min_frac_th_cld) THEN !Thermal and environnement
     1664
     1665! unlike in the other cloudth routine,
     1666! We consider distributions of the saturation deficit WRT liquid
     1667! at positive AND negative celcius temperature
     1668! subsequently, cloud fraction corresponds to the part of the pdf corresponding
     1669! to superstauration with respect to liquid WHATEVER the temperature
    17151670
    17161671! Environment:
    17171672
    17181673
    1719         IF (Tbefenv(ind1) .GE. RTT) THEN
    1720                Lv=RLVTT
    1721         ELSE
    1722                Lv=RLSTT
    1723         ENDIF
    1724        
    1725 
    1726         alenv=(0.622*Lv*zqsatenv(ind1))/(rdd*zthl(ind1,ind2)**2)     !qsl, p84
    1727         aenv=1./(1.+(alenv*Lv/cppd))                                      !al, p84
    1728         senv=aenv*(po(ind1)-zqsatenv(ind1))                          !s, p84
     1674        alenv=(0.622*RLVTT*qslenv(ind1))/(rd*thetal_env(ind1)**2)     
     1675        aenv=1./(1.+(alenv*RLVTT/rcpd))                             
     1676        senv=aenv*(qt_env(ind1)-qslenv(ind1))                           
    17291677     
    1730         ! For MPCs:
    1731         IF (mpc_bl_points(ind1,ind2) .EQ. 1) THEN
    1732         alenvl=(0.622*RLVTT*qslenv(ind1))/(rdd*zthl(ind1,ind2)**2)     
    1733         aenvl=1./(1.+(alenv*Lv/cppd))         
    1734         senvl=aenvl*(po(ind1)-qslenv(ind1))   
    1735         senvi=senv
    1736         ENDIF
    1737 
    17381678
    17391679! Thermals:
    17401680
    17411681
    1742         IF (Tbefth(ind1) .GE. RTT) THEN
    1743             Lv=RLVTT
    1744         ELSE
    1745             Lv=RLSTT
    1746         ENDIF
    1747 
    1748        
    1749         alth=(0.622*Lv*zqsatth(ind1))/(rdd*ztla(ind1,ind2)**2)       
    1750         ath=1./(1.+(alth*Lv/cppd))                                                         
    1751         sth=ath*(zqta(ind1,ind2)-zqsatth(ind1))                     
    1752 
    1753        ! For MPCs:
    1754         IF (mpc_bl_points(ind1,ind2) .GT. 0) THEN
    1755          althi=alth
    1756          althl=(0.622*RLVTT*qslth(ind1))/(rdd*ztla(ind1,ind2)**2)                   
    1757          athl=1./(1.+(alth*RLVTT/cppd))
    1758          athi=alth     
    1759          sthl=athl*(zqta(ind1,ind2)-qslth(ind1))   
    1760          sthi=athi*(zqta(ind1,ind2)-qsith(ind1,ind2)) 
    1761          sthil=athi*(zqta(ind1,ind2)-qslth(ind1))
    1762         ENDIF     
    1763 
    1764 
    1765 
    1766 !-------------------------------------------------------------------------------
    1767 !  Version 3: Changes by J. Jouhaud; condensation for q > -delta s
    1768 !  Rq: in this subroutine, we assume iflag_clouth_vert .EQ. 3
    1769 !-------------------------------------------------------------------------------
    1770 
    1771         IF (mpc_bl_points(ind1,ind2) .EQ. 0) THEN ! No BL MPC
    1772 
    1773        ! Standard deviation of the distributions
    1774 
    1775            sigma1s_fraca = (sigma1s_factor**0.5)*(fraca(ind1,ind2)**sigma1s_power) / &
    1776            &                (1-fraca(ind1,ind2))*((sth-senv)**2)**0.5
     1682        alth=(0.622*RLVTT*qslth(ind1))/(rd*thetal_th(ind1)**2)       
     1683        ath=1./(1.+(alth*RLVTT/rcpd))                                                         
     1684        sth=ath*(qt_th(ind1)-qslth(ind1))                     
     1685
     1686
     1687!       IF (mpc_bl_points(ind1,ind2) .EQ. 0) THEN ! No BL MPC
     1688
     1689
     1690! Standard deviation of the distributions
     1691
     1692           sigma1s_fraca = (sigma1s_factor**0.5)*(frac_th(ind1)**sigma1s_power) / &
     1693           &                (1-frac_th(ind1))*((sth-senv)**2)**0.5
    17771694
    17781695           IF (cloudth_ratqsmin>0.) THEN
    1779              sigma1s_ratqs = cloudth_ratqsmin*po(ind1)
     1696             sigma1s_ratqs = cloudth_ratqsmin*qt(ind1)
    17801697           ELSE
    1781              sigma1s_ratqs = ratqs(ind1,ind2)*po(ind1)
     1698             sigma1s_ratqs = ratqs(ind1)*qt(ind1)
    17821699           ENDIF
    17831700 
    17841701           sigma1s = sigma1s_fraca + sigma1s_ratqs
    17851702           IF (iflag_ratqs.eq.11) then
    1786               sigma1s = ratqs(ind1,ind2)*po(ind1)*aenv
     1703              sigma1s = ratqs(ind1)*qt(ind1)*aenv
    17871704           ENDIF
    1788            sigma2s=(sigma2s_factor*(((sth-senv)**2)**0.5)/((fraca(ind1,ind2)+0.02)**sigma2s_power))+0.002*zqta(ind1,ind2)
     1705           sigma2s=(sigma2s_factor*(((sth-senv)**2)**0.5)/((frac_th(ind1)+0.02)**sigma2s_power))+0.002*qt_th(ind1)
    17891706
    17901707
    … …  
    18011718           exp_xth2 = exp(-1.*xth2**2)
    18021719     
    1803       !surface CF
    1804 
    1805            cth(ind1,ind2)=0.5*(1.-1.*erf(xth1))
    1806            cenv(ind1,ind2)=0.5*(1.-1.*erf(xenv1))
    1807            ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)
    1808 
    1809 
    1810       !volume CF and condensed water
     1720!surface CF
     1721
     1722           cth(ind1)=0.5*(1.-1.*erf(xth1))
     1723           cenv(ind1)=0.5*(1.-1.*erf(xenv1))
     1724           ctot(ind1)=frac_th(ind1)*cth(ind1)+(1.-1.*frac_th(ind1))*cenv(ind1)
     1725
     1726
     1727!volume CF, condensed water and ice fraction
    18111728
    18121729            !environnement
    … …  
    18161733
    18171734            IF (deltasenv .LT. 1.e-10) THEN
    1818               qcenv(ind1,ind2)=IntJ
    1819               cenv_vol(ind1,ind2)=IntJ_CF
     1735              qcenv(ind1)=IntJ
     1736              cenv_vol(ind1)=IntJ_CF
    18201737            ELSE
    18211738              IntI1=(((senv+deltasenv)**2+(sigma1s)**2)/(8*deltasenv))*(erf(xenv2)-erf(xenv1))
    … …  
    18241741              IntI1_CF=((senv+deltasenv)*(erf(xenv2)-erf(xenv1)))/(4*deltasenv)
    18251742              IntI3_CF=(sqrt2*sigma1s*(exp_xenv1-exp_xenv2))/(4*sqrtpi*deltasenv)
    1826               qcenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
    1827               cenv_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF
     1743              qcenv(ind1)=IntJ+IntI1+IntI2+IntI3
     1744              cenv_vol(ind1)=IntJ_CF+IntI1_CF+IntI3_CF
     1745              IF (Tbefenv(ind1) .LT. temp_nowater) THEN
     1746                  ! freeze all droplets in cirrus temperature regime
     1747                  icefracenv(ind1)=1.
     1748              ENDIF
    18281749            ENDIF
    18291750             
    … …  
    18361757     
    18371758            IF (deltasth .LT. 1.e-10) THEN
    1838               qcth(ind1,ind2)=IntJ
    1839               cth_vol(ind1,ind2)=IntJ_CF
     1759              qcth(ind1)=IntJ
     1760              cth_vol(ind1)=IntJ_CF
    18401761            ELSE
    18411762              IntI1=(((sth+deltasth)**2+(sigma2s)**2)/(8*deltasth))*(erf(xth2)-erf(xth1))
    … …  
    18441765              IntI1_CF=((sth+deltasth)*(erf(xth2)-erf(xth1)))/(4*deltasth)
    18451766              IntI3_CF=(sqrt2*sigma2s*(exp_xth1-exp_xth2))/(4*sqrtpi*deltasth)
    1846               qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
    1847               cth_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF
     1767              qcth(ind1)=IntJ+IntI1+IntI2+IntI3
     1768              cth_vol(ind1)=IntJ_CF+IntI1_CF+IntI3_CF
     1769              IF (Tbefth(ind1) .LT. temp_nowater) THEN
     1770                  ! freeze all droplets in cirrus temperature regime
     1771                  icefracth(ind1)=1.
     1772              ENDIF
     1773
    18481774            ENDIF
    18491775
    1850               qctot(ind1,ind2)=fraca(ind1,ind2)*qcth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qcenv(ind1,ind2)
    1851               ctot_vol(ind1,ind2)=fraca(ind1,ind2)*cth_vol(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv_vol(ind1,ind2)
    1852              
    1853 
    1854             IF (cenv(ind1,ind2).LT.1.e-10.or.cth(ind1,ind2).LT.1.e-10) THEN
    1855                 ctot(ind1,ind2)=0.
    1856                 ctot_vol(ind1,ind2)=0.
    1857                 qcloud(ind1)=zqsatenv(ind1)
     1776            qctot(ind1)=frac_th(ind1)*qcth(ind1)+(1.-1.*frac_th(ind1))*qcenv(ind1)
     1777            ctot_vol(ind1)=frac_th(ind1)*cth_vol(ind1)+(1.-1.*frac_th(ind1))*cenv_vol(ind1)
     1778
     1779            IF (cenv(ind1).LT.min_neb_th.and.cth(ind1).LT.min_neb_th) THEN
     1780                ctot(ind1)=0.
     1781                ctot_vol(ind1)=0.
     1782                qcloud(ind1)=qslenv(ind1)
    18581783                qincloud(ind1)=0.
     1784                icefrac(ind1)=0.
    18591785            ELSE               
    1860                 qcloud(ind1)=qctot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
    1861                 qincloud(ind1)=qctot(ind1,ind2)/ctot(ind1,ind2)
     1786                qcloud(ind1)=qctot(ind1)/ctot(ind1)+qslenv(ind1)
     1787                qincloud(ind1)=qctot(ind1)/ctot(ind1)
     1788                IF (qctot(ind1) .GT. 0) THEN
     1789                  icefrac(ind1)=(frac_th(ind1)*qcth(ind1)*icefracth(ind1)+(1.-1.*frac_th(ind1))*qcenv(ind1)*icefracenv(ind1))/qctot(ind1)
     1790                  icefrac(ind1)=max(min(1.,icefrac(ind1)), 0.)
     1791                ENDIF
    18621792            ENDIF
    18631793
    18641794
    1865         ELSE ! mpc_bl_points>0
    1866 
    1867             ! Treat boundary layer mixed phase clouds
    1868            
    1869             ! thermals
    1870             !=========
    1871 
    1872             ! ice phase
    1873             !...........
    1874 
    1875             qiceth=0.
    1876             deltazlev_mpc=dz(ind1,:)
    1877             temp_mpc=ztla(ind1,:)*zpspsk(ind1,:)
    1878             pres_mpc=pplay(ind1,:)
    1879             fraca_mpc=fraca(ind1,:)
    1880             snowf_mpc=snowflux(ind1,:)
    1881             iflag_topthermals=0
    1882             IF ((mpc_bl_points(ind1,ind2) .EQ. 1) .AND. (mpc_bl_points(ind1,ind2+1) .EQ. 0))  THEN
    1883                 iflag_topthermals = 1
    1884             ELSE IF ((mpc_bl_points(ind1,ind2) .EQ. 1) .AND. (mpc_bl_points(ind1,ind2+1) .EQ. 1) &
    1885                     .AND. (mpc_bl_points(ind1,ind2+2) .EQ. 0) ) THEN
    1886                 iflag_topthermals = 2
    1887             ELSE
    1888                 iflag_topthermals = 0
    1889             ENDIF
    1890 
    1891             CALL ICE_MPC_BL_CLOUDS(ind1,ind2,klev,Ni,Ei,C_cap,d_top,iflag_topthermals,temp_mpc,pres_mpc,zqta(ind1,:), &
    1892                                    qsith(ind1,:),qlth(ind1,:),deltazlev_mpc,wiceth(ind1,:),fraca_mpc,qith(ind1,:))
    1893 
    1894             ! qmax calculation
    1895             sigma2s=(sigma2s_factor*((MAX((sthl-senvl),0.)**2)**0.5)/((fraca(ind1,ind2)+0.02)**sigma2s_power))+0.002*zqta(ind1,ind2)
    1896             deltasth=athl*vert_alpha_th*sigma2s
    1897             xth1=-(sthl+deltasth)/(sqrt(2.)*sigma2s)
    1898             xth2=-(sthl-deltasth)/(sqrt(2.)*sigma2s)
    1899             exp_xth1 = exp(-1.*xth1**2)
    1900             exp_xth2 = exp(-1.*xth2**2)
    1901             IntJ=0.5*sthl*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2
    1902             IntJ_CF=0.5*(1.-1.*erf(xth2))
    1903             IntI1=(((sthl+deltasth)**2+(sigma2s)**2)/(8*deltasth))*(erf(xth2)-erf(xth1))
    1904             IntI2=(sigma2s**2/(4*deltasth*sqrtpi))*(xth1*exp_xth1-xth2*exp_xth2)
    1905             IntI3=((sqrt2*sigma2s*(sthl+deltasth))/(4*sqrtpi*deltasth))*(exp_xth1-exp_xth2)
    1906             IntI1_CF=((sthl+deltasth)*(erf(xth2)-erf(xth1)))/(4*deltasth)
    1907             IntI3_CF=(sqrt2*sigma2s*(exp_xth1-exp_xth2))/(4*sqrtpi*deltasth)
    1908             qmax=MAX(IntJ+IntI1+IntI2+IntI3,0.)
    1909 
    1910 
    1911             ! Liquid phase
    1912             !................
    1913             ! We account for the effect of ice crystals in thermals on sthl
    1914             ! and on the width of the distribution
    1915 
    1916             sthlc=sthl*1./(1.+C_mpc*qith(ind1,ind2))  &
    1917                 + (1.-1./(1.+C_mpc*qith(ind1,ind2))) * athl*(qsith(ind1,ind2)-qslth(ind1)) 
    1918 
    1919             sigma2sc=(sigma2s_factor*((MAX((sthlc-senvl),0.)**2)**0.5) &
    1920                  /((fraca(ind1,ind2)+0.02)**sigma2s_power)) &
    1921                  +0.002*zqta(ind1,ind2)
    1922             deltasthc=athl*vert_alpha_th*sigma2sc
    1923      
    1924            
    1925             xth1=-(sthlc+deltasthc)/(sqrt(2.)*sigma2sc)
    1926             xth2=-(sthlc-deltasthc)/(sqrt(2.)*sigma2sc)           
    1927             exp_xth1 = exp(-1.*xth1**2)
    1928             exp_xth2 = exp(-1.*xth2**2)
    1929             IntJ=0.5*sthlc*(1-erf(xth2))+(sigma2sc/sqrt2pi)*exp_xth2
    1930             IntJ_CF=0.5*(1.-1.*erf(xth2))
    1931             IntI1=(((sthlc+deltasthc)**2+(sigma2sc)**2)/(8*deltasthc))*(erf(xth2)-erf(xth1))
    1932             IntI2=(sigma2sc**2/(4*deltasthc*sqrtpi))*(xth1*exp_xth1-xth2*exp_xth2)
    1933             IntI3=((sqrt2*sigma2sc*(sthlc+deltasthc))/(4*sqrtpi*deltasthc))*(exp_xth1-exp_xth2)
    1934             IntI1_CF=((sthlc+deltasthc)*(erf(xth2)-erf(xth1)))/(4*deltasthc)
    1935             IntI3_CF=(sqrt2*sigma2sc*(exp_xth1-exp_xth2))/(4*sqrtpi*deltasthc)
    1936             qliqth=IntJ+IntI1+IntI2+IntI3
    1937            
    1938             qlth(ind1,ind2)=MAX(0.,qliqth)
    1939 
    1940             ! Condensed water
    1941            
    1942             qcth(ind1,ind2)=qlth(ind1,ind2)+qith(ind1,ind2)
    1943 
    1944 
    1945             ! consistency with subgrid distribution
    1946            
    1947              IF ((qcth(ind1,ind2) .GT. qmax) .AND. (qcth(ind1,ind2) .GT. 0)) THEN
    1948                  fraci=qith(ind1,ind2)/qcth(ind1,ind2)
    1949                  qcth(ind1,ind2)=qmax
    1950                  qith(ind1,ind2)=fraci*qmax
    1951                  qlth(ind1,ind2)=(1.-fraci)*qmax
    1952              ENDIF
    1953 
    1954             ! Cloud Fraction
    1955             !...............
    1956             ! calculation of qbase which is the value of the water vapor within mixed phase clouds
    1957             ! such that the total water in cloud = qbase+qliqth+qiceth
    1958             ! sbase is the value of s such that int_sbase^\intfy s ds = cloud fraction
    1959             ! sbase and qbase calculation (note that sbase is wrt liq so negative)
    1960             ! look for an approximate solution with iteration
    1961            
    1962             ttarget=qcth(ind1,ind2)
    1963             mini= athl*(qsith(ind1,ind2)-qslth(ind1))
    1964             maxi= 0. !athl*(3.*sqrt(sigma2s))
    1965             niter=20
    1966             pas=(maxi-mini)/niter
    1967             stmp=mini
    1968             sbase=stmp
    1969             coutref=1.E6
    1970             DO iter=1,niter
    1971                 cout=ABS(sigma2s/SQRT(2.*RPI)*EXP(-((sthl-stmp)/sigma2s)**2)+(sthl-stmp)/SQRT(2.)*(1.-erf(-(sthl-stmp)/sigma2s)) &
    1972                      + stmp/2.*(1.-erf(-(sthl-stmp)/sigma2s)) -ttarget)
    1973                IF (cout .LT. coutref) THEN
    1974                      sbase=stmp
    1975                      coutref=cout
    1976                 ELSE
    1977                      stmp=stmp+pas
    1978                 ENDIF
    1979             ENDDO
    1980             qbase=MAX(0., sbase/athl+qslth(ind1))
    1981 
    1982             ! surface cloud fraction in thermals
    1983             cth(ind1,ind2)=0.5*(1.-erf((sbase-sthl)/sqrt(2.)/sigma2s))
    1984             cth(ind1,ind2)=MIN(MAX(cth(ind1,ind2),0.),1.)
    1985 
    1986 
    1987             !volume cloud fraction in thermals
    1988             !to be checked
    1989             xth1=-(sthl+deltasth-sbase)/(sqrt(2.)*sigma2s)
    1990             xth2=-(sthl-deltasth-sbase)/(sqrt(2.)*sigma2s)           
    1991             exp_xth1 = exp(-1.*xth1**2)
    1992             exp_xth2 = exp(-1.*xth2**2)
    1993 
    1994             IntJ=0.5*sthl*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2
    1995             IntJ_CF=0.5*(1.-1.*erf(xth2))
    1996      
    1997             IF (deltasth .LT. 1.e-10) THEN
    1998               cth_vol(ind1,ind2)=IntJ_CF
    1999             ELSE
    2000               IntI1=(((sthl+deltasth-sbase)**2+(sigma2s)**2)/(8*deltasth))*(erf(xth2)-erf(xth1))
    2001               IntI2=(sigma2s**2/(4*deltasth*sqrtpi))*(xth1*exp_xth1-xth2*exp_xth2)
    2002               IntI3=((sqrt2*sigma2s*(sthl+deltasth))/(4*sqrtpi*deltasth))*(exp_xth1-exp_xth2)
    2003               IntI1_CF=((sthl-sbase+deltasth)*(erf(xth2)-erf(xth1)))/(4*deltasth)
    2004               IntI3_CF=(sqrt2*sigma2s*(exp_xth1-exp_xth2))/(4*sqrtpi*deltasth)
    2005               cth_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF
    2006             ENDIF
    2007               cth_vol(ind1,ind2)=MIN(MAX(0.,cth_vol(ind1,ind2)),1.)
    2008 
    2009 
    2010 
    2011             ! Environment
    2012             !=============
    2013             ! In the environment/downdrafts, only liquid clouds
    2014             ! See Shupe et al. 2008, JAS
    2015 
    2016             ! standard deviation of the distribution in the environment
    2017             sth=sthl
    2018             senv=senvl
    2019             sigma1s_fraca = (sigma1s_factor**0.5)*(fraca(ind1,ind2)**sigma1s_power) / &
    2020                           &                (1-fraca(ind1,ind2))*(MAX((sth-senv),0.)**2)**0.5
    2021             ! for mixed phase clouds, there is no contribution from large scale ratqs to the distribution
    2022             ! in the environement
    2023 
    2024             sigma1s_ratqs=1E-10
    2025             IF (cloudth_ratqsmin>0.) THEN
    2026                 sigma1s_ratqs = cloudth_ratqsmin*po(ind1)
    2027             ENDIF
    2028 
    2029             sigma1s = sigma1s_fraca + sigma1s_ratqs
    2030             IF (iflag_ratqs.eq.11) then
    2031                sigma1s = ratqs(ind1,ind2)*po(ind1)*aenv
    2032             ENDIF
    2033             IF (iflag_ratqs.eq.11) then
    2034               sigma1s = ratqs(ind1,ind2)*po(ind1)*aenvl
    2035             ENDIF
    2036             deltasenv=aenvl*vert_alpha*sigma1s
    2037             xenv1=-(senvl+deltasenv)/(sqrt(2.)*sigma1s)
    2038             xenv2=-(senvl-deltasenv)/(sqrt(2.)*sigma1s)
    2039             exp_xenv1 = exp(-1.*xenv1**2)
    2040             exp_xenv2 = exp(-1.*xenv2**2)
    2041 
    2042             !surface CF
    2043             cenv(ind1,ind2)=0.5*(1.-1.*erf(xenv1))
    2044 
    2045             !volume CF and condensed water
    2046             IntJ=0.5*senvl*(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp_xenv2
    2047             IntJ_CF=0.5*(1.-1.*erf(xenv2))
    2048 
    2049             IF (deltasenv .LT. 1.e-10) THEN
    2050               qcenv(ind1,ind2)=IntJ
    2051               cenv_vol(ind1,ind2)=IntJ_CF
    2052             ELSE
    2053               IntI1=(((senvl+deltasenv)**2+(sigma1s)**2)/(8*deltasenv))*(erf(xenv2)-erf(xenv1))
    2054               IntI2=(sigma1s**2/(4*deltasenv*sqrtpi))*(xenv1*exp_xenv1-xenv2*exp_xenv2)
    2055               IntI3=((sqrt2*sigma1s*(senv+deltasenv))/(4*sqrtpi*deltasenv))*(exp_xenv1-exp_xenv2)
    2056               IntI1_CF=((senvl+deltasenv)*(erf(xenv2)-erf(xenv1)))/(4*deltasenv)
    2057               IntI3_CF=(sqrt2*sigma1s*(exp_xenv1-exp_xenv2))/(4*sqrtpi*deltasenv)
    2058               qcenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3 ! only liquid water in environment
    2059               cenv_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF
    2060             ENDIF
    2061 
    2062             qcenv(ind1,ind2)=MAX(qcenv(ind1,ind2),0.)
    2063             cenv_vol(ind1,ind2)=MIN(MAX(cenv_vol(ind1,ind2),0.),1.)
    2064 
    2065 
    2066            
    2067             ! Thermals + environment
    2068             ! =======================
    2069             ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)
    2070             qctot(ind1,ind2)=fraca(ind1,ind2)*qcth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qcenv(ind1,ind2)
    2071             ctot_vol(ind1,ind2)=fraca(ind1,ind2)*cth_vol(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv_vol(ind1,ind2)
    2072             IF (qcth(ind1,ind2) .GT. 0) THEN
    2073                icefrac(ind1,ind2)=fraca(ind1,ind2)*qith(ind1,ind2) &
    2074                     /(fraca(ind1,ind2)*qcth(ind1,ind2) &
    2075                     +(1.-1.*fraca(ind1,ind2))*qcenv(ind1,ind2))
    2076                 icefrac(ind1,ind2)=MAX(MIN(1.,icefrac(ind1,ind2)),0.)
    2077             ELSE
    2078                 icefrac(ind1,ind2)=0.
    2079             ENDIF
    2080 
    2081             IF (cenv(ind1,ind2).LT.1.e-10.or.cth(ind1,ind2).LT.1.e-10) THEN
    2082                 ctot(ind1,ind2)=0.
    2083                 ctot_vol(ind1,ind2)=0.
    2084                 qincloud(ind1)=0.
    2085                 qcloud(ind1)=zqsatenv(ind1)
    2086             ELSE               
    2087                 qcloud(ind1)=fraca(ind1,ind2)*(qcth(ind1,ind2)/cth(ind1,ind2)+qbase) &
    2088                             +(1.-1.*fraca(ind1,ind2))*(qcenv(ind1,ind2)/cenv(ind1,ind2)+qslenv(ind1))
    2089                 qincloud(ind1)=MAX(fraca(ind1,ind2)*(qcth(ind1,ind2)/cth(ind1,ind2)) &
    2090                             +(1.-1.*fraca(ind1,ind2))*(qcenv(ind1,ind2)/cenv(ind1,ind2)),0.)
    2091             ENDIF
    2092 
    2093         ENDIF ! mpc_bl_points
    2094 
     1795!        ELSE ! mpc_bl_points>0
    20951796
    20961797    ELSE  ! gaussian for environment only
    20971798
    20981799     
    2099 
    2100 
    2101         IF (Tbefenvonly(ind1) .GE. RTT) THEN
    2102                 Lv=RLVTT
    2103         ELSE
    2104                 Lv=RLSTT
    2105         ENDIF
    2106        
    2107 
    2108         zthl(ind1,ind2)=temp(ind1,ind2)*(101325./paprs(ind1,ind2))**(rdd/cppd)
    2109         alenv=(0.622*Lv*zqsatenvonly(ind1))/(rdd*zthl(ind1,ind2)**2)
    2110         aenv=1./(1.+(alenv*Lv/cppd))
    2111         senv=aenv*(po(ind1)-zqsatenvonly(ind1))
     1800        alenv=(0.622*RLVTT*qslenv(ind1))/(rd*thetal_env(ind1)**2)
     1801        aenv=1./(1.+(alenv*RLVTT/rcpd))
     1802        senv=aenv*(qt_env(ind1)-qslenv(ind1))
    21121803        sth=0.
    2113      
    2114         sigma1s=ratqs(ind1,ind2)*zqenvonly(ind1)
     1804        sigma1s=ratqs(ind1)*qt_env(ind1)
    21151805        sigma2s=0.
    21161806
    2117         sqrt2pi=sqrt(2.*pi)
    21181807        xenv=senv/(sqrt(2.)*sigma1s)
    2119         ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv))
    2120         ctot_vol(ind1,ind2)=ctot(ind1,ind2)
    2121         qctot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))
    2122      
    2123         IF (ctot(ind1,ind2).LT.1.e-3) THEN
    2124           ctot(ind1,ind2)=0.
    2125           qcloud(ind1)=zqsatenvonly(ind1)
     1808        cenv(ind1)=0.5*(1.-1.*erf(xenv))
     1809        ctot(ind1)=0.5*(1.+1.*erf(xenv))
     1810        ctot_vol(ind1)=ctot(ind1)
     1811        qctot(ind1)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1))
     1812
     1813
     1814        IF (ctot(ind1).LT.min_neb_th) THEN
     1815          ctot(ind1)=0.
     1816          qcloud(ind1)=qslenv(ind1)
    21261817          qincloud(ind1)=0.
    21271818        ELSE
    2128           qcloud(ind1)=qctot(ind1,ind2)/ctot(ind1,ind2)+zqsatenvonly(ind1)
    2129           qincloud(ind1)=MAX(qctot(ind1,ind2)/ctot(ind1,ind2),0.)
     1819          qcloud(ind1)=qctot(ind1)/ctot(ind1)+qslenv(ind1)
     1820          qincloud(ind1)=MAX(qctot(ind1)/ctot(ind1),0.)
    21301821        ENDIF
    21311822 
    … …  
    21341825
    21351826    ! Outputs used to check the PDFs
    2136     cloudth_senv(ind1,ind2) = senv
    2137     cloudth_sth(ind1,ind2) = sth
    2138     cloudth_sigmaenv(ind1,ind2) = sigma1s
    2139     cloudth_sigmath(ind1,ind2) = sigma2s
     1827    cloudth_senv(ind1) = senv
     1828    cloudth_sth(ind1) = sth
     1829    cloudth_sigmaenv(ind1) = sigma1s
     1830    cloudth_sigmath(ind1) = sigma2s
    21401831
    21411832
    21421833    ENDDO       !loop on klon
    21431834
    2144     ! Calcule ice fall velocity in thermals
    2145 
    2146     CALL FALLICE_VELOCITY(klon,qith(:,ind2),Tbefth(:),rhoth(:),paprs(:,ind2),falseklon(:),wiceth(:,ind2))
    21471835
    21481836RETURN
    … …  
    21501838
    21511839END SUBROUTINE cloudth_mpc
     1840
     1841
     1842
     1843!        ELSE ! mpc_bl_points>0
     1844!
     1845!            ! Treat boundary layer mixed phase clouds
     1846!           
     1847!            ! thermals
     1848!            !=========
     1849!
     1850!           ! ice phase
     1851!           !...........
     1852!
     1853!            qiceth=0.
     1854!            deltazlev_mpc=dz(ind1,:)
     1855!            temp_mpc=ztla(ind1,:)*zpspsk(ind1,:)
     1856!            pres_mpc=pplay(ind1,:)
     1857!            fraca_mpc=fraca(ind1,:)
     1858!            snowf_mpc=snowflux(ind1,:)
     1859!            iflag_topthermals=0
     1860!            IF ((mpc_bl_points(ind1,ind2) .EQ. 1) .AND. (mpc_bl_points(ind1,ind2+1) .EQ. 0))  THEN
     1861!                iflag_topthermals = 1
     1862!            ELSE IF ((mpc_bl_points(ind1,ind2) .EQ. 1) .AND. (mpc_bl_points(ind1,ind2+1) .EQ. 1) &
     1863!                    .AND. (mpc_bl_points(ind1,ind2+2) .EQ. 0) ) THEN
     1864!                iflag_topthermals = 2
     1865!            ELSE
     1866!                iflag_topthermals = 0
     1867!            ENDIF
     1868!
     1869!            CALL ICE_MPC_BL_CLOUDS(ind1,ind2,klev,Ni,Ei,C_cap,d_top,iflag_topthermals,temp_mpc,pres_mpc,zqta(ind1,:), &
     1870!                                   qsith(ind1,:),qlth(ind1,:),deltazlev_mpc,wiceth(ind1,:),fraca_mpc,qith(ind1,:))
     1871!
     1872!            ! qmax calculation
     1873!            sigma2s=(sigma2s_factor*((MAX((sthl-senvl),0.)**2)**0.5)/((fraca(ind1,ind2)+0.02)**sigma2s_power))+0.002*zqta(ind1,ind2)
     1874!            deltasth=athl*vert_alpha_th*sigma2s
     1875!            xth1=-(sthl+deltasth)/(sqrt(2.)*sigma2s)
     1876!            xth2=-(sthl-deltasth)/(sqrt(2.)*sigma2s)
     1877!            exp_xth1 = exp(-1.*xth1**2)
     1878!            exp_xth2 = exp(-1.*xth2**2)
     1879!            IntJ=0.5*sthl*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2
     1880!            IntJ_CF=0.5*(1.-1.*erf(xth2))
     1881!            IntI1=(((sthl+deltasth)**2+(sigma2s)**2)/(8*deltasth))*(erf(xth2)-erf(xth1))
     1882!            IntI2=(sigma2s**2/(4*deltasth*sqrtpi))*(xth1*exp_xth1-xth2*exp_xth2)
     1883!            IntI3=((sqrt2*sigma2s*(sthl+deltasth))/(4*sqrtpi*deltasth))*(exp_xth1-exp_xth2)
     1884!            IntI1_CF=((sthl+deltasth)*(erf(xth2)-erf(xth1)))/(4*deltasth)
     1885!            IntI3_CF=(sqrt2*sigma2s*(exp_xth1-exp_xth2))/(4*sqrtpi*deltasth)
     1886!            qmax=MAX(IntJ+IntI1+IntI2+IntI3,0.)
     1887!
     1888!
     1889!            ! Liquid phase
     1890!            !................
     1891!            ! We account for the effect of ice crystals in thermals on sthl
     1892!            ! and on the width of the distribution
     1893!
     1894!            sthlc=sthl*1./(1.+C_mpc*qith(ind1,ind2))  &
     1895!                + (1.-1./(1.+C_mpc*qith(ind1,ind2))) * athl*(qsith(ind1,ind2)-qslth(ind1)) 
     1896!
     1897!            sigma2sc=(sigma2s_factor*((MAX((sthlc-senvl),0.)**2)**0.5) &
     1898!                 /((fraca(ind1,ind2)+0.02)**sigma2s_power)) &
     1899!                 +0.002*zqta(ind1,ind2)
     1900!            deltasthc=athl*vert_alpha_th*sigma2sc
     1901!     
     1902!           
     1903!            xth1=-(sthlc+deltasthc)/(sqrt(2.)*sigma2sc)
     1904!            xth2=-(sthlc-deltasthc)/(sqrt(2.)*sigma2sc)           
     1905!            exp_xth1 = exp(-1.*xth1**2)
     1906!            exp_xth2 = exp(-1.*xth2**2)
     1907!            IntJ=0.5*sthlc*(1-erf(xth2))+(sigma2sc/sqrt2pi)*exp_xth2
     1908!            IntJ_CF=0.5*(1.-1.*erf(xth2))
     1909!            IntI1=(((sthlc+deltasthc)**2+(sigma2sc)**2)/(8*deltasthc))*(erf(xth2)-erf(xth1))
     1910!            IntI2=(sigma2sc**2/(4*deltasthc*sqrtpi))*(xth1*exp_xth1-xth2*exp_xth2)
     1911!            IntI3=((sqrt2*sigma2sc*(sthlc+deltasthc))/(4*sqrtpi*deltasthc))*(exp_xth1-exp_xth2)
     1912!            IntI1_CF=((sthlc+deltasthc)*(erf(xth2)-erf(xth1)))/(4*deltasthc)
     1913!            IntI3_CF=(sqrt2*sigma2sc*(exp_xth1-exp_xth2))/(4*sqrtpi*deltasthc)
     1914!            qliqth=IntJ+IntI1+IntI2+IntI3
     1915!           
     1916!            qlth(ind1,ind2)=MAX(0.,qliqth)
     1917!
     1918!            ! Condensed water
     1919!           
     1920!            qcth(ind1,ind2)=qlth(ind1,ind2)+qith(ind1,ind2)
     1921!
     1922!
     1923!            ! consistency with subgrid distribution
     1924!           
     1925!             IF ((qcth(ind1,ind2) .GT. qmax) .AND. (qcth(ind1,ind2) .GT. 0)) THEN
     1926!                 fraci=qith(ind1,ind2)/qcth(ind1,ind2)
     1927!                 qcth(ind1,ind2)=qmax
     1928!                 qith(ind1,ind2)=fraci*qmax
     1929!                 qlth(ind1,ind2)=(1.-fraci)*qmax
     1930!             ENDIF
     1931!
     1932!            ! Cloud Fraction
     1933!            !...............
     1934!            ! calculation of qbase which is the value of the water vapor within mixed phase clouds
     1935!            ! such that the total water in cloud = qbase+qliqth+qiceth
     1936!            ! sbase is the value of s such that int_sbase^\intfy s ds = cloud fraction
     1937!            ! sbase and qbase calculation (note that sbase is wrt liq so negative)
     1938!            ! look for an approximate solution with iteration
     1939!           
     1940!            ttarget=qcth(ind1,ind2)
     1941!            mini= athl*(qsith(ind1,ind2)-qslth(ind1))
     1942!            maxi= 0. !athl*(3.*sqrt(sigma2s))
     1943!            niter=20
     1944!            pas=(maxi-mini)/niter
     1945!            stmp=mini
     1946!            sbase=stmp
     1947!            coutref=1.E6
     1948!            DO iter=1,niter
     1949!                cout=ABS(sigma2s/SQRT(2.*RPI)*EXP(-((sthl-stmp)/sigma2s)**2)+(sthl-stmp)/SQRT(2.)*(1.-erf(-(sthl-stmp)/sigma2s)) &
     1950!                     + stmp/2.*(1.-erf(-(sthl-stmp)/sigma2s)) -ttarget)
     1951!               IF (cout .LT. coutref) THEN
     1952!                     sbase=stmp
     1953!                     coutref=cout
     1954!                ELSE
     1955!                     stmp=stmp+pas
     1956!                ENDIF
     1957!            ENDDO
     1958!            qbase=MAX(0., sbase/athl+qslth(ind1))
     1959!
     1960!            ! surface cloud fraction in thermals
     1961!            cth(ind1,ind2)=0.5*(1.-erf((sbase-sthl)/sqrt(2.)/sigma2s))
     1962!            cth(ind1,ind2)=MIN(MAX(cth(ind1,ind2),0.),1.)
     1963!
     1964!
     1965!            !volume cloud fraction in thermals
     1966!            !to be checked
     1967!            xth1=-(sthl+deltasth-sbase)/(sqrt(2.)*sigma2s)
     1968!            xth2=-(sthl-deltasth-sbase)/(sqrt(2.)*sigma2s)           
     1969!            exp_xth1 = exp(-1.*xth1**2)
     1970!            exp_xth2 = exp(-1.*xth2**2)
     1971!
     1972!            IntJ=0.5*sthl*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2
     1973!            IntJ_CF=0.5*(1.-1.*erf(xth2))
     1974!     
     1975!            IF (deltasth .LT. 1.e-10) THEN
     1976!              cth_vol(ind1,ind2)=IntJ_CF
     1977!            ELSE
     1978!              IntI1=(((sthl+deltasth-sbase)**2+(sigma2s)**2)/(8*deltasth))*(erf(xth2)-erf(xth1))
     1979!              IntI2=(sigma2s**2/(4*deltasth*sqrtpi))*(xth1*exp_xth1-xth2*exp_xth2)
     1980!              IntI3=((sqrt2*sigma2s*(sthl+deltasth))/(4*sqrtpi*deltasth))*(exp_xth1-exp_xth2)
     1981!              IntI1_CF=((sthl-sbase+deltasth)*(erf(xth2)-erf(xth1)))/(4*deltasth)
     1982!              IntI3_CF=(sqrt2*sigma2s*(exp_xth1-exp_xth2))/(4*sqrtpi*deltasth)
     1983!              cth_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF
     1984!            ENDIF
     1985!              cth_vol(ind1,ind2)=MIN(MAX(0.,cth_vol(ind1,ind2)),1.)
     1986!
     1987!
     1988!
     1989!            ! Environment
     1990!            !=============
     1991!            ! In the environment/downdrafts, only liquid clouds
     1992!            ! See Shupe et al. 2008, JAS
     1993!
     1994!            ! standard deviation of the distribution in the environment
     1995!            sth=sthl
     1996!            senv=senvl
     1997!            sigma1s_fraca = (sigma1s_factor**0.5)*(fraca(ind1,ind2)**sigma1s_power) / &
     1998!                          &                (1-fraca(ind1,ind2))*(MAX((sth-senv),0.)**2)**0.5
     1999!            ! for mixed phase clouds, there is no contribution from large scale ratqs to the distribution
     2000!            ! in the environement
     2001!
     2002!            sigma1s_ratqs=1E-10
     2003!            IF (cloudth_ratqsmin>0.) THEN
     2004!                sigma1s_ratqs = cloudth_ratqsmin*po(ind1)
     2005!            ENDIF
     2006!
     2007!            sigma1s = sigma1s_fraca + sigma1s_ratqs
     2008!            IF (iflag_ratqs.eq.11) then
     2009!               sigma1s = ratqs(ind1,ind2)*po(ind1)*aenv
     2010!            ENDIF
     2011!            IF (iflag_ratqs.eq.11) then
     2012!              sigma1s = ratqs(ind1,ind2)*po(ind1)*aenvl
     2013!            ENDIF
     2014!            deltasenv=aenvl*vert_alpha*sigma1s
     2015!            xenv1=-(senvl+deltasenv)/(sqrt(2.)*sigma1s)
     2016!            xenv2=-(senvl-deltasenv)/(sqrt(2.)*sigma1s)
     2017!            exp_xenv1 = exp(-1.*xenv1**2)
     2018!            exp_xenv2 = exp(-1.*xenv2**2)
     2019!
     2020!            !surface CF
     2021!            cenv(ind1,ind2)=0.5*(1.-1.*erf(xenv1))
     2022!
     2023!            !volume CF and condensed water
     2024!            IntJ=0.5*senvl*(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp_xenv2
     2025!            IntJ_CF=0.5*(1.-1.*erf(xenv2))
     2026!
     2027!            IF (deltasenv .LT. 1.e-10) THEN
     2028!              qcenv(ind1,ind2)=IntJ
     2029!              cenv_vol(ind1,ind2)=IntJ_CF
     2030!            ELSE
     2031!              IntI1=(((senvl+deltasenv)**2+(sigma1s)**2)/(8*deltasenv))*(erf(xenv2)-erf(xenv1))
     2032!              IntI2=(sigma1s**2/(4*deltasenv*sqrtpi))*(xenv1*exp_xenv1-xenv2*exp_xenv2)
     2033!              IntI3=((sqrt2*sigma1s*(senv+deltasenv))/(4*sqrtpi*deltasenv))*(exp_xenv1-exp_xenv2)
     2034!              IntI1_CF=((senvl+deltasenv)*(erf(xenv2)-erf(xenv1)))/(4*deltasenv)
     2035!              IntI3_CF=(sqrt2*sigma1s*(exp_xenv1-exp_xenv2))/(4*sqrtpi*deltasenv)
     2036!              qcenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3 ! only liquid water in environment
     2037!              cenv_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF
     2038!            ENDIF
     2039!
     2040!            qcenv(ind1,ind2)=MAX(qcenv(ind1,ind2),0.)
     2041!            cenv_vol(ind1,ind2)=MIN(MAX(cenv_vol(ind1,ind2),0.),1.)
     2042!
     2043!
     2044!           
     2045!            ! Thermals + environment
     2046!            ! =======================
     2047!            ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)
     2048!            qctot(ind1,ind2)=fraca(ind1,ind2)*qcth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qcenv(ind1,ind2)
     2049!            ctot_vol(ind1,ind2)=fraca(ind1,ind2)*cth_vol(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv_vol(ind1,ind2)
     2050!            IF (qcth(ind1,ind2) .GT. 0) THEN
     2051!               icefrac(ind1,ind2)=fraca(ind1,ind2)*qith(ind1,ind2) &
     2052!                    /(fraca(ind1,ind2)*qcth(ind1,ind2) &
     2053!                    +(1.-1.*fraca(ind1,ind2))*qcenv(ind1,ind2))
     2054!                icefrac(ind1,ind2)=MAX(MIN(1.,icefrac(ind1,ind2)),0.)
     2055!            ELSE
     2056!                icefrac(ind1,ind2)=0.
     2057!            ENDIF
     2058!
     2059!            IF (cenv(ind1,ind2).LT.1.e-10.or.cth(ind1,ind2).LT.1.e-10) THEN
     2060!                ctot(ind1,ind2)=0.
     2061!                ctot_vol(ind1,ind2)=0.
     2062!                qincloud(ind1)=0.
     2063!                qcloud(ind1)=zqsatenv(ind1)
     2064!            ELSE               
     2065!                qcloud(ind1)=fraca(ind1,ind2)*(qcth(ind1,ind2)/cth(ind1,ind2)+qbase) &
     2066!                            +(1.-1.*fraca(ind1,ind2))*(qcenv(ind1,ind2)/cenv(ind1,ind2)+qslenv(ind1))
     2067!                qincloud(ind1)=MAX(fraca(ind1,ind2)*(qcth(ind1,ind2)/cth(ind1,ind2)) &
     2068!                            +(1.-1.*fraca(ind1,ind2))*(qcenv(ind1,ind2)/cenv(ind1,ind2)),0.)
     2069!            ENDIF
     2070!
     2071!        ENDIF ! mpc_bl_points
     2072
     2073
    21522074
    21532075!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  • LMDZ6/trunk/libf/phylmd/lmdz_lscp.F90

    r4879 r4910  
    1212     radocond, radicefrac, rain, snow,                  &
    1313     frac_impa, frac_nucl, beta,                        &
    14      prfl, psfl, rhcl, zqta, fraca,                     &
    15      ztv, zpspsk, ztla, zthl, iflag_cld_th,             &
     14     prfl, psfl, rhcl, qta, fraca,                     &
     15     tv, pspsk, tla, thl, iflag_cld_th,             &
    1616     iflag_ice_thermo, ok_ice_sursat, qsatl, qsats,     &
    1717     distcltop,temp_cltop,                              &
    … …  
    104104USE lmdz_lscp_ini, ONLY : prt_level, lunout
    105105USE lmdz_lscp_ini, ONLY : seuil_neb, niter_lscp, iflag_evap_prec, t_coup, DDT0, ztfondue, rain_int_min
    106 USE lmdz_lscp_ini, ONLY : iflag_mpc_bl, ok_radocond_snow, a_tr_sca, cld_expo_con, cld_expo_lsc
     106USE lmdz_lscp_ini, ONLY : ok_radocond_snow, a_tr_sca, cld_expo_con, cld_expo_lsc
    107107USE lmdz_lscp_ini, ONLY : iflag_cloudth_vert, iflag_rain_incloud_vol, iflag_t_glace, t_glace_min
    108108USE lmdz_lscp_ini, ONLY : coef_eva, coef_sub,cld_tau_lsc, cld_tau_con, cld_lc_lsc, cld_lc_con
    109109USE lmdz_lscp_ini, ONLY : iflag_bergeron, iflag_fisrtilp_qsat, iflag_vice, cice_velo, dice_velo
    110 USE lmdz_lscp_ini, ONLY : iflag_autoconversion, ffallv_con, ffallv_lsc
     110USE lmdz_lscp_ini, ONLY : iflag_autoconversion, ffallv_con, ffallv_lsc, min_frac_th_cld
    111111USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG
    112112USE lmdz_lscp_ini, ONLY : ok_poprecip, ok_bug_fonte_lscp
    … …  
    138138  !Inputs associated with thermal plumes
    139139
    140   REAL, DIMENSION(klon,klev),      INTENT(IN)   :: ztv             ! virtual potential temperature [K]
    141   REAL, DIMENSION(klon,klev),      INTENT(IN)   :: zqta            ! specific humidity within thermals [kg/kg]
    142   REAL, DIMENSION(klon,klev),      INTENT(IN)   :: fraca           ! fraction of thermals within the mesh [-]
    143   REAL, DIMENSION(klon,klev),      INTENT(IN)   :: zpspsk          ! exner potential (p/100000)**(R/cp)
    144   REAL, DIMENSION(klon,klev),      INTENT(IN)   :: ztla            ! liquid temperature within thermals [K]
     140  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: tv             ! virtual potential temperature [K]
     141  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: qta            ! specific humidity within thermals [kg/kg]
     142  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: fraca          ! fraction of thermals within the mesh [-]
     143  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: pspsk          ! exner potential (p/100000)**(R/cp)
     144  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: tla            ! liquid temperature within thermals [K]
    145145
    146146  ! INPUT/OUTPUT variables
    147147  !------------------------
    148148 
    149   REAL, DIMENSION(klon,klev),      INTENT(INOUT)   :: zthl         ! liquid potential temperature [K]
     149  REAL, DIMENSION(klon,klev),      INTENT(INOUT)   :: thl          ! liquid potential temperature [K]
    150150  REAL, DIMENSION(klon,klev),      INTENT(INOUT)   :: ratqs        ! function of pressure that sets the large-scale
    151151
    … …  
    236236  REAL, DIMENSION(klon) :: Zpdf_a,zpdf_b,zpdf_e1,zpdf_e2
    237237  REAL :: erf
    238   REAL, DIMENSION(klon,klev) :: icefrac_mpc
     238  REAL, DIMENSION(klon) :: zfice_th
    239239  REAL, DIMENSION(klon) :: qcloud, qincloud_mpc
    240240  REAL, DIMENSION(klon) :: zrfl, zrfln
    … …  
    713713
    714714        IF (iflag_cld_th.GE.5) THEN
    715 
    716              IF (iflag_mpc_bl .LT. 1) THEN
    717 
    718                 IF (iflag_cloudth_vert.LE.2) THEN
    719 
    720                     CALL cloudth(klon,klev,k,ztv,                  &
    721                          zq,zqta,fraca,                            &
    722                          qcloud,ctot,zpspsk,paprs,pplay,ztla,zthl, &
     715               ! Cloud cover and content in meshes affected by shallow convection,
     716               ! are retrieved from a bi-gaussian distribution of the saturation deficit
     717               ! following Jam et al. 2013
     718
     719               IF (iflag_cloudth_vert.LE.2) THEN
     720                  ! Old version of Arnaud Jam
     721
     722                    CALL cloudth(klon,klev,k,tv,                  &
     723                         zq,qta,fraca,                            &
     724                         qcloud,ctot,pspsk,paprs,pplay,tla,thl, &
    723725                         ratqs,zqs,temp,                              &
    724726                         cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
    … …  
    726728
    727729                ELSEIF (iflag_cloudth_vert.GE.3 .AND. iflag_cloudth_vert.LE.5) THEN
    728 
    729                     CALL cloudth_v3(klon,klev,k,ztv,                        &
    730                          zq,zqta,fraca,                                     &
    731                          qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
     730                   ! Default version of Arnaud Jam
     731                       
     732                    CALL cloudth_v3(klon,klev,k,tv,                        &
     733                         zq,qta,fraca,                                     &
     734                         qcloud,ctot,ctot_vol,pspsk,paprs,pplay,tla,thl, &
    732735                         ratqs,zqs,temp, &
    733736                         cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
    734737
    735                 !Jean Jouhaud's version, Decembre 2018
    736                 !-------------------------------------
    737738
    738739                ELSEIF (iflag_cloudth_vert.EQ.6) THEN
    739 
    740                     CALL cloudth_v6(klon,klev,k,ztv,                        &
    741                          zq,zqta,fraca,                                     &
    742                          qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
     740                   ! Jean Jouhaud's version, with specific separation between surface and volume
     741                   ! cloud fraction Decembre 2018
     742
     743                    CALL cloudth_v6(klon,klev,k,tv,                        &
     744                         zq,qta,fraca,                                     &
     745                         qcloud,ctot,ctot_vol,pspsk,paprs,pplay,tla,thl, &
    743746                         ratqs,zqs,temp, &
    744747                         cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
    745748
    746                 ENDIF
    747 
    748          ELSE
    749                 ! cloudth_v3 + cold microphysical considerations
    750                 ! + stationary mixed-phase cloud model
    751 
    752                     CALL cloudth_mpc(klon,klev,k,iflag_mpc_bl,mpc_bl_points, &
    753                          temp,ztv,zq,zqta,fraca, zpspsk,                      &
    754                          paprs,pplay,ztla,zthl,ratqs,zqs,psfl,             &
    755                          qcloud,qincloud_mpc,icefrac_mpc,ctot,ctot_vol,    &
    756                          cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
    757 
    758          ENDIF ! iflag_mpc_bl
     749                ELSEIF (iflag_cloudth_vert .EQ. 7) THEN
     750                   ! Updated version of Arnaud Jam (correction by E. Vignon) + adapted treatment
     751                   ! for boundary-layer mixed phase clouds following Vignon et al. 
     752                    CALL cloudth_mpc(klon,klev,k,mpc_bl_points,zt,zq,qta(:,k),fraca(:,k), &
     753                                     pspsk(:,k),paprs(:,k+1),paprs(:,k),pplay(:,k), tla(:,k), &
     754                                     ratqs(:,k),qcloud,qincloud_mpc,zfice_th,ctot(:,k),ctot_vol(:,k), &
     755                                     cloudth_sth(:,k),cloudth_senv(:,k),cloudth_sigmath(:,k),cloudth_sigmaenv(:,k))
     756
     757                ENDIF
     758
    759759
    760760                DO i=1,klon
    … …  
    774774           
    775775                ! lognormal distribution when no thermals
    776             lognormale = fraca(:,k) < 1e-10
     776            lognormale = fraca(:,k) < min_frac_th_cld
    777777
    778778        ELSE
    … …  
    934934            temp_cltop(:,k)=temp_cltop1D(:)
    935935        ENDIF   
    936 
    937936        ! Partition function in stratiform clouds (will be overwritten in boundary-layer MPCs)
    938         CALL icefrac_lscp(klon,zt,iflag_ice_thermo,distcltop1D,temp_cltop1D,zfice, dzfice)
     937        CALL icefrac_lscp(klon,zt,iflag_ice_thermo,distcltop1D,temp_cltop1D,zfice,dzfice)
    939938
    940939
    … …  
    943942
    944943            IF (mpc_bl_points(i,k) .GT. 0) THEN
     944               
    945945                zcond(i) = MAX(0.0,qincloud_mpc(i))*rneb(i,k)
    946946                ! following line is very strange and probably wrong
    … …  
    948948                ! water vapor update and partition function if thermals
    949949                zq(i) = zq(i) - zcond(i)       
    950                 zfice(i)=icefrac_mpc(i,k)
     950                zfice(i)=zfice_th(i)
    951951
    952952            ELSE
    … …  
    13141314    ENDIF
    13151315
     1316
    13161317    ENDIF ! ok_poprecip
    13171318

  • LMDZ6/trunk/libf/phylmd/lmdz_lscp_ini.F90

    r4898 r4910  
    1212  !$OMP THREADPRIVATE(ok_bug_fonte_lscp)
    1313
    14   REAL, SAVE, PROTECTED :: seuil_neb=0.001      ! cloud fraction threshold: a cloud really exists when exceeded
     14  REAL, SAVE, PROTECTED :: seuil_neb=0.001      ! cloud fraction threshold: a cloud can precipitate when exceeded
    1515  !$OMP THREADPRIVATE(seuil_neb)
     16
     17  REAL, SAVE, PROTECTED :: min_neb_th=1e-10      ! a cloud produced by bi-gaussian really exists when exceeded
     18  !$OMP THREADPRIVATE(min_neb_th)
     19
     20  REAL, SAVE, PROTECTED :: min_frac_thermals=1.e-10   ! minimum thermals fraction for use of bigaussian distribution
     21  !$OMP THREADPRIVATE(min_frac_thermals)
    1622
    1723  INTEGER, SAVE :: lunout, prt_level            ! Logical unit number and level for standard output
    … …  
    4349  !$OMP THREADPRIVATE(a_tr_sca)
    4450 
    45   INTEGER, SAVE, PROTECTED ::  iflag_mpc_bl=0   ! flag to activate boundary layer mixed phase cloud param
    46   !$OMP THREADPRIVATE(iflag_mpc_bl)
    47  
     51  REAL, SAVE, PROTECTED ::  min_frac_th_cld=1.e-10   ! minimum thermal fraction to compute a thermal cloud fraction
     52  !$OMP THREADPRIVATE(min_frac_th_cld)
     53
    4854  LOGICAL, SAVE, PROTECTED :: ok_radocond_snow=.false. ! take into account the mass of ice precip in the cloud ice content seen by radiation
    4955  !$OMP THREADPRIVATE(ok_radocond_snow)
    … …  
    260266    CALL getin_p('iflag_evap_prec',iflag_evap_prec)
    261267    CALL getin_p('seuil_neb',seuil_neb)
    262     CALL getin_p('iflag_mpc_bl',iflag_mpc_bl)
    263268    CALL getin_p('ok_radocond_snow',ok_radocond_snow)
    264269    CALL getin_p('t_glace_max',t_glace_max)
    … …  
    320325    WRITE(lunout,*) 'lscp_ini, iflag_evap_prec:', iflag_evap_prec
    321326    WRITE(lunout,*) 'lscp_ini, seuil_neb:', seuil_neb
    322     WRITE(lunout,*) 'lscp_ini, iflag_mpc_bl:', iflag_mpc_bl
    323327    WRITE(lunout,*) 'lscp_ini, ok_radocond_snow:', ok_radocond_snow
    324328    WRITE(lunout,*) 'lscp_ini, t_glace_max:', t_glace_max
Note: See TracChangeset for help on using the changeset viewer.