Changeset 3999


Ignore:
Timestamp:
Nov 5, 2021, 12:40:08 PM (3 years ago)
Author:
evignon
Message:

commission de la nouvelle routine de condensation
grande echelle simplifiee (lscp, version epuree de fistilp)
et du schema de nuages de phase mixte (en developpement)

La routine lscp n'est active que sous flag
ok_new_lscp=y

Location:
LMDZ6/trunk/libf/phylmd
Files:
2 added
8 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/trunk/libf/phylmd/clesphys.h

    r3435 r3999  
    9393       LOGICAL :: adjust_tropopause
    9494       LOGICAL :: ok_daily_climoz
     95       LOGICAL :: ok_new_lscp
    9596! flag to bypass or not the phytrac module
    9697       INTEGER :: iflag_phytrac
    … …  
    141142     &     , ok_chlorophyll,ok_conserv_q, adjust_tropopause             &
    142143     &     , ok_daily_climoz, ok_all_xml, ok_lwoff                      &
    143      &     , iflag_phytrac
     144     &     , iflag_phytrac, ok_new_lscp
    144145     
    145146       save /clesphys/

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

    r3570 r3999  
    655655      REAL zqs(ngrid), qcloud(ngrid)
    656656      REAL erf
    657 
    658657
    659658
    … …  
    911910      END DO
    912911
    913 
    914912!------------------------------------------------------------------------------
    915913! Initialize
    916914!------------------------------------------------------------------------------
     915
    917916      sigma1(:,:)=0.
    918917      sigma2(:,:)=0.
    … …  
    10131012!zqsatth = qsat thermals
    10141013!ztla = Tl thermals
    1015 
    10161014!------------------------------------------------------------------------------
    10171015! s standard deviation
    … …  
    12171215      else  ! gaussienne environnement seule
    12181216     
     1217
    12191218      zqenv(ind1)=po(ind1)
    12201219      Tbef=t(ind1,ind2)
    … …  
    15341533
    15351534END SUBROUTINE cloudth_v6
     1535
     1536
     1537
     1538
     1539
     1540!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     1541SUBROUTINE cloudth_mpc(klon,klev,ind2,mpc_bl_points,                     &
     1542&           temp,ztv,po,zqta,fraca,zpspsk,paprs,pplay,ztla,zthl,            &
     1543&           ratqs,zqs,snowflux,qcloud,qincloud,icefrac,ctot,ctot_vol)
     1544!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     1545! Author : Arnaud Octavio Jam (LMD/CNRS), Etienne Vignon (LMDZ/CNRS)
     1546! Date: Adapted from cloudth_vert_v3 in 2021
     1547! Aim : computes qc and rneb in thermals with cold microphysical considerations
     1548!       + for mixed phase boundary layer clouds, calculate ql and qi from
     1549!       a stationary MPC model
     1550! IMPORTANT NOTE: we assume iflag_clouth_vert=3
     1551!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     1552
     1553
     1554      USE ioipsl_getin_p_mod, ONLY : getin_p
     1555      USE phys_output_var_mod, ONLY : cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv
     1556      USE lscp_tools_mod, ONLY: CALC_QSAT_ECMWF, ICEFRAC_LSCP
     1557      USE phys_local_var_mod, ONLY : qlth, qith
     1558
     1559      IMPLICIT NONE
     1560
     1561#include "YOMCST.h"
     1562#include "YOETHF.h"
     1563#include "FCTTRE.h"
     1564#include "thermcell.h"
     1565#include "nuage.h"
     1566     
     1567
     1568!------------------------------------------------------------------------------
     1569! Declaration
     1570!------------------------------------------------------------------------------
     1571
     1572! INPUT/OUTPUT
     1573
     1574      INTEGER, INTENT(IN)                         :: klon,klev,ind2
     1575      INTEGER, DIMENSION(klon,klev), INTENT(INOUT)   ::  mpc_bl_points ! 1 where BL MPC, 0 otherwise
     1576
     1577      REAL, DIMENSION(klon,klev), INTENT(IN)      ::  temp          ! Temperature [K]
     1578      REAL, DIMENSION(klon,klev), INTENT(IN)      ::  ztv           ! Virtual potential temp [K]
     1579      REAL, DIMENSION(klon),      INTENT(IN)      ::  po            ! specific humidity [kg/kg]
     1580      REAL, DIMENSION(klon,klev), INTENT(IN)      ::  zqta          ! specific humidity within thermals [kg/kg]
     1581      REAL, DIMENSION(klon,klev+1), INTENT(IN)    ::  fraca         ! Fraction of the mesh covered by thermals [0-1]
     1582      REAL, DIMENSION(klon,klev), INTENT(IN)      ::  zpspsk
     1583      REAL, DIMENSION(klon,klev+1), INTENT(IN)    ::  paprs         ! Pressure at layer interfaces [Pa]
     1584      REAL, DIMENSION(klon,klev), INTENT(IN)      ::  pplay         ! Pressure at the center of layers [Pa]
     1585      REAL, DIMENSION(klon,klev), INTENT(IN)      ::  ztla          ! Liquid temp [K]
     1586      REAL, DIMENSION(klon,klev), INTENT(INOUT)      ::  zthl       ! Liquid potential temp [K]
     1587      REAL, DIMENSION(klon,klev), INTENT(IN)      ::  ratqs         ! Parameter that determines the width of the total water distrib.
     1588      REAL, DIMENSION(klon),      INTENT(IN)      ::  zqs           ! Saturation specific humidity in the mesh [kg/kg]
     1589      REAL, DIMENSION(klon,klev+1), INTENT(IN)    ::  snowflux      ! snow flux at the interface of the layer [kg/m2/s]
     1590
     1591
     1592      REAL, DIMENSION(klon,klev), INTENT(OUT)      ::  ctot         ! Cloud fraction [0-1]
     1593      REAL, DIMENSION(klon,klev), INTENT(OUT)      ::  ctot_vol     ! Volume cloud fraction [0-1]
     1594      REAL, DIMENSION(klon),      INTENT(OUT)      ::  qcloud       ! In cloud total water content [kg/kg]
     1595      REAL, DIMENSION(klon),      INTENT(OUT)      ::  qincloud       ! In cloud condensed water content [kg/kg]
     1596      REAL, DIMENSION(klon,klev), INTENT(OUT)      ::  icefrac      ! Fraction of ice in clouds [0-1]
     1597
     1598
     1599! LOCAL VARIABLES
     1600
     1601      INTEGER itap,ind1,l,ig,iter,k
     1602      LOGICAL flag_topthermals
     1603
     1604
     1605      REAL zqsatth(klon,klev), zqsatenv(klon,klev)
     1606      REAL sigma1(klon,klev)                                                         
     1607      REAL sigma2(klon,klev)
     1608      REAL qcth(klon,klev)
     1609      REAL qcenv(klon,klev)
     1610      REAL qctot(klon,klev)
     1611      REAL cth(klon,klev)
     1612      REAL cenv(klon,klev)   
     1613      REAL cth_vol(klon,klev)
     1614      REAL cenv_vol(klon,klev)
     1615      REAL rneb(klon,klev)
     1616      REAL zqenv(klon)   
     1617
     1618      REAL qsatmmussig1,qsatmmussig2,sqrtpi,sqrt2,sqrt2pi,pi
     1619      REAL rdd,cppd,Lv
     1620      REAL alth,alenv,ath,aenv
     1621      REAL sth,senv,sigma1s,sigma2s,sigma1s_fraca,sigma1s_ratqs
     1622      REAL inverse_rho,beta,a_Brooks,b_Brooks,A_Maj_Brooks,Dx_Brooks,f_Brooks
     1623      REAL xth,xenv,exp_xenv1,exp_xenv2,exp_xth1,exp_xth2
     1624      REAL xth1,xth2,xenv1,xenv2,deltasth, deltasenv
     1625      REAL IntJ,IntI1,IntI2,IntI3,IntJ_CF,IntI1_CF,IntI3_CF,coeffqlenv,coeffqlth
     1626      REAL Tbef,zdelta,qsatbef,zcor
     1627      REAL qlbef,dqsatdt
     1628      REAL erf
     1629      REAL zpdf_sig(klon),zpdf_k(klon),zpdf_delta(klon)
     1630      REAL zpdf_a(klon),zpdf_b(klon),zpdf_e1(klon),zpdf_e2(klon)
     1631      REAL rhodz(klon,klev)
     1632      REAL zrho(klon,klev)
     1633      REAL dz(klon,klev)
     1634      REAL qslth, qsith, qslenv, alenvl, aenvl
     1635      REAL sthi, sthl, althl, athl
     1636      REAL senvi, senvl, qbase, sbase, qliqth, qiceth
     1637      REAL qimax, ttarget, stmp, cout, coutref
     1638      REAL maxi, mini, pas, temp_lim
     1639      REAL deltazlev_mpc(klev),qth_mpc(klev), temp_mpc(klev), pres_mpc(klev), fraca_mpc(klev+1), snowf_mpc(klev+1)
     1640
     1641      INTEGER, SAVE :: niter=20
     1642
     1643      ! Modifty the saturation deficit PDF in thermals
     1644      ! in the presence of ice crystals
     1645      REAL,SAVE :: C_mpc
     1646      !$OMP THREADPRIVATE(C_mpc)
     1647      ! Change the width of the PDF used for vertical subgrid scale heterogeneity
     1648      ! (J Jouhaud, JL Dufresne, JB Madeleine)
     1649      REAL,SAVE :: vert_alpha, vert_alpha_th
     1650      !$OMP THREADPRIVATE(vert_alpha, vert_alpha_th)
     1651      REAL,SAVE :: sigma1s_factor=1.1
     1652      REAL,SAVE :: sigma1s_power=0.6
     1653      REAL,SAVE :: sigma2s_factor=0.09
     1654      REAL,SAVE :: sigma2s_power=0.5
     1655      REAL,SAVE :: cloudth_ratqsmin=-1.
     1656      !$OMP THREADPRIVATE(sigma1s_factor,sigma1s_power,sigma2s_factor,sigma2s_power,cloudth_ratqsmin)
     1657      INTEGER, SAVE :: iflag_cloudth_vert_noratqs=0
     1658      !$OMP THREADPRIVATE(iflag_cloudth_vert_noratqs)
     1659      LOGICAL, SAVE :: firstcall = .TRUE.
     1660      !$OMP THREADPRIVATE(firstcall)
     1661
     1662      CHARACTER (len = 80) :: abort_message
     1663      CHARACTER (len = 20) :: routname = 'cloudth_mpc'
     1664
     1665
     1666!------------------------------------------------------------------------------
     1667! Initialisation
     1668!------------------------------------------------------------------------------
     1669
     1670
     1671! Few initial checksS
     1672
     1673      IF (iflag_cloudth_vert.NE.3) THEN
     1674         abort_message = 'clouth_mpc cannot be used if iflag_cloudth_vert .NE. 3'
     1675         CALL abort_physic(routname,abort_message,1)
     1676      ENDIF
     1677
     1678      DO k = 1,klev
     1679      DO ind1 = 1, klon
     1680        rhodz(ind1,k) = (paprs(ind1,k)-paprs(ind1,k+1))/rg !kg/m2
     1681        zrho(ind1,k) = pplay(ind1,k)/temp(ind1,k)/rd !kg/m3
     1682        dz(ind1,k) = rhodz(ind1,k)/zrho(ind1,k) !m : epaisseur de la couche en metre
     1683      END DO
     1684      END DO
     1685
     1686
     1687      sigma1(:,:)=0.
     1688      sigma2(:,:)=0.
     1689      qcth(:,:)=0.
     1690      qcenv(:,:)=0. 
     1691      qctot(:,:)=0.
     1692      qlth(:,ind2)=0.
     1693      qith(:,ind2)=0.
     1694      rneb(:,:)=0.
     1695      qcloud(:)=0.
     1696      cth(:,:)=0.
     1697      cenv(:,:)=0.
     1698      ctot(:,:)=0.
     1699      cth_vol(:,:)=0.
     1700      cenv_vol(:,:)=0.
     1701      ctot_vol(:,:)=0.
     1702      qsatmmussig1=0.
     1703      qsatmmussig2=0.
     1704      rdd=287.04
     1705      cppd=1005.7
     1706      pi=3.14159
     1707      sqrt2pi=sqrt(2.*pi)
     1708      sqrt2=sqrt(2.)
     1709      sqrtpi=sqrt(pi)
     1710      icefrac(:,ind2)=0.
     1711
     1712
     1713
     1714      IF (firstcall) THEN
     1715
     1716        vert_alpha=0.5
     1717        CALL getin_p('cloudth_vert_alpha',vert_alpha)
     1718        WRITE(*,*) 'cloudth_vert_alpha = ', vert_alpha
     1719        ! The factor used for the thermal is equal to that of the environment
     1720        !   if nothing is explicitly specified in the def file
     1721        vert_alpha_th=vert_alpha
     1722        CALL getin_p('cloudth_vert_alpha_th',vert_alpha_th)
     1723        WRITE(*,*) 'cloudth_vert_alpha_th = ', vert_alpha_th
     1724        ! Factor used in the calculation of sigma1s
     1725        CALL getin_p('cloudth_sigma1s_factor',sigma1s_factor)
     1726        WRITE(*,*) 'cloudth_sigma1s_factor = ', sigma1s_factor
     1727        ! Power used in the calculation of sigma1s
     1728        CALL getin_p('cloudth_sigma1s_power',sigma1s_power)
     1729        WRITE(*,*) 'cloudth_sigma1s_power = ', sigma1s_power
     1730        ! Factor used in the calculation of sigma2s
     1731        CALL getin_p('cloudth_sigma2s_factor',sigma2s_factor)
     1732        WRITE(*,*) 'cloudth_sigma2s_factor = ', sigma2s_factor
     1733        ! Power used in the calculation of sigma2s
     1734        CALL getin_p('cloudth_sigma2s_power',sigma2s_power)
     1735        WRITE(*,*) 'cloudth_sigma2s_power = ', sigma2s_power
     1736        ! Minimum value for the environmental air subgrid water distrib
     1737        CALL getin_p('cloudth_ratqsmin',cloudth_ratqsmin)
     1738        WRITE(*,*) 'cloudth_ratqsmin = ', cloudth_ratqsmin
     1739        ! Remove the dependency to ratqs from the variance of the vertical PDF
     1740        CALL getin_p('iflag_cloudth_vert_noratqs',iflag_cloudth_vert_noratqs)
     1741        WRITE(*,*) 'iflag_cloudth_vert_noratqs = ', iflag_cloudth_vert_noratqs
     1742        ! Modifies the PDF in thermals when ice crystals are present
     1743        C_mpc=1.e2
     1744        CALL getin_p('C_mpc',C_mpc)
     1745        WRITE(*,*) 'C_mpc = ', C_mpc
     1746
     1747        firstcall=.FALSE.
     1748
     1749      ENDIF
     1750
     1751
     1752
     1753!-------------------------------------------------------------------------------
     1754! Identify grid points with potential mixed-phase conditions
     1755!-------------------------------------------------------------------------------
     1756
     1757      temp_lim=RTT-40.0
     1758
     1759      DO ind1=1,klon
     1760            IF ((temp(ind1,ind2) .LT. RTT) .AND. (temp(ind1,ind2) .GT. temp_lim) &
     1761            .AND. (iflag_mpc_bl .GE. 2) .AND. (ind2<=klev-2)  &
     1762            .AND. (ztv(ind1,1).GT.ztv(ind1,2)) .AND.(fraca(ind1,ind2).GT.1.e-10)) THEN
     1763                mpc_bl_points(ind1,ind2)=1
     1764            ELSE
     1765                mpc_bl_points(ind1,ind2)=0
     1766            ENDIF
     1767      ENDDO
     1768
     1769
     1770!-------------------------------------------------------------------------------
     1771! Thermal fraction calculation and standard deviation of the distribution
     1772!------------------------------------------------------------------------------- 
     1773
     1774    DO ind1=1,klon
     1775
     1776
     1777    IF ((ztv(ind1,1).GT.ztv(ind1,2)).AND.(fraca(ind1,ind2).GT.1.e-10)) THEN !Thermal and environnement
     1778
     1779
     1780! Environment:
     1781
     1782        zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2)) !qt = a*qtth + (1-a)*qtenv
     1783        Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
     1784           
     1785        CALL CALC_QSAT_ECMWF(Tbef,0.,paprs(ind1,ind2),RTT,0,.false.,qsatbef,dqsatdt)
     1786        zqsatenv(ind1,ind2)=qsatbef
     1787
     1788        IF (Tbef .GE. RTT) THEN
     1789               Lv=RLVTT
     1790        ELSE
     1791               Lv=RLSTT
     1792        ENDIF
     1793       
     1794
     1795        alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)     !qsl, p84
     1796        aenv=1./(1.+(alenv*Lv/cppd))                                      !al, p84
     1797        senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))                          !s, p84
     1798     
     1799        ! For MPCs:
     1800        IF (mpc_bl_points(ind1,ind2) .EQ. 1) THEN
     1801        CALL CALC_QSAT_ECMWF(Tbef,0.,paprs(ind1,ind2),RTT,1,.false.,qslenv,dqsatdt)
     1802        alenvl=(0.622*RLVTT*qslenv)/(rdd*zthl(ind1,ind2)**2)     
     1803        aenvl=1./(1.+(alenv*Lv/cppd))         
     1804        senvl=aenvl*(po(ind1)-qslenv)   
     1805        ENDIF
     1806
     1807
     1808! Thermals:
     1809
     1810        Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
     1811        CALL CALC_QSAT_ECMWF(Tbef,0.,paprs(ind1,ind2),RTT,0,.false.,qsatbef,dqsatdt)
     1812        zqsatth(ind1,ind2)=qsatbef
     1813
     1814        IF (Tbef .GE. RTT) THEN
     1815            Lv=RLVTT
     1816        ELSE
     1817            Lv=RLSTT
     1818        ENDIF
     1819
     1820       
     1821        alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2)       
     1822        ath=1./(1.+(alth*Lv/cppd))                                                         
     1823        sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2))                     
     1824
     1825       ! For MPCs:
     1826        IF (mpc_bl_points(ind1,ind2) .GT. 0) THEN
     1827         CALL CALC_QSAT_ECMWF(Tbef,0.,paprs(ind1,ind2),RTT,1,.false.,qslth,dqsatdt)
     1828         CALL CALC_QSAT_ECMWF(Tbef,0.,paprs(ind1,ind2),RTT,2,.false.,qsith,dqsatdt)
     1829         althl=(0.622*RLVTT*qslth)/(rdd*ztla(ind1,ind2)**2)                   
     1830         athl=1./(1.+(alth*RLVTT/cppd))                                       
     1831         sthl=athl*(zqta(ind1,ind2)-qslth)   
     1832         sthi=athl*(zqta(ind1,ind2)-qsith) 
     1833        ENDIF     
     1834
     1835
     1836
     1837!-------------------------------------------------------------------------------
     1838!  Version 3: Changes by J. Jouhaud; condensation for q > -delta s
     1839!  Rq: in this subroutine, we assume iflag_clouth_vert .EQ. 3
     1840!-------------------------------------------------------------------------------
     1841
     1842        IF (mpc_bl_points(ind1,ind2) .EQ. 0) THEN ! No BL MPC
     1843
     1844       ! Standard deviation of the distributions
     1845
     1846           sigma1s_fraca = (sigma1s_factor**0.5)*(fraca(ind1,ind2)**sigma1s_power) / &
     1847           &                (1-fraca(ind1,ind2))*((sth-senv)**2)**0.5
     1848
     1849           IF (cloudth_ratqsmin>0.) THEN
     1850             sigma1s_ratqs = cloudth_ratqsmin*po(ind1)
     1851           ELSE
     1852             sigma1s_ratqs = ratqs(ind1,ind2)*po(ind1)
     1853           ENDIF
     1854 
     1855           sigma1s = sigma1s_fraca + sigma1s_ratqs
     1856           sigma2s=(sigma2s_factor*(((sth-senv)**2)**0.5)/((fraca(ind1,ind2)+0.02)**sigma2s_power))+0.002*zqta(ind1,ind2)
     1857
     1858
     1859           deltasenv=aenv*vert_alpha*sigma1s
     1860           deltasth=ath*vert_alpha_th*sigma2s
     1861
     1862           xenv1=-(senv+deltasenv)/(sqrt(2.)*sigma1s)
     1863           xenv2=-(senv-deltasenv)/(sqrt(2.)*sigma1s)
     1864           exp_xenv1 = exp(-1.*xenv1**2)
     1865           exp_xenv2 = exp(-1.*xenv2**2)
     1866           xth1=-(sth+deltasth)/(sqrt(2.)*sigma2s)
     1867           xth2=-(sth-deltasth)/(sqrt(2.)*sigma2s)
     1868           exp_xth1 = exp(-1.*xth1**2)
     1869           exp_xth2 = exp(-1.*xth2**2)
     1870     
     1871      !surface CF
     1872
     1873           cth(ind1,ind2)=0.5*(1.-1.*erf(xth1))
     1874           cenv(ind1,ind2)=0.5*(1.-1.*erf(xenv1))
     1875           ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)
     1876
     1877
     1878      !volume CF and condensed water
     1879
     1880            !environnement
     1881
     1882            IntJ=0.5*senv*(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp_xenv2
     1883            IntJ_CF=0.5*(1.-1.*erf(xenv2))
     1884
     1885            IF (deltasenv .LT. 1.e-10) THEN
     1886              qcenv(ind1,ind2)=IntJ
     1887              cenv_vol(ind1,ind2)=IntJ_CF
     1888            ELSE
     1889              IntI1=(((senv+deltasenv)**2+(sigma1s)**2)/(8*deltasenv))*(erf(xenv2)-erf(xenv1))
     1890              IntI2=(sigma1s**2/(4*deltasenv*sqrtpi))*(xenv1*exp_xenv1-xenv2*exp_xenv2)
     1891              IntI3=((sqrt2*sigma1s*(senv+deltasenv))/(4*sqrtpi*deltasenv))*(exp_xenv1-exp_xenv2)
     1892              IntI1_CF=((senv+deltasenv)*(erf(xenv2)-erf(xenv1)))/(4*deltasenv)
     1893              IntI3_CF=(sqrt2*sigma1s*(exp_xenv1-exp_xenv2))/(4*sqrtpi*deltasenv)
     1894              qcenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
     1895              cenv_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF
     1896            ENDIF
     1897             
     1898
     1899
     1900            !thermals
     1901
     1902            IntJ=0.5*sth*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2
     1903            IntJ_CF=0.5*(1.-1.*erf(xth2))
     1904     
     1905            IF (deltasth .LT. 1.e-10) THEN
     1906              qcth(ind1,ind2)=IntJ
     1907              cth_vol(ind1,ind2)=IntJ_CF
     1908            ELSE
     1909              IntI1=(((sth+deltasth)**2+(sigma2s)**2)/(8*deltasth))*(erf(xth2)-erf(xth1))
     1910              IntI2=(sigma2s**2/(4*deltasth*sqrtpi))*(xth1*exp_xth1-xth2*exp_xth2)
     1911              IntI3=((sqrt2*sigma2s*(sth+deltasth))/(4*sqrtpi*deltasth))*(exp_xth1-exp_xth2)
     1912              IntI1_CF=((sth+deltasth)*(erf(xth2)-erf(xth1)))/(4*deltasth)
     1913              IntI3_CF=(sqrt2*sigma2s*(exp_xth1-exp_xth2))/(4*sqrtpi*deltasth)
     1914              qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
     1915              cth_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF
     1916            ENDIF
     1917
     1918              qctot(ind1,ind2)=fraca(ind1,ind2)*qcth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qcenv(ind1,ind2)
     1919              ctot_vol(ind1,ind2)=fraca(ind1,ind2)*cth_vol(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv_vol(ind1,ind2)
     1920             
     1921
     1922            IF (cenv(ind1,ind2).LT.1.e-10.or.cth(ind1,ind2).LT.1.e-10) THEN
     1923                ctot(ind1,ind2)=0.
     1924                ctot_vol(ind1,ind2)=0.
     1925                qcloud(ind1)=zqsatenv(ind1,ind2)
     1926                qincloud(ind1)=0.
     1927            ELSE               
     1928                qcloud(ind1)=qctot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
     1929                qincloud(ind1)=qctot(ind1,ind2)/ctot(ind1,ind2)
     1930            ENDIF
     1931
     1932
     1933        ELSE ! mpc_bl_points>0
     1934
     1935            ! Treat boundary layer mixed phase clouds
     1936           
     1937            ! thermals
     1938            !=========
     1939
     1940            ! ice phase
     1941            !...........
     1942
     1943            deltazlev_mpc=dz(ind1,:)
     1944            temp_mpc=ztla(ind1,:)*zpspsk(ind1,:)
     1945            pres_mpc=pplay(ind1,:)
     1946            fraca_mpc=fraca(ind1,:)
     1947            snowf_mpc=snowflux(ind1,:)
     1948            qth_mpc=zqta(ind1,:)
     1949            flag_topthermals=.FALSE.
     1950            IF ((mpc_bl_points(ind1,ind2) .EQ. 1) .AND. (mpc_bl_points(ind1,ind2+1) .EQ. 0))  THEN
     1951                flag_topthermals = .TRUE.
     1952            ENDIF
     1953
     1954            CALL ICE_MPC_BL_CLOUDS(ind1,ind2,klev,flag_topthermals,temp_mpc,pres_mpc,qth_mpc,qlth(ind1,:),qith(ind1,:),deltazlev_mpc,snowf_mpc,fraca_mpc,qiceth)
     1955
     1956
     1957
     1958            ! We account for the effect of ice crystals in thermals on sthl
     1959            ! and on the width of the distribution
     1960
     1961            sthl=sthl*1./(1.+C_mpc*qiceth)  &
     1962                + (1.-1./(1.+C_mpc*qiceth)) * athl*(zqta(ind1,ind2)-(qsith+qiceth)) 
     1963
     1964            sthi=sthi*1./(1.+C_mpc*qiceth)  &
     1965                + (1.-1./(1.+C_mpc*qiceth)) * athl*(zqta(ind1,ind2)-(qsith+qiceth)) 
     1966
     1967           ! standard deviation of the water distribution in thermals
     1968            sth=sthl
     1969            senv=senvl
     1970            sigma2s=(sigma2s_factor*((MAX((sth-senv),0.)**2)**0.5)/((fraca(ind1,ind2)+0.02)**sigma2s_power))+0.002*zqta(ind1,ind2)
     1971            deltasth=athl*vert_alpha_th*sigma2s
     1972     
     1973            ! Liquid phase
     1974            !.............
     1975            xth1=-(sthl+deltasth)/(sqrt(2.)*sigma2s)
     1976            xth2=-(sthl-deltasth)/(sqrt(2.)*sigma2s)           
     1977            exp_xth1 = exp(-1.*xth1**2)
     1978            exp_xth2 = exp(-1.*xth2**2)
     1979            IntJ=0.5*sthl*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2
     1980            IntJ_CF=0.5*(1.-1.*erf(xth2))
     1981            IntI1=(((sthl+deltasth)**2+(sigma2s)**2)/(8*deltasth))*(erf(xth2)-erf(xth1))
     1982            IntI2=(sigma2s**2/(4*deltasth*sqrtpi))*(xth1*exp_xth1-xth2*exp_xth2)
     1983            IntI3=((sqrt2*sigma2s*(sthl+deltasth))/(4*sqrtpi*deltasth))*(exp_xth1-exp_xth2)
     1984            IntI1_CF=((sthl+deltasth)*(erf(xth2)-erf(xth1)))/(4*deltasth)
     1985            IntI3_CF=(sqrt2*sigma2s*(exp_xth1-exp_xth2))/(4*sqrtpi*deltasth)
     1986            qliqth=IntJ+IntI1+IntI2+IntI3
     1987
     1988            ! qimax calculation
     1989            xth1=-(sthi+deltasth)/(sqrt(2.)*sigma2s)
     1990            xth2=-(sthi-deltasth)/(sqrt(2.)*sigma2s)           
     1991            exp_xth1 = exp(-1.*xth1**2)
     1992            exp_xth2 = exp(-1.*xth2**2)
     1993            IntJ=0.5*sthi*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2
     1994            IntJ_CF=0.5*(1.-1.*erf(xth2))
     1995            IntI1=(((sthi+deltasth)**2+(sigma2s)**2)/(8*deltasth))*(erf(xth2)-erf(xth1))
     1996            IntI2=(sigma2s**2/(4*deltasth*sqrtpi))*(xth1*exp_xth1-xth2*exp_xth2)
     1997            IntI3=((sqrt2*sigma2s*(sthi+deltasth))/(4*sqrtpi*deltasth))*(exp_xth1-exp_xth2)
     1998            IntI1_CF=((sthi+deltasth)*(erf(xth2)-erf(xth1)))/(4*deltasth)
     1999            IntI3_CF=(sqrt2*sigma2s*(exp_xth1-exp_xth2))/(4*sqrtpi*deltasth)
     2000            qimax=IntJ+IntI1+IntI2+IntI3
     2001            qimax=qimax-qliqth
     2002
     2003
     2004            ! Condensed water
     2005            ! Guarantee the consistency between qiceth and the subgrid scale PDF of total water
     2006            qlth(ind1,ind2)=MAX(0.,qliqth)
     2007            qith(ind1,ind2)=MAX(0.,MIN(qiceth,qimax))
     2008            qcth(ind1,ind2)=qlth(ind1,ind2)+qith(ind1,ind2)
     2009
     2010            ! calculation of qbase which is the value of the water vapor within mixed phase clouds
     2011            ! such that the total water in cloud = qbase+qliqth+qiceth
     2012            ! sbase is the value of s such that int_sbase^\intfy s ds = cloud fraction
     2013            ! sbase and qbase calculation (note that sbase is wrt liq so negative)
     2014            ! look for an approximate solution with iteration
     2015           
     2016            ttarget=qcth(ind1,ind2)
     2017            mini=athl*(qsith-qslth)
     2018            maxi=0.
     2019            pas=(maxi-mini)/niter
     2020            stmp=mini
     2021            sbase=stmp
     2022            coutref=1.E6
     2023            DO iter=1,niter
     2024                cout=ABS(sigma2s/SQRT(2.*RPI)*EXP(-((sthl-stmp)/sigma2s)**2)+(sthl-stmp)/SQRT(2.)*(1.-erf(-(sthl-stmp)/sigma2s)) &
     2025                     + stmp/2.*(1.-erf(-(sthl-stmp)/sigma2s)) -ttarget)
     2026               IF (cout .LT. coutref) THEN
     2027                     sbase=stmp
     2028                     coutref=cout
     2029                ELSE
     2030                     stmp=stmp+pas
     2031                ENDIF
     2032            ENDDO
     2033            qbase=MAX(0., sbase/athl+qslth)
     2034
     2035            ! surface cloud fraction in thermals
     2036            cth(ind1,ind2)=0.5*(1.-erf((sbase-sthl)/sqrt(2.)/sigma2s))
     2037            cth(ind1,ind2)=MIN(MAX(cth(ind1,ind2),0.),1.)
     2038
     2039
     2040            !volume cloud fraction in thermals
     2041            !to be checked
     2042            xth1=-(sthl+deltasth-sbase)/(sqrt(2.)*sigma2s)
     2043            xth2=-(sthl-deltasth-sbase)/(sqrt(2.)*sigma2s)           
     2044            exp_xth1 = exp(-1.*xth1**2)
     2045            exp_xth2 = exp(-1.*xth2**2)
     2046
     2047            IntJ=0.5*sthl*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2
     2048            IntJ_CF=0.5*(1.-1.*erf(xth2))
     2049     
     2050            IF (deltasth .LT. 1.e-10) THEN
     2051              cth_vol(ind1,ind2)=IntJ_CF
     2052            ELSE
     2053              IntI1=(((sthl+deltasth-sbase)**2+(sigma2s)**2)/(8*deltasth))*(erf(xth2)-erf(xth1))
     2054              IntI2=(sigma2s**2/(4*deltasth*sqrtpi))*(xth1*exp_xth1-xth2*exp_xth2)
     2055              IntI3=((sqrt2*sigma2s*(sth+deltasth))/(4*sqrtpi*deltasth))*(exp_xth1-exp_xth2)
     2056              IntI1_CF=((sthl-sbase+deltasth)*(erf(xth2)-erf(xth1)))/(4*deltasth)
     2057              IntI3_CF=(sqrt2*sigma2s*(exp_xth1-exp_xth2))/(4*sqrtpi*deltasth)
     2058              cth_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF
     2059            ENDIF
     2060              cth_vol(ind1,ind2)=MIN(MAX(0.,cth_vol(ind1,ind2)),1.)
     2061
     2062            ! Environment
     2063            !=============
     2064            ! In the environment/downdrafts, only liquid clouds
     2065            ! See Shupe et al. 2008, JAS
     2066
     2067            ! standard deviation of the distribution in the environment
     2068            sth=sthl
     2069            senv=senvl
     2070            sigma1s_fraca = (sigma1s_factor**0.5)*(fraca(ind1,ind2)**sigma1s_power) / &
     2071                          &                (1-fraca(ind1,ind2))*(MAX((sth-senv),0.)**2)**0.5
     2072            ! for mixed phase clouds, there is no contribution from large scale ratqs to the distribution
     2073            ! in the environement
     2074
     2075            sigma1s_ratqs=1E-10
     2076            IF (cloudth_ratqsmin>0.) THEN
     2077                sigma1s_ratqs = cloudth_ratqsmin*po(ind1)
     2078            ENDIF
     2079
     2080            sigma1s = sigma1s_fraca + sigma1s_ratqs
     2081            deltasenv=aenvl*vert_alpha*sigma1s
     2082            xenv1=-(senvl+deltasenv)/(sqrt(2.)*sigma1s)
     2083            xenv2=-(senvl-deltasenv)/(sqrt(2.)*sigma1s)
     2084            exp_xenv1 = exp(-1.*xenv1**2)
     2085            exp_xenv2 = exp(-1.*xenv2**2)
     2086
     2087            !surface CF
     2088            cenv(ind1,ind2)=0.5*(1.-1.*erf(xenv1))
     2089
     2090            !volume CF and condensed water
     2091            IntJ=0.5*senvl*(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp_xenv2
     2092            IntJ_CF=0.5*(1.-1.*erf(xenv2))
     2093
     2094            IF (deltasenv .LT. 1.e-10) THEN
     2095              qcenv(ind1,ind2)=IntJ
     2096              cenv_vol(ind1,ind2)=IntJ_CF
     2097            ELSE
     2098              IntI1=(((senvl+deltasenv)**2+(sigma1s)**2)/(8*deltasenv))*(erf(xenv2)-erf(xenv1))
     2099              IntI2=(sigma1s**2/(4*deltasenv*sqrtpi))*(xenv1*exp_xenv1-xenv2*exp_xenv2)
     2100              IntI3=((sqrt2*sigma1s*(senv+deltasenv))/(4*sqrtpi*deltasenv))*(exp_xenv1-exp_xenv2)
     2101              IntI1_CF=((senvl+deltasenv)*(erf(xenv2)-erf(xenv1)))/(4*deltasenv)
     2102              IntI3_CF=(sqrt2*sigma1s*(exp_xenv1-exp_xenv2))/(4*sqrtpi*deltasenv)
     2103              qcenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3 ! only liquid water in environment
     2104              cenv_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF
     2105            ENDIF
     2106
     2107            qcenv(ind1,ind2)=MAX(qcenv(ind1,ind2),0.)
     2108            cenv_vol(ind1,ind2)=MIN(MAX(cenv_vol(ind1,ind2),0.),1.)
     2109
     2110
     2111           
     2112            ! Thermals + environment
     2113            ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)
     2114            qctot(ind1,ind2)=fraca(ind1,ind2)*qcth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qcenv(ind1,ind2)
     2115            ctot_vol(ind1,ind2)=fraca(ind1,ind2)*cth_vol(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv_vol(ind1,ind2)
     2116            IF (qcth(ind1,ind2) .GT. 0) THEN
     2117                icefrac(ind1,ind2)=fraca(ind1,ind2)*qith(ind1,ind2)/(fraca(ind1,ind2)*qcth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qcenv(ind1,ind2))
     2118                icefrac(ind1,ind2)=MAX(MIN(1.,icefrac(ind1,ind2)),0.)
     2119            ELSE
     2120                icefrac(ind1,ind2)=0.
     2121            ENDIF
     2122
     2123            IF (cenv(ind1,ind2).LT.1.e-10.or.cth(ind1,ind2).LT.1.e-10) THEN
     2124                ctot(ind1,ind2)=0.
     2125                ctot_vol(ind1,ind2)=0.
     2126                qincloud(ind1)=0.
     2127                qcloud(ind1)=zqsatenv(ind1,ind2)
     2128            ELSE               
     2129                qcloud(ind1)=fraca(ind1,ind2)*(qcth(ind1,ind2)/cth(ind1,ind2)+qbase) &
     2130                            +(1.-1.*fraca(ind1,ind2))*(qcenv(ind1,ind2)/cenv(ind1,ind2)+qslenv)
     2131                qincloud(ind1)=MAX(fraca(ind1,ind2)*(qcth(ind1,ind2)/cth(ind1,ind2)) &
     2132                            +(1.-1.*fraca(ind1,ind2))*(qcenv(ind1,ind2)/cenv(ind1,ind2)),0.)
     2133            ENDIF
     2134
     2135        ENDIF ! mpc_bl_points
     2136
     2137
     2138    ELSE  ! gaussian for environment only
     2139
     2140     
     2141        zqenv(ind1)=po(ind1)
     2142        Tbef=temp(ind1,ind2)
     2143
     2144        CALL CALC_QSAT_ECMWF(Tbef,0.,paprs(ind1,ind2),RTT,0,.false.,qsatbef,dqsatdt)
     2145        zqsatenv(ind1,ind2)=qsatbef
     2146
     2147        IF (Tbef .GE. RTT) THEN
     2148                Lv=RLVTT
     2149        ELSE
     2150                Lv=RLSTT
     2151        ENDIF
     2152       
     2153
     2154        zthl(ind1,ind2)=temp(ind1,ind2)*(101325./paprs(ind1,ind2))**(rdd/cppd)
     2155        alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)
     2156        aenv=1./(1.+(alenv*Lv/cppd))
     2157        senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
     2158        sth=0.
     2159     
     2160        sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
     2161        sigma2s=0.
     2162
     2163        sqrt2pi=sqrt(2.*pi)
     2164        xenv=senv/(sqrt(2.)*sigma1s)
     2165        ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv))
     2166        ctot_vol(ind1,ind2)=ctot(ind1,ind2)
     2167        qctot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))
     2168     
     2169        IF (ctot(ind1,ind2).LT.1.e-3) THEN
     2170          ctot(ind1,ind2)=0.
     2171          qcloud(ind1)=zqsatenv(ind1,ind2)
     2172          qincloud(ind1)=0.
     2173        ELSE
     2174          qcloud(ind1)=qctot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
     2175          qincloud(ind1)=MAX(qctot(ind1,ind2)/ctot(ind1,ind2),0.)
     2176        ENDIF
     2177 
     2178
     2179    ENDIF       ! From the separation (thermal/envrionnement) and (environnement only,) l.335 et l.492
     2180
     2181    ! Outputs used to check the PDFs
     2182    cloudth_senv(ind1,ind2) = senv
     2183    cloudth_sth(ind1,ind2) = sth
     2184    cloudth_sigmaenv(ind1,ind2) = sigma1s
     2185    cloudth_sigmath(ind1,ind2) = sigma2s
     2186
     2187
     2188    ENDDO       !loop on klon
     2189
     2190RETURN
     2191
     2192
     2193END SUBROUTINE cloudth_mpc
     2194
     2195!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     2196SUBROUTINE ICE_MPC_BL_CLOUDS(ind1,ind2,klev,flag_topthermals,temp,pres,qth,qlth,qith,deltazlev,snowf,fraca,qi)
     2197
     2198! parameterization of ice for boundary
     2199! layer mixed-phase clouds assuming a stationary system
     2200!
     2201! Note that vapor deposition on ice crystals and riming of liquid droplets
     2202! depend on the ice number concentration Ni
     2203! One could assume that Ni depends on qi, e.g.,  Ni=beta*(qi*rho)**xi
     2204! and use values from Hong et al. 2004, MWR for instance
     2205! One may also estimate Ni as a function of T, as in Meyers 1922 or Fletcher 1962
     2206! One could also think of a more complex expression of Ni;
     2207! function of qi, T, the concentration in aerosols or INP ..
     2208! Here we prefer fixing Ni to a tuning parameter
     2209! By default we take 2.0L-1=2.0e3m-3, median value from measured vertical profiles near Svalbard
     2210! in Mioche et al. 2017
     2211!
     2212!
     2213! References:
     2214!------------
     2215! This parameterization is thoroughly described in Vignon et al.
     2216!
     2217! More specifically
     2218! for the Water vapor deposition process:
     2219!
     2220! Rotstayn, L. D., 1997: A physically based scheme for the treat-
     2221! ment of stratiform cloudfs and precipitation in large-scale
     2222! models. I: Description and evaluation of the microphysical
     2223! processes. Quart. J. Roy. Meteor. Soc., 123, 1227–1282.
     2224!
     2225!  Morrison, H., and A. Gettelman, 2008: A new two-moment bulk
     2226!  stratiform cloud microphysics scheme in the NCAR Com-
     2227!  munity Atmosphere Model (CAM3). Part I: Description and
     2228!  numerical tests. J. Climate, 21, 3642–3659
     2229!
     2230! for the Riming process:
     2231!
     2232! Rutledge, S. A., and P. V. Hobbs, 1983: The mesoscale and micro-
     2233! scale structure and organization of clouds and precipitation in
     2234! midlatitude cyclones. VII: A model for the ‘‘seeder-feeder’’
     2235! process in warm-frontal rainbands. J. Atmos. Sci., 40, 1185–1206
     2236!
     2237! Thompson, G., R. M. Rasmussen, and K. Manning, 004: Explicit
     2238! forecasts of winter precipitation using an improved bulk
     2239! microphysics scheme. Part I: Description and sensitivityThompson, G., R. M. Rasmussen, and K. Manning, 2004: Explicit
     2240! forecasts of winter precipitation using an improved bulk
     2241! microphysics scheme. Part I: Description and sensitivity analysis. Mon. Wea. Rev., 132, 519–542
     2242!
     2243! For the formation of clouds by thermals:
     2244!
     2245! Rio, C., & Hourdin, F. (2008). A thermal plume model for the convective boundary layer : Representation of cumulus clouds. Journal of
     2246! the Atmospheric Sciences, 65, 407–425.
     2247!
     2248! Jam, A., Hourdin, F., Rio, C., & Couvreux, F. (2013). Resolved versus parametrized boundary-layer plumes. Part III: Derivation of a
     2249! statistical scheme for cumulus clouds. Boundary-layer Meteorology, 147, 421–441. https://doi.org/10.1007/s10546-012-9789-3
     2250!
     2251!
     2252!
     2253! Contact: Etienne Vignon, etienne.vignon@lmd.ipsl.fr
     2254!=============================================================================
     2255
     2256    USE lscp_tools_mod, ONLY: CALC_QSAT_ECMWF
     2257    USE ioipsl_getin_p_mod, ONLY : getin_p
     2258    USE phys_state_var_mod, ONLY : fm_therm, detr_therm, entr_therm
     2259
     2260    IMPLICIT none
     2261
     2262
     2263    INCLUDE "YOMCST.h"
     2264    INCLUDE "nuage.h"
     2265
     2266    INTEGER, INTENT(IN) :: ind1,ind2, klev ! horizontal and vertical indices and dimensions
     2267    LOGICAL, INTENT(IN) :: flag_topthermals ! uppermost layer of thermals ?
     2268    REAL,  DIMENSION(klev), INTENT(IN) :: temp       ! temperature [K] within thermals
     2269    REAL,  DIMENSION(klev), INTENT(IN) :: pres       ! pressure [Pa]
     2270    REAL,  DIMENSION(klev), INTENT(IN) :: qth        ! mean specific water content in thermals [kg/kg]
     2271    REAL,  DIMENSION(klev), INTENT(IN) :: qlth       ! condensed liquid water in thermals, approximated value [kg/kg]
     2272    REAL,  DIMENSION(klev), INTENT(IN) :: qith       ! condensed ice water , thermals [kg/kg]
     2273    REAL,  DIMENSION(klev), INTENT(IN) :: deltazlev  ! layer thickness [m]
     2274    REAL,  DIMENSION(klev+1), INTENT(IN) :: snowf      ! snow flux at the upper inferface
     2275    REAL,  DIMENSION(klev+1), INTENT(IN) :: fraca      ! fraction of the mesh covered by thermals
     2276
     2277    REAL,  INTENT(OUT) :: qi        ! ice cloud specific content [kg/kg]
     2278
     2279
     2280    REAL, SAVE    :: Ni,  C_cap, Ei, d_top
     2281    !$OMP THREADPRIVATE(Ni, C_cap,Ei, d_top)
     2282    LOGICAL, SAVE :: firstcall = .TRUE.
     2283    !$OMP THREADPRIVATE(firstcall)
     2284
     2285
     2286    INTEGER ind2p1,ind2p2
     2287    REAL rho(klev)
     2288    REAL unsurtaudet, unsurtaustardep, unsurtaurim
     2289    REAL qsl, qsi, dqs, AA, BB, Ka, Dv, rhoi
     2290    REAL  p0, t0
     2291    REAL alpha, flux_term
     2292    REAL det_term, precip_term, rim_term, dep_term
     2293
     2294 
     2295    IF (firstcall) THEN
     2296        Ni=2.0e3
     2297        CALL getin_p('Ni', Ni)
     2298        WRITE(*,*) 'Ni = ', Ni
     2299
     2300        Ei=0.5
     2301        CALL getin_p('Ei', Ei)
     2302        WRITE(*,*) 'Ei = ', Ei
     2303
     2304        C_cap=0.5
     2305        CALL getin_p('C_cap', C_cap)
     2306        WRITE(*,*) 'C_cap = ', C_cap
     2307
     2308        d_top=0.8
     2309        CALL getin_p('d_top', d_top)
     2310        WRITE(*,*) 'd_top = ', d_top
     2311
     2312
     2313        firstcall=.FALSE.
     2314    ENDIF
     2315
     2316
     2317    ind2p1=ind2+1
     2318    ind2p2=ind2+2
     2319
     2320    ! Liquid water content:
     2321    !=====================
     2322    ! the liquid water content is not calculated in this routine
     2323
     2324    ! Ice water content
     2325    ! ==================
     2326
     2327    rho=pres/temp/RD  ! air density kg/m3
     2328
     2329    Ka=2.4e-2      ! thermal conductivity of the air, SI
     2330    p0=101325.0    ! ref pressure
     2331    T0=273.15      ! ref temp
     2332    rhoi=500.0     ! cloud ice density following Reisner et al. 1998
     2333    alpha=700.     ! fallvelocity param
     2334
     2335
     2336    IF (flag_topthermals) THEN ! uppermost thermals level, solve a third order polynomial with Cardan's method
     2337
     2338    Dv=0.0001*0.211*(p0/pres(ind2))*((temp(ind2)/T0)**1.94) ! water vapor diffusivity in air, SI
     2339
     2340    ! Detrainment term:
     2341    unsurtaudet=detr_therm(ind1,ind2)/rho(ind2)/deltazlev(ind2)
     2342 
     2343    ! vertical flux
     2344   
     2345    flux_term=d_top*fm_therm(ind1,ind2)/deltazlev(ind2)
     2346
     2347    ! Deposition term
     2348    CALL CALC_QSAT_ECMWF(temp(ind2),0.,pres(ind2),RTT,2,.false.,qsi,dqs)
     2349    CALL CALC_QSAT_ECMWF(temp(ind2),0.,pres(ind2),RTT,1,.false.,qsl,dqs)
     2350    AA=RLSTT/Ka/temp(ind2)*(RLSTT/RV/temp(ind2)-1.)
     2351    BB=1./(rho(ind2)*Dv*qsi)
     2352    unsurtaustardep=C_cap*(Ni**0.66)*(qsl-qsi)/qsi*4.*RPI/(AA+BB)*(6.*rho(ind2)/rhoi/RPI/Gamma(4.))**(0.33)
     2353
     2354    ! Riming term  neglected at this level
     2355    !unsurtaurim=rho(ind2)*alpha*3./rhoi/2.*Ei*qlth(ind2)*((p0/pres(ind2))**0.4)
     2356
     2357    qi=rho(ind2)*unsurtaustardep/MAX((rho(ind2)*unsurtaudet-flux_term),1E-12)
     2358    qi=MAX(qi,0.)**(3./2.)
     2359
     2360    ELSE ! other levels, estimate qi(k) from variables at k+1 and k+2
     2361
     2362    Dv=0.0001*0.211*(p0/pres(ind2p1))*((temp(ind2p1)/T0)**1.94) ! water vapor diffusivity in air, SI
     2363
     2364    ! Detrainment term:
     2365
     2366    unsurtaudet=detr_therm(ind1,ind2p1)/rho(ind2p1)/deltazlev(ind2p1)
     2367    det_term=-unsurtaudet*qith(ind2p1)*rho(ind2p1)
     2368   
     2369   
     2370    ! Deposition term
     2371
     2372    CALL CALC_QSAT_ECMWF(temp(ind2p1),0.,pres(ind2p1),RTT,2,.false.,qsi,dqs)
     2373    CALL CALC_QSAT_ECMWF(temp(ind2p1),0.,pres(ind2p1),RTT,1,.false.,qsl,dqs)
     2374    AA=RLSTT/Ka/temp(ind2p1)*(RLSTT/RV/temp(ind2p1)-1.)
     2375    BB=1./(rho(ind2p1)*Dv*qsi)
     2376    unsurtaustardep=C_cap*(Ni**0.66)*(qsl-qsi)/qsi*4.*RPI/(AA+BB)*(6.*rho(ind2p1)/rhoi/RPI/Gamma(4.))**(0.33)
     2377    dep_term=rho(ind2p1)*(qith(ind2p1)**0.33)*unsurtaustardep
     2378 
     2379    ! Riming term
     2380
     2381    unsurtaurim=rho(ind2p1)*alpha*3./rhoi/2.*Ei*qlth(ind2p1)*((p0/pres(ind2p1))**0.4)
     2382    rim_term=rho(ind2p1)*qith(ind2p1)*unsurtaurim
     2383
     2384    ! Precip term
     2385
     2386    !precip_term=-1./deltazlev(ind2p1)*(fraca(ind2p2)*snowf(ind2p2)-fraca(ind2p1)*snowf(ind2p1))
     2387    ! We assume that there is no solid precipitation outside thermals (so no multiplication by fraca)
     2388    precip_term=-1./deltazlev(ind2p1)*(snowf(ind2p2)-snowf(ind2p1))
     2389
     2390    ! Calculation in a top-to-bottom loop
     2391
     2392    IF (fm_therm(ind1,ind2p1) .GT. 0.) THEN
     2393       qi= 1./fm_therm(ind1,ind2p1)* &
     2394          (deltazlev(ind2p1)*(-rim_term-dep_term-det_term-precip_term) + &
     2395          fm_therm(ind1,ind2p2)*(qith(ind2p1)))
     2396    ELSE
     2397       qi=0.
     2398    ENDIF
     2399
     2400    ENDIF ! flag_topthermals
     2401 
     2402    qi=MAX(0.,qi)
     2403
     2404    RETURN
     2405
     2406END SUBROUTINE ICE_MPC_BL_CLOUDS
     2407
     2408
     2409
     2410
    15362411END MODULE cloudth_mod
    15372412

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

    r3989 r3999  
    182182    INTEGER,SAVE :: iflag_cloudth_vert_omp
    183183    INTEGER,SAVE :: iflag_rain_incloud_vol_omp
     184    INTEGER,SAVE :: iflag_vice_omp
    184185    REAL,SAVE :: rad_froid_omp, rad_chau1_omp, rad_chau2_omp
    185186    REAL,SAVE :: t_glace_min_omp, t_glace_max_omp
    186187    REAL,SAVE :: exposant_glace_omp
     188    INTEGER,SAVE :: iflag_gammasat_omp, iflag_mpc_bl_omp
    187189    REAL,SAVE :: rei_min_omp, rei_max_omp
    188190    INTEGER,SAVE :: iflag_sic_omp, iflag_inertie_omp
    … …  
    251253    LOGICAL, SAVE :: adjust_tropopause_omp
    252254    LOGICAL, SAVE :: ok_daily_climoz_omp
     255    LOGICAL, SAVE :: ok_new_lscp_omp
     256    LOGICAL, SAVE :: ok_icefra_lscp_omp
    253257
    254258    INTEGER, INTENT(OUT):: read_climoz ! read ozone climatology, OpenMP shared
    … …  
    14331437
    14341438    !
     1439    !Config Key  = iflag_gammasat
     1440    !Config Desc = 
     1441    !Config Def  = 0
     1442    !Config Help =
     1443    !
     1444    iflag_gammasat_omp=0
     1445    CALL getin('iflag_gammasat',iflag_gammasat_omp)
     1446
     1447
     1448    !
     1449    !Config Key  = iflag_mpc_bl
     1450    !Config Desc =
     1451    !Config Def  = 0
     1452    !Config Help =
     1453    !
     1454    iflag_mpc_bl_omp=0
     1455    CALL getin('iflag_mpc_bl',iflag_mpc_bl_omp)
     1456
     1457
     1458
     1459    !
    14351460    !Config Key  = iflag_t_glace
    14361461    !Config Desc = 
    … …  
    14581483    iflag_rain_incloud_vol_omp = 0
    14591484    CALL getin('iflag_rain_incloud_vol',iflag_rain_incloud_vol_omp)
     1485
     1486    !
     1487    !Config Key  = iflag_vice
     1488    !Config Desc = 
     1489    !Config Def  = 0
     1490    !Config Help =
     1491    !
     1492    iflag_vice_omp = 0
     1493    CALL getin('iflag_vice',iflag_vice_omp)
     1494
     1495
    14601496
    14611497    !
    … …  
    23132349    !Config Help = .FALSE. ensure much fewer (no calendar dependency)
    23142350    !  and lighter monthly climoz files, inetrpolated in time at gcm run time.
    2315     !
     2351   
     2352    ok_new_lscp_omp = .FALSE.
     2353    CALL getin('ok_new_lscp', ok_new_lscp_omp)
     2354    !
     2355    !Config Key  = ok_new_lscp_omp
     2356    !Config Desc = new cloud scheme ith ice and mixed phase (Etienne and JB)
     2357    !Config Def  = .FALSE.
     2358    !Config Help = ...
     2359
     2360
     2361
     2362    ok_icefra_lscp_omp = .FALSE.
     2363    CALL getin('ok_icefra_lscp', ok_icefra_lscp_omp)
     2364    !
     2365    !Config Key  = ok_icefra_lscp_omp
     2366    !Config Desc = ice fraction in radiation from lscp
     2367    !Config Def  = .FALSE.
     2368    !Config Help = ...
     2369
     2370
     2371
     2372
    23162373    ecrit_LES_omp = 1./8.
    23172374    CALL getin('ecrit_LES', ecrit_LES_omp)
    … …  
    24272484    t_glace_max = t_glace_max_omp
    24282485    exposant_glace = exposant_glace_omp
     2486    iflag_gammasat=iflag_gammasat_omp
     2487    iflag_mpc_bl=iflag_mpc_bl_omp
    24292488    iflag_t_glace = iflag_t_glace_omp
    24302489    iflag_cloudth_vert=iflag_cloudth_vert_omp
    24312490    iflag_rain_incloud_vol=iflag_rain_incloud_vol_omp
     2491    iflag_vice=iflag_vice_omp
    24322492    iflag_ice_thermo = iflag_ice_thermo_omp
    24332493    rei_min = rei_min_omp
    … …  
    26312691    carbon_cycle_rad = carbon_cycle_rad_omp
    26322692    level_coupling_esm = level_coupling_esm_omp
     2693    ok_new_lscp = ok_new_lscp_omp
     2694    ok_icefra_lscp=ok_icefra_lscp_omp
    26332695    read_fco2_ocean_cor = read_fco2_ocean_cor_omp
    26342696    var_fco2_ocean_cor = var_fco2_ocean_cor_omp
    … …  
    28492911    WRITE(lunout,*) ' t_glace_max = ',t_glace_max
    28502912    WRITE(lunout,*) ' exposant_glace = ',exposant_glace
     2913    WRITE(lunout,*) ' iflag_gammasat = ',iflag_gammasat
     2914    WRITE(lunout,*) ' iflag_mpc_bl = ',iflag_mpc_bl
    28512915    WRITE(lunout,*) ' iflag_t_glace = ',iflag_t_glace
    28522916    WRITE(lunout,*) ' iflag_cloudth_vert = ',iflag_cloudth_vert
    28532917    WRITE(lunout,*) ' iflag_rain_incloud_vol = ',iflag_rain_incloud_vol
     2918    WRITE(lunout,*) ' iflag_vice = ',iflag_vice
    28542919    WRITE(lunout,*) ' iflag_ice_thermo = ',iflag_ice_thermo
    28552920    WRITE(lunout,*) ' rei_min = ',rei_min
    … …  
    29643029    WRITE(lunout,*) ' adjust_tropopause = ', adjust_tropopause
    29653030    WRITE(lunout,*) ' ok_daily_climoz = ',ok_daily_climoz
     3031    WRITE(lunout,*) ' ok_new_lscp = ', ok_new_lscp
     3032    WRITE(lunout,*) ' ok_icefra_lscp = ', ok_icefra_lscp
    29663033    WRITE(lunout,*) ' read_climoz = ', read_climoz
    29673034    WRITE(lunout,*) ' carbon_cycle_tr = ', carbon_cycle_tr

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

    r3281 r3999  
    11! $Id$
    22
    3 SUBROUTINE newmicro(flag_aerosol, ok_cdnc, bl95_b0, bl95_b1, paprs, pplay, t, pqlwp, pclc, &
     3SUBROUTINE newmicro(flag_aerosol, ok_cdnc, bl95_b0, bl95_b1, paprs, pplay, t, pqlwp, picefra, pclc, &
    44    pcltau, pclemi, pch, pcl, pcm, pct, pctlwp, xflwp, xfiwp, xflwc, xfiwc, &
    55    mass_solu_aero, mass_solu_aero_pi, pcldtaupi, re, fl, reliq, reice, &
    … …  
    99  USE phys_local_var_mod, ONLY: scdnc, cldncl, reffclwtop, lcc, reffclws, &
    1010      reffclwc, cldnvi, lcc3d, lcc3dcon, lcc3dstra, icc3dcon, icc3dstra,  &
    11       zfice, dNovrN
     11      zfice, dNovrN, ptconv
    1212  USE phys_state_var_mod, ONLY: rnebcon, clwcon
    1313  USE icefrac_lsc_mod ! computes ice fraction (JBM 3/14)
    1414  USE ioipsl_getin_p_mod, ONLY : getin_p
    1515  USE print_control_mod, ONLY: lunout
     16  USE lscp_tools_mod, only: icefrac_lscp
     17
    1618
    1719
    … …  
    3133  ! pqlwp---input-R-eau liquide nuageuse dans l'atmosphere dans la partie
    3234  ! nuageuse (kg/kg)
     35  ! picefra--input-R-fraction de glace dans les nuages
    3336  ! pclc----input-R-couverture nuageuse pour le rayonnement (0 a 1)
    3437  ! mass_solu_aero-----input-R-total mass concentration for all soluble
    … …  
    5861  include "radepsi.h"
    5962  include "radopt.h"
     63  include "clesphys.h"
    6064
    6165  ! choix de l'hypothese de recouvrement nuageuse via radopt.h (IM, 19.07.2016)
    … …  
    8185  REAL t(klon, klev)
    8286  REAL pclc(klon, klev)
    83   REAL pqlwp(klon, klev)
     87  REAL pqlwp(klon, klev), picefra(klon,klev)
    8488  REAL pcltau(klon, klev)
    8589  REAL pclemi(klon, klev)
    … …  
    148152  ! jq-end
    149153  ! IM cf. CR:parametres supplementaires
     154  REAL dzfice(klon,klev)
    150155  REAL zclear(klon)
    151156  REAL zcloud(klon)
    … …  
    229234  ELSE ! of IF (iflag_t_glace.EQ.0)
    230235    DO k = 1, klev
    231         CALL icefrac_lsc(klon,t(:,k),pplay(:,k)/paprs(:,1),zfice(:,k))
    232  
    233 
    234         ! JBM: icefrac_lsc is now contained icefrac_lsc_mod
     236
     237! JBM: icefrac_lsc is now contained icefrac_lsc_mod
    235238!       zfice(i, k) = icefrac_lsc(t(i,k), t_glace_min, &
    236239!                                 t_glace_max, exposant_glace)
    237       DO i = 1, klon
     240
     241      IF (ok_new_lscp) THEN
     242          CALL icefrac_lscp(klon,t(:,k),pplay(:,k)/paprs(:,1),zfice(:,k),dzfice(:,k))
     243      ELSE
     244          CALL icefrac_lsc(klon,t(:,k),pplay(:,k)/paprs(:,1),zfice(:,k))
     245      ENDIF
     246
     247      DO i = 1, klon
     248       
     249        IF ((.NOT. ptconv(i,k)) .AND. ok_new_lscp .AND. ok_icefra_lscp) THEN
     250        ! EV: take the ice fraction directly from the lscp code
     251        ! consistent only for non convective grid points
     252        ! critical for mixed phase clouds
     253            zfice(i,k)=picefra(i,k)
     254        ENDIF
     255
    238256        ! -layer calculation
    239257        rhodz(i, k) = (paprs(i,k)-paprs(i,k+1))/rg ! kg/m2

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

    r2346 r3999  
    11! $Id$
    22
    3 SUBROUTINE nuage(paprs, pplay, t, pqlwp, pclc, pcltau, pclemi, pch, pcl, pcm, &
     3SUBROUTINE nuage(paprs, pplay, t, pqlwp,picefra, pclc, pcltau, pclemi, pch, pcl, pcm, &
    44    pct, pctlwp, ok_aie, mass_solu_aero, mass_solu_aero_pi, bl95_b0, bl95_b1, &
    55    cldtaupi, re, fl)
    66  USE dimphy
     7  USE lscp_tools_mod, only: icefrac_lscp
    78  USE icefrac_lsc_mod ! computes ice fraction (JBM 3/14)
     9  USE phys_local_var_mod, ONLY: ptconv
    810  IMPLICIT NONE
    911  ! ======================================================================
    … …  
    1416  ! t-------input-R-temperature
    1517  ! pqlwp---input-R-eau liquide nuageuse dans l'atmosphere (kg/kg)
     18  ! picefra--inout-R-fraction de glace dans les nuages (-)
    1619  ! pclc----input-R-couverture nuageuse pour le rayonnement (0 a 1)
    1720  ! ok_aie--input-L-apply aerosol indirect effect or not
    … …  
    3639  include "YOMCST.h"
    3740  include "nuage.h" ! JBM 3/14
     41  include "clesphys.h"
    3842
    3943  REAL paprs(klon, klev+1), pplay(klon, klev)
    … …  
    4145
    4246  REAL pclc(klon, klev)
    43   REAL pqlwp(klon, klev)
     47  REAL pqlwp(klon, klev), picefra(klon,klev)
    4448  REAL pcltau(klon, klev), pclemi(klon, klev)
    4549
    … …  
    8993
    9094  REAL cldtaupi(klon, klev) ! pre-industrial cloud opt thickness for diag
     95  REAl dzfice(klon)
    9196  ! jq-end
    9297
    … …  
    106111!       zfice(i) = icefrac_lsc(t(i,k), t_glace_min, &
    107112!                           t_glace_max, exposant_glace)
    108          CALL icefrac_lsc(klon,t(:,k),pplay(:,k)/paprs(:,1),zfice(:))
     113        IF (ok_new_lscp) THEN
     114            CALL icefrac_lscp(klon,t(:,k),pplay(:,k)/paprs(:,1),zfice(:),dzfice(:))
     115        ELSE
     116            CALL icefrac_lsc(klon,t(:,k),pplay(:,k)/paprs(:,1),zfice(:))
     117
     118        ENDIF
     119
     120        IF ((.NOT. ptconv(i,k)) .AND. ok_new_lscp .AND. ok_icefra_lscp) THEN
     121        ! EV: take the ice fraction directly from the lscp code
     122        ! consistent only for non convective grid points
     123        ! critical for mixed phase clouds
     124            DO i=1,klon
     125            zfice(i)=picefra(i,k)
     126            ENDDO
     127        ENDIF
     128
     129
    109130     ENDIF
    110131

  • LMDZ6/trunk/libf/phylmd/nuage.h

    r2945 r3999  
    1111      INTEGER iflag_t_glace, iflag_cloudth_vert, iflag_cld_cv
    1212      INTEGER iflag_rain_incloud_vol
     13   
     14      INTEGER iflag_mpc_bl, iflag_gammasat, iflag_vice
     15      LOGICAL ok_icefra_lscp
    1316
    1417      common /nuagecom/ rad_froid,rad_chau1, rad_chau2,t_glace_max,     &
    … …  
    1720     &                  tmax_fonte_cv,                                  &
    1821     &                  iflag_t_glace,iflag_cloudth_vert,iflag_cld_cv,  &
    19      &                  iflag_rain_incloud_vol
     22     &                  iflag_rain_incloud_vol,                         &
     23     &                  ok_icefra_lscp,                                 &
     24     &                  iflag_mpc_bl, iflag_gammasat, iflag_vice   
    2025!$OMP THREADPRIVATE(/nuagecom/)

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

    r3956 r3999  
    434434      REAL,ALLOCATABLE,SAVE,DIMENSION(:,:) :: cldemi, cldfra, cldtau, fiwc, fl, re, flwc
    435435!$OMP THREADPRIVATE(cldemi, cldfra, cldtau, fiwc, fl, re, flwc)
     436      REAL,ALLOCATABLE,SAVE,DIMENSION(:,:) :: qlth, qith
     437!$OMP THREADPRIVATE(qlth, qith)
    436438      REAL,ALLOCATABLE,SAVE,DIMENSION(:,:) :: ref_liq, ref_ice, theta, zphi
    437439!$OMP THREADPRIVATE(ref_liq, ref_ice, theta, zphi)
    … …  
    730732      ALLOCATE(rain_lsc(klon))
    731733      ALLOCATE(rain_num(klon))
    732 !
     734      ALLOCATE(qlth(klon,klev), qith(klon,klev))
     735      !
    733736      ALLOCATE(sens_x(klon), sens_w(klon))
    734737      ALLOCATE(zxfluxlat_x(klon), zxfluxlat_w(klon))
    … …  
    10291032      DEALLOCATE(rain_lsc)
    10301033      DEALLOCATE(rain_num)
     1034      DEALLOCATE(qlth, qith)
    10311035!
    10321036      DEALLOCATE(sens_x, sens_w)

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

    r3989 r3999  
    7474    USE VERTICAL_LAYERS_MOD, ONLY: aps,bps, ap, bp
    7575    USE write_field_phy
     76    USE lscp_mod, ONLY : lscp
    7677
    7778    !USE cmp_seri_mod
    … …  
    960961    !IM cf. AM 081204 BEG
    961962    LOGICAL ptconvth(klon,klev)
     963
     964    REAL picefra(klon,klev)
    962965    !IM cf. AM 081204 END
    963966    !
    … …  
    35103513       print *,'itap, ->fisrtilp ',itap
    35113514    ENDIF
    3512     !
     3515
     3516    picefra(:,:)=0.
     3517
     3518    IF (ok_new_lscp) THEN
     3519
     3520    CALL lscp(phys_tstep,paprs,pplay, &
     3521         t_seri, q_seri,ptconv,ratqs, &
     3522         d_t_lsc, d_q_lsc, d_ql_lsc, d_qi_lsc, rneb, cldliq, picefra, &
     3523         rain_lsc, snow_lsc, &
     3524         pfrac_impa, pfrac_nucl, pfrac_1nucl, &
     3525         frac_impa, frac_nucl, beta_prec_fisrt, &
     3526         prfl, psfl, rhcl,  &
     3527         zqasc, fraca,ztv,zpspsk,ztla,zthl,iflag_cld_th, &
     3528         iflag_ice_thermo)
     3529
     3530    ELSE
    35133531    CALL fisrtilp(phys_tstep,paprs,pplay, &
    35143532         t_seri, q_seri,ptconv,ratqs, &
    … …  
    35203538         zqasc, fraca,ztv,zpspsk,ztla,zthl,iflag_cld_th, &
    35213539         iflag_ice_thermo)
    3522     !
     3540    ENDIF
    35233541    WHERE (rain_lsc < 0) rain_lsc = 0.
    35243542    WHERE (snow_lsc < 0) snow_lsc = 0.
    … …  
    40304048          ENDIF
    40314049          CALL newmicro (flag_aerosol, ok_cdnc, bl95_b0, bl95_b1, &
    4032                paprs, pplay, t_seri, cldliq, cldfra, &
     4050               paprs, pplay, t_seri, cldliq, picefra, cldfra, &
    40334051               cldtau, cldemi, cldh, cldl, cldm, cldt, cldq, &
    40344052               flwp, fiwp, flwc, fiwc, &
    … …  
    40384056       ELSE
    40394057          CALL nuage (paprs, pplay, &
    4040                t_seri, cldliq, cldfra, cldtau, cldemi, &
     4058               t_seri, cldliq, picefra, cldfra, cldtau, cldemi, &
    40414059               cldh, cldl, cldm, cldt, cldq, &
    40424060               ok_aie, &
Note: See TracChangeset for help on using the changeset viewer.