Changeset 4368 for LMDZ6/branches/Ocean_skin/libf/phylmd/cloudth_mod.F90
- Timestamp:
- Dec 6, 2022, 12:01:16 AM (22 months ago)
- Location:
- LMDZ6/branches/Ocean_skin
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
libf/phylmd/cloudth_mod.F90 (modified) (47 diffs)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/branches/Ocean_skin
- Property svn:mergeinfo changed
-
LMDZ6/branches/Ocean_skin/libf/phylmd/cloudth_mod.F90
r4013 r4368 24 24 #include "YOETHF.h" 25 25 #include "FCTTRE.h" 26 #include "thermcell.h"27 26 #include "nuage.h" 28 27 … … 269 268 #include "YOETHF.h" 270 269 #include "FCTTRE.h" 271 #include "thermcell.h"272 270 #include "nuage.h" 273 271 … … 609 607 #include "YOETHF.h" 610 608 #include "FCTTRE.h" 611 #include "thermcell.h"612 609 #include "nuage.h" 613 610 … … 833 830 #include "YOETHF.h" 834 831 #include "FCTTRE.h" 835 #include "thermcell.h"836 832 #include "nuage.h" 837 833 … … 1150 1146 1151 1147 ELSE IF (iflag_cloudth_vert == 5) THEN 1152 sigma1s=(0.71794+0.000498239*dz(ind1,ind2))*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*(((sth-senv)**2)**0.5)+ratqs(ind1,ind2)*po(ind1) !Environment 1148 sigma1s=(0.71794+0.000498239*dz(ind1,ind2))*(fraca(ind1,ind2)**0.5) & 1149 /(1-fraca(ind1,ind2))*(((sth-senv)**2)**0.5) & 1150 +ratqs(ind1,ind2)*po(ind1) !Environment 1153 1151 sigma2s=(0.03218+0.000092655*dz(ind1,ind2))/((fraca(ind1,ind2)+0.02)**0.5)*(((sth-senv)**2)**0.5)+0.002*zqta(ind1,ind2) !Thermals 1154 1152 !sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1) … … 1295 1293 #include "YOETHF.h" 1296 1294 #include "FCTTRE.h" 1297 #include "thermcell.h"1298 1295 #include "nuage.h" 1299 1296 … … 1419 1416 !-------------------------------------------- 1420 1417 sigma_th=(0.03218+0.000092655*dz(ind1,ind2))/((fraca(ind1,ind2)+0.01)**0.5)*(((sth-senv)**2)**0.5)+0.002*zqta(ind1,ind2) 1421 sigma_env=(0.71794+0.000498239*dz(ind1,ind2))*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*(((sth-senv)**2)**0.5)+ratqs(ind1,ind2)*po(ind1) 1418 sigma_env=(0.71794+0.000498239*dz(ind1,ind2))*(fraca(ind1,ind2)**0.5) & 1419 /(1-fraca(ind1,ind2))*(((sth-senv)**2)**0.5) & 1420 +ratqs(ind1,ind2)*po(ind1) 1422 1421 xth=sth/(sqrt2*sigma_th) 1423 1422 xenv=senv/(sqrt2*sigma_env) … … 1539 1538 1540 1539 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1541 SUBROUTINE cloudth_mpc(klon,klev,ind2,mpc_bl_points, &1540 SUBROUTINE cloudth_mpc(klon,klev,ind2,mpc_bl_points, & 1542 1541 & temp,ztv,po,zqta,fraca,zpspsk,paprs,pplay,ztla,zthl, & 1543 1542 & ratqs,zqs,snowflux,qcloud,qincloud,icefrac,ctot,ctot_vol) … … 1554 1553 USE ioipsl_getin_p_mod, ONLY : getin_p 1555 1554 USE phys_output_var_mod, ONLY : cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv 1556 USE lscp_tools_mod, ONLY: CALC_QSAT_ECMWF, ICEFRAC_LSCP1557 USE phys_local_var_mod, ONLY : qlth, qith 1555 USE lscp_tools_mod, ONLY: CALC_QSAT_ECMWF, FALLICE_VELOCITY 1556 USE phys_local_var_mod, ONLY : qlth, qith, qsith, wiceth 1558 1557 1559 1558 IMPLICIT NONE … … 1562 1561 #include "YOETHF.h" 1563 1562 #include "FCTTRE.h" 1564 #include "thermcell.h"1565 1563 #include "nuage.h" 1566 1564 … … 1573 1571 1574 1572 INTEGER, INTENT(IN) :: klon,klev,ind2 1575 INTEGER, DIMENSION(klon,klev), INTENT(INOUT) :: mpc_bl_points ! 1 where BL MPC, 0 otherwise1573 1576 1574 1577 1575 REAL, DIMENSION(klon,klev), INTENT(IN) :: temp ! Temperature [K] … … 1589 1587 REAL, DIMENSION(klon,klev+1), INTENT(IN) :: snowflux ! snow flux at the interface of the layer [kg/m2/s] 1590 1588 1591 1589 INTEGER, DIMENSION(klon,klev), INTENT(INOUT) :: mpc_bl_points ! grid points with convective (thermals) mixed phase clouds 1590 1592 1591 REAL, DIMENSION(klon,klev), INTENT(OUT) :: ctot ! Cloud fraction [0-1] 1593 1592 REAL, DIMENSION(klon,klev), INTENT(OUT) :: ctot_vol ! Volume cloud fraction [0-1] … … 1600 1599 1601 1600 INTEGER itap,ind1,l,ig,iter,k 1602 LOGICAL flag_topthermals1603 1604 1605 REAL zqsatth(klon ,klev), zqsatenv(klon,klev)1601 INTEGER iflag_topthermals, niter 1602 LOGICAL falseklon(klon) 1603 1604 REAL zqsatth(klon), zqsatenv(klon), zqsatenvonly(klon) 1606 1605 REAL sigma1(klon,klev) 1607 1606 REAL sigma2(klon,klev) … … 1614 1613 REAL cenv_vol(klon,klev) 1615 1614 REAL rneb(klon,klev) 1616 REAL zqenv(klon) 1617 1615 REAL zqenv(klon), zqth(klon), zqenvonly(klon) 1618 1616 REAL qsatmmussig1,qsatmmussig2,sqrtpi,sqrt2,sqrt2pi,pi 1619 1617 REAL rdd,cppd,Lv … … 1624 1622 REAL xth1,xth2,xenv1,xenv2,deltasth, deltasenv 1625 1623 REAL IntJ,IntI1,IntI2,IntI3,IntJ_CF,IntI1_CF,IntI3_CF,coeffqlenv,coeffqlth 1626 REAL Tbef,zdelta,qsatbef,zcor 1627 REAL qlbef,dqsatdt 1624 REAL zdelta,qsatbef,zcor 1625 REAL Tbefth(klon), Tbefenv(klon), Tbefenvonly(klon), rhoth(klon) 1626 REAL qlbef 1627 REAL dqsatenv(klon), dqsatth(klon), dqsatenvonly(klon) 1628 1628 REAL erf 1629 1629 REAL zpdf_sig(klon),zpdf_k(klon),zpdf_delta(klon) … … 1632 1632 REAL zrho(klon,klev) 1633 1633 REAL dz(klon,klev) 1634 REAL qslth, qsith, qslenv, alenvl, aenvl 1635 REAL sthi, sthl, althl, athl 1634 REAL qslenv(klon), qslth(klon) 1635 REAL alenvl, aenvl 1636 REAL sthi, sthl, sthil, althl, athl, althi, athi, sthlc, deltasthc, sigma2sc 1636 1637 REAL senvi, senvl, qbase, sbase, qliqth, qiceth 1637 REAL q imax, ttarget, stmp, cout, coutref1638 REAL qmax, ttarget, stmp, cout, coutref, fraci 1638 1639 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 1640 REAL deltazlev_mpc(klev), temp_mpc(klev), pres_mpc(klev), fraca_mpc(klev+1), snowf_mpc(klev+1) 1642 1641 1643 1642 ! Modifty the saturation deficit PDF in thermals … … 1645 1644 REAL,SAVE :: C_mpc 1646 1645 !$OMP THREADPRIVATE(C_mpc) 1646 REAL, SAVE :: Ni,C_cap,Ei,d_top 1647 !$OMP THREADPRIVATE(Ni,C_cap,Ei,d_top) 1647 1648 ! Change the width of the PDF used for vertical subgrid scale heterogeneity 1648 1649 ! (J Jouhaud, JL Dufresne, JB Madeleine) … … 1692 1693 qlth(:,ind2)=0. 1693 1694 qith(:,ind2)=0. 1695 wiceth(:,ind2)=0. 1694 1696 rneb(:,:)=0. 1695 1697 qcloud(:)=0. … … 1700 1702 cenv_vol(:,:)=0. 1701 1703 ctot_vol(:,:)=0. 1704 falseklon(:)=.false. 1702 1705 qsatmmussig1=0. 1703 1706 qsatmmussig2=0. … … 1709 1712 sqrtpi=sqrt(pi) 1710 1713 icefrac(:,ind2)=0. 1714 althl=0. 1715 althi=0. 1716 athl=0. 1717 athi=0. 1718 senvl=0. 1719 senvi=0. 1720 sthi=0. 1721 sthl=0. 1722 sthil=0. 1711 1723 1712 1724 … … 1741 1753 WRITE(*,*) 'iflag_cloudth_vert_noratqs = ', iflag_cloudth_vert_noratqs 1742 1754 ! Modifies the PDF in thermals when ice crystals are present 1743 C_mpc=1.e 21755 C_mpc=1.e4 1744 1756 CALL getin_p('C_mpc',C_mpc) 1745 1757 WRITE(*,*) 'C_mpc = ', C_mpc 1758 Ni=2.0e3 1759 CALL getin_p('Ni', Ni) 1760 WRITE(*,*) 'Ni = ', Ni 1761 Ei=0.5 1762 CALL getin_p('Ei', Ei) 1763 WRITE(*,*) 'Ei = ', Ei 1764 C_cap=0.5 1765 CALL getin_p('C_cap', C_cap) 1766 WRITE(*,*) 'C_cap = ', C_cap 1767 d_top=1.2 1768 CALL getin_p('d_top', d_top) 1769 WRITE(*,*) 'd_top = ', d_top 1746 1770 1747 1771 firstcall=.FALSE. … … 1759 1783 DO ind1=1,klon 1760 1784 IF ((temp(ind1,ind2) .LT. RTT) .AND. (temp(ind1,ind2) .GT. temp_lim) & 1761 .AND. (iflag_mpc_bl .GE. 2) .AND. (ind2<=klev-2) &1785 .AND. (iflag_mpc_bl .GE. 1) .AND. (ind2<=klev-2) & 1762 1786 .AND. (ztv(ind1,1).GT.ztv(ind1,2)) .AND.(fraca(ind1,ind2).GT.1.e-10)) THEN 1763 1787 mpc_bl_points(ind1,ind2)=1 … … 1772 1796 !------------------------------------------------------------------------------- 1773 1797 1774 DO ind1=1,klon 1798 ! calculation of temperature, humidity and saturation specific humidity 1799 1800 Tbefenv(:)=zthl(:,ind2)*zpspsk(:,ind2) 1801 Tbefth(:)=ztla(:,ind2)*zpspsk(:,ind2) 1802 Tbefenvonly(:)=temp(:,ind2) 1803 rhoth(:)=paprs(:,ind2)/Tbefth(:)/RD 1804 1805 zqenv(:)=(po(:)-fraca(:,ind2)*zqta(:,ind2))/(1.-fraca(:,ind2)) !qt = a*qtth + (1-a)*qtenv 1806 zqth(:)=zqta(:,ind2) 1807 zqenvonly(:)=po(:) 1808 1809 1810 CALL CALC_QSAT_ECMWF(klon,Tbefenvonly,zqenvonly,paprs(:,ind2),RTT,0,.false.,zqsatenvonly,dqsatenv) 1811 1812 CALL CALC_QSAT_ECMWF(klon,Tbefenv,zqenv,paprs(:,ind2),RTT,0,.false.,zqsatenv,dqsatenv) 1813 CALL CALC_QSAT_ECMWF(klon,Tbefenv,zqenv,paprs(:,ind2),RTT,1,.false.,qslenv,dqsatenv) 1814 1815 CALL CALC_QSAT_ECMWF(klon,Tbefth,zqth,paprs(:,ind2),RTT,1,.false.,qslth,dqsatth) 1816 CALL CALC_QSAT_ECMWF(klon,Tbefth,zqth,paprs(:,ind2),RTT,2,.false.,qsith(:,ind2),dqsatth) 1817 CALL CALC_QSAT_ECMWF(klon,Tbefth,zqth,paprs(:,ind2),RTT,0,.false.,zqsatth,dqsatth) 1818 1819 1820 DO ind1=1,klon 1775 1821 1776 1822 … … 1780 1826 ! Environment: 1781 1827 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 1828 1829 IF (Tbefenv(ind1) .GE. RTT) THEN 1789 1830 Lv=RLVTT 1790 1831 ELSE … … 1793 1834 1794 1835 1795 alenv=(0.622*Lv*zqsatenv(ind1 ,ind2))/(rdd*zthl(ind1,ind2)**2) !qsl, p841836 alenv=(0.622*Lv*zqsatenv(ind1))/(rdd*zthl(ind1,ind2)**2) !qsl, p84 1796 1837 aenv=1./(1.+(alenv*Lv/cppd)) !al, p84 1797 senv=aenv*(po(ind1)-zqsatenv(ind1 ,ind2)) !s, p841838 senv=aenv*(po(ind1)-zqsatenv(ind1)) !s, p84 1798 1839 1799 1840 ! For MPCs: 1800 1841 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) 1842 alenvl=(0.622*RLVTT*qslenv(ind1))/(rdd*zthl(ind1,ind2)**2) 1803 1843 aenvl=1./(1.+(alenv*Lv/cppd)) 1804 senvl=aenvl*(po(ind1)-qslenv) 1844 senvl=aenvl*(po(ind1)-qslenv(ind1)) 1845 senvi=senv 1805 1846 ENDIF 1806 1847 … … 1808 1849 ! Thermals: 1809 1850 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 1851 1852 IF (Tbefth(ind1) .GE. RTT) THEN 1815 1853 Lv=RLVTT 1816 1854 ELSE … … 1819 1857 1820 1858 1821 alth=(0.622*Lv*zqsatth(ind1 ,ind2))/(rdd*ztla(ind1,ind2)**2)1859 alth=(0.622*Lv*zqsatth(ind1))/(rdd*ztla(ind1,ind2)**2) 1822 1860 ath=1./(1.+(alth*Lv/cppd)) 1823 sth=ath*(zqta(ind1,ind2)-zqsatth(ind1 ,ind2))1861 sth=ath*(zqta(ind1,ind2)-zqsatth(ind1)) 1824 1862 1825 1863 ! For MPCs: 1826 1864 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) 1865 althi=alth 1866 althl=(0.622*RLVTT*qslth(ind1))/(rdd*ztla(ind1,ind2)**2) 1867 athl=1./(1.+(alth*RLVTT/cppd)) 1868 athi=alth 1869 sthl=athl*(zqta(ind1,ind2)-qslth(ind1)) 1870 sthi=athi*(zqta(ind1,ind2)-qsith(ind1,ind2)) 1871 sthil=athi*(zqta(ind1,ind2)-qslth(ind1)) 1833 1872 ENDIF 1834 1873 … … 1923 1962 ctot(ind1,ind2)=0. 1924 1963 ctot_vol(ind1,ind2)=0. 1925 qcloud(ind1)=zqsatenv(ind1 ,ind2)1964 qcloud(ind1)=zqsatenv(ind1) 1926 1965 qincloud(ind1)=0. 1927 1966 ELSE … … 1941 1980 !........... 1942 1981 1982 qiceth=0. 1943 1983 deltazlev_mpc=dz(ind1,:) 1944 1984 temp_mpc=ztla(ind1,:)*zpspsk(ind1,:) … … 1946 1986 fraca_mpc=fraca(ind1,:) 1947 1987 snowf_mpc=snowflux(ind1,:) 1948 qth_mpc=zqta(ind1,:) 1949 flag_topthermals=.FALSE. 1988 iflag_topthermals=0 1950 1989 IF ((mpc_bl_points(ind1,ind2) .EQ. 1) .AND. (mpc_bl_points(ind1,ind2+1) .EQ. 0)) THEN 1951 flag_topthermals = .TRUE. 1990 iflag_topthermals = 1 1991 ELSE IF ((mpc_bl_points(ind1,ind2) .EQ. 1) .AND. (mpc_bl_points(ind1,ind2+1) .EQ. 1) & 1992 .AND. (mpc_bl_points(ind1,ind2+2) .EQ. 0) ) THEN 1993 iflag_topthermals = 2 1994 ELSE 1995 iflag_topthermals = 0 1952 1996 ENDIF 1953 1997 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) 1998 CALL ICE_MPC_BL_CLOUDS(ind1,ind2,klev,Ni,Ei,C_cap,d_top,iflag_topthermals,temp_mpc,pres_mpc,zqta(ind1,:), & 1999 qsith(ind1,:),qlth(ind1,:),deltazlev_mpc,wiceth(ind1,:),fraca_mpc,qith(ind1,:)) 2000 2001 ! qmax calculation 2002 sigma2s=(sigma2s_factor*((MAX((sthl-senvl),0.)**2)**0.5)/((fraca(ind1,ind2)+0.02)**sigma2s_power))+0.002*zqta(ind1,ind2) 1971 2003 deltasth=athl*vert_alpha_th*sigma2s 1972 1973 ! Liquid phase1974 !.............1975 2004 xth1=-(sthl+deltasth)/(sqrt(2.)*sigma2s) 1976 xth2=-(sthl-deltasth)/(sqrt(2.)*sigma2s) 2005 xth2=-(sthl-deltasth)/(sqrt(2.)*sigma2s) 1977 2006 exp_xth1 = exp(-1.*xth1**2) 1978 2007 exp_xth2 = exp(-1.*xth2**2) … … 1984 2013 IntI1_CF=((sthl+deltasth)*(erf(xth2)-erf(xth1)))/(4*deltasth) 1985 2014 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) 2015 qmax=MAX(IntJ+IntI1+IntI2+IntI3,0.) 2016 2017 2018 ! Liquid phase 2019 !................ 2020 ! We account for the effect of ice crystals in thermals on sthl 2021 ! and on the width of the distribution 2022 2023 sthlc=sthl*1./(1.+C_mpc*qith(ind1,ind2)) & 2024 + (1.-1./(1.+C_mpc*qith(ind1,ind2))) * athl*(qsith(ind1,ind2)-qslth(ind1)) 2025 2026 sigma2sc=(sigma2s_factor*((MAX((sthlc-senvl),0.)**2)**0.5) & 2027 /((fraca(ind1,ind2)+0.02)**sigma2s_power)) & 2028 +0.002*zqta(ind1,ind2) 2029 deltasthc=athl*vert_alpha_th*sigma2sc 2030 2031 2032 xth1=-(sthlc+deltasthc)/(sqrt(2.)*sigma2sc) 2033 xth2=-(sthlc-deltasthc)/(sqrt(2.)*sigma2sc) 1991 2034 exp_xth1 = exp(-1.*xth1**2) 1992 2035 exp_xth2 = exp(-1.*xth2**2) 1993 IntJ=0.5*sth i*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth22036 IntJ=0.5*sthlc*(1-erf(xth2))+(sigma2sc/sqrt2pi)*exp_xth2 1994 2037 IntJ_CF=0.5*(1.-1.*erf(xth2)) 1995 IntI1=(((sth i+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=((sth i+deltasth)*(erf(xth2)-erf(xth1)))/(4*deltasth)1999 IntI3_CF=(sqrt2*sigma2s *(exp_xth1-exp_xth2))/(4*sqrtpi*deltasth)2000 q imax=IntJ+IntI1+IntI2+IntI32001 qimax=qimax-qliqth2002 2038 IntI1=(((sthlc+deltasthc)**2+(sigma2sc)**2)/(8*deltasthc))*(erf(xth2)-erf(xth1)) 2039 IntI2=(sigma2sc**2/(4*deltasthc*sqrtpi))*(xth1*exp_xth1-xth2*exp_xth2) 2040 IntI3=((sqrt2*sigma2sc*(sthlc+deltasthc))/(4*sqrtpi*deltasthc))*(exp_xth1-exp_xth2) 2041 IntI1_CF=((sthlc+deltasthc)*(erf(xth2)-erf(xth1)))/(4*deltasthc) 2042 IntI3_CF=(sqrt2*sigma2sc*(exp_xth1-exp_xth2))/(4*sqrtpi*deltasthc) 2043 qliqth=IntJ+IntI1+IntI2+IntI3 2044 2045 qlth(ind1,ind2)=MAX(0.,qliqth) 2003 2046 2004 2047 ! 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)) 2048 2008 2049 qcth(ind1,ind2)=qlth(ind1,ind2)+qith(ind1,ind2) 2009 2050 2051 2052 ! consistency with subgrid distribution 2053 2054 IF ((qcth(ind1,ind2) .GT. qmax) .AND. (qcth(ind1,ind2) .GT. 0)) THEN 2055 fraci=qith(ind1,ind2)/qcth(ind1,ind2) 2056 qcth(ind1,ind2)=qmax 2057 qith(ind1,ind2)=fraci*qmax 2058 qlth(ind1,ind2)=(1.-fraci)*qmax 2059 ENDIF 2060 2061 ! Cloud Fraction 2062 !............... 2010 2063 ! calculation of qbase which is the value of the water vapor within mixed phase clouds 2011 2064 ! such that the total water in cloud = qbase+qliqth+qiceth … … 2015 2068 2016 2069 ttarget=qcth(ind1,ind2) 2017 mini=athl*(qsith-qslth) 2018 maxi=0. 2070 mini= athl*(qsith(ind1,ind2)-qslth(ind1)) 2071 maxi= 0. !athl*(3.*sqrt(sigma2s)) 2072 niter=20 2019 2073 pas=(maxi-mini)/niter 2020 2074 stmp=mini … … 2031 2085 ENDIF 2032 2086 ENDDO 2033 qbase=MAX(0., sbase/athl+qslth )2087 qbase=MAX(0., sbase/athl+qslth(ind1)) 2034 2088 2035 2089 ! surface cloud fraction in thermals … … 2053 2107 IntI1=(((sthl+deltasth-sbase)**2+(sigma2s)**2)/(8*deltasth))*(erf(xth2)-erf(xth1)) 2054 2108 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)2109 IntI3=((sqrt2*sigma2s*(sthl+deltasth))/(4*sqrtpi*deltasth))*(exp_xth1-exp_xth2) 2056 2110 IntI1_CF=((sthl-sbase+deltasth)*(erf(xth2)-erf(xth1)))/(4*deltasth) 2057 2111 IntI3_CF=(sqrt2*sigma2s*(exp_xth1-exp_xth2))/(4*sqrtpi*deltasth) … … 2059 2113 ENDIF 2060 2114 cth_vol(ind1,ind2)=MIN(MAX(0.,cth_vol(ind1,ind2)),1.) 2115 2116 2061 2117 2062 2118 ! Environment … … 2110 2166 2111 2167 2112 ! Thermals + environment 2168 ! Thermals + environment 2169 ! ======================= 2113 2170 ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2) 2114 2171 qctot(ind1,ind2)=fraca(ind1,ind2)*qcth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qcenv(ind1,ind2) 2115 2172 ctot_vol(ind1,ind2)=fraca(ind1,ind2)*cth_vol(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv_vol(ind1,ind2) 2116 2173 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)) 2174 icefrac(ind1,ind2)=fraca(ind1,ind2)*qith(ind1,ind2) & 2175 /(fraca(ind1,ind2)*qcth(ind1,ind2) & 2176 +(1.-1.*fraca(ind1,ind2))*qcenv(ind1,ind2)) 2118 2177 icefrac(ind1,ind2)=MAX(MIN(1.,icefrac(ind1,ind2)),0.) 2119 2178 ELSE … … 2125 2184 ctot_vol(ind1,ind2)=0. 2126 2185 qincloud(ind1)=0. 2127 qcloud(ind1)=zqsatenv(ind1 ,ind2)2186 qcloud(ind1)=zqsatenv(ind1) 2128 2187 ELSE 2129 2188 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 )2189 +(1.-1.*fraca(ind1,ind2))*(qcenv(ind1,ind2)/cenv(ind1,ind2)+qslenv(ind1)) 2131 2190 qincloud(ind1)=MAX(fraca(ind1,ind2)*(qcth(ind1,ind2)/cth(ind1,ind2)) & 2132 2191 +(1.-1.*fraca(ind1,ind2))*(qcenv(ind1,ind2)/cenv(ind1,ind2)),0.) … … 2139 2198 2140 2199 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 2200 2201 2202 IF (Tbefenvonly(ind1) .GE. RTT) THEN 2148 2203 Lv=RLVTT 2149 2204 ELSE … … 2153 2208 2154 2209 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)2210 alenv=(0.622*Lv*zqsatenvonly(ind1))/(rdd*zthl(ind1,ind2)**2) 2156 2211 aenv=1./(1.+(alenv*Lv/cppd)) 2157 senv=aenv*(po(ind1)-zqsatenv (ind1,ind2))2212 senv=aenv*(po(ind1)-zqsatenvonly(ind1)) 2158 2213 sth=0. 2159 2214 2160 sigma1s=ratqs(ind1,ind2)*zqenv (ind1)2215 sigma1s=ratqs(ind1,ind2)*zqenvonly(ind1) 2161 2216 sigma2s=0. 2162 2217 … … 2169 2224 IF (ctot(ind1,ind2).LT.1.e-3) THEN 2170 2225 ctot(ind1,ind2)=0. 2171 qcloud(ind1)=zqsatenv (ind1,ind2)2226 qcloud(ind1)=zqsatenvonly(ind1) 2172 2227 qincloud(ind1)=0. 2173 2228 ELSE 2174 qcloud(ind1)=qctot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv (ind1,ind2)2229 qcloud(ind1)=qctot(ind1,ind2)/ctot(ind1,ind2)+zqsatenvonly(ind1) 2175 2230 qincloud(ind1)=MAX(qctot(ind1,ind2)/ctot(ind1,ind2),0.) 2176 2231 ENDIF … … 2186 2241 2187 2242 2188 ENDDO !loop on klon 2243 ENDDO !loop on klon 2244 2245 ! Calcule ice fall velocity in thermals 2246 2247 CALL FALLICE_VELOCITY(klon,qith(:,ind2),Tbefth(:),rhoth(:),paprs(:,ind2),falseklon(:),wiceth(:,ind2)) 2189 2248 2190 2249 RETURN … … 2194 2253 2195 2254 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 2196 SUBROUTINE ICE_MPC_BL_CLOUDS(ind1,ind2,klev, flag_topthermals,temp,pres,qth,qlth,qith,deltazlev,snowf,fraca,qi)2255 SUBROUTINE ICE_MPC_BL_CLOUDS(ind1,ind2,klev,Ni,Ei,C_cap,d_top,iflag_topthermals,temp,pres,qth,qsith,qlth,deltazlev,vith,fraca,qith) 2197 2256 2198 2257 ! parameterization of ice for boundary … … 2254 2313 !============================================================================= 2255 2314 2256 USE lscp_tools_mod, ONLY: CALC_QSAT_ECMWF2257 2315 USE ioipsl_getin_p_mod, ONLY : getin_p 2258 2316 USE phys_state_var_mod, ONLY : fm_therm, detr_therm, entr_therm … … 2260 2318 IMPLICIT none 2261 2319 2262 2263 2320 INCLUDE "YOMCST.h" 2264 2321 INCLUDE "nuage.h" 2265 2322 2266 INTEGER, INTENT(IN) :: ind1,ind2, klev ! horizontal and vertical indices and dimensions 2267 LOGICAL, INTENT(IN) :: flag_topthermals ! uppermost layer of thermals ? 2323 INTEGER, INTENT(IN) :: ind1,ind2, klev ! horizontal and vertical indices and dimensions 2324 INTEGER, INTENT(IN) :: iflag_topthermals ! uppermost layer of thermals ? 2325 REAL, INTENT(IN) :: Ni ! ice number concentration [m-3] 2326 REAL, INTENT(IN) :: Ei ! ice-droplet collision efficiency 2327 REAL, INTENT(IN) :: C_cap ! ice crystal capacitance 2328 REAL, INTENT(IN) :: d_top ! cloud-top ice crystal mixing parameter 2268 2329 REAL, DIMENSION(klev), INTENT(IN) :: temp ! temperature [K] within thermals 2269 2330 REAL, DIMENSION(klev), INTENT(IN) :: pres ! pressure [Pa] 2270 2331 REAL, DIMENSION(klev), INTENT(IN) :: qth ! mean specific water content in thermals [kg/kg] 2332 REAL, DIMENSION(klev), INTENT(IN) :: qsith ! saturation specific humidity wrt ice in thermals [kg/kg] 2271 2333 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 2334 REAL, DIMENSION(klev), INTENT(IN) :: deltazlev ! layer thickness [m] 2274 REAL, DIMENSION(klev +1), INTENT(IN) :: snowf ! snow flux at the upper inferface2335 REAL, DIMENSION(klev), INTENT(IN) :: vith ! ice crystal fall velocity [m/s] 2275 2336 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) 2337 REAL, DIMENSION(klev), INTENT(INOUT) :: qith ! condensed ice water , thermals [kg/kg] 2284 2338 2285 2339 … … 2287 2341 REAL rho(klev) 2288 2342 REAL unsurtaudet, unsurtaustardep, unsurtaurim 2289 REAL q sl, qsi, dqs, AA, BB, Ka, Dv, rhoi2290 REAL p0, t02343 REAL qi, AA, BB, Ka, Dv, rhoi 2344 REAL p0, t0, fp1, fp2 2291 2345 REAL alpha, flux_term 2292 2346 REAL det_term, precip_term, rim_term, dep_term 2293 2347 2294 2295 IF (firstcall) THEN2296 Ni=2.0e32297 CALL getin_p('Ni', Ni)2298 WRITE(*,*) 'Ni = ', Ni2299 2300 Ei=0.52301 CALL getin_p('Ei', Ei)2302 WRITE(*,*) 'Ei = ', Ei2303 2304 C_cap=0.52305 CALL getin_p('C_cap', C_cap)2306 WRITE(*,*) 'C_cap = ', C_cap2307 2308 d_top=0.82309 CALL getin_p('d_top', d_top)2310 WRITE(*,*) 'd_top = ', d_top2311 2312 2313 firstcall=.FALSE.2314 ENDIF2315 2316 2348 2317 2349 ind2p1=ind2+1 2318 2350 ind2p2=ind2+2 2319 2320 ! Liquid water content:2321 !=====================2322 ! the liquid water content is not calculated in this routine2323 2324 ! Ice water content2325 ! ==================2326 2351 2327 2352 rho=pres/temp/RD ! air density kg/m3 … … 2334 2359 2335 2360 2336 IF ( flag_topthermals) THEN ! uppermost thermals level, solve a third order polynomial with Cardan's method2361 IF (iflag_topthermals .GT. 0) THEN ! uppermost thermals levels 2337 2362 2338 2363 Dv=0.0001*0.211*(p0/pres(ind2))*((temp(ind2)/T0)**1.94) ! water vapor diffusivity in air, SI 2339 2364 2340 2365 ! Detrainment term: 2366 2341 2367 unsurtaudet=detr_therm(ind1,ind2)/rho(ind2)/deltazlev(ind2) 2342 2368 2343 ! vertical flux2344 2345 flux_term=d_top*fm_therm(ind1,ind2)/deltazlev(ind2)2346 2369 2347 2370 ! 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 2371 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) 2372 BB=1./(rho(ind2)*Dv*qsith(ind2)) 2373 unsurtaustardep=C_cap*(Ni**0.66)*(qth(ind2)-qsith(ind2))/qsith(ind2) & 2374 *4.*RPI/(AA+BB)*(6.*rho(ind2)/rhoi/RPI/Gamma(4.))**(0.33) 2353 2375 2354 2376 ! Riming term neglected at this level 2355 2377 !unsurtaurim=rho(ind2)*alpha*3./rhoi/2.*Ei*qlth(ind2)*((p0/pres(ind2))**0.4) 2356 2378 2357 qi= rho(ind2)*unsurtaustardep/MAX((rho(ind2)*unsurtaudet-flux_term),1E-12)2379 qi=fraca(ind2)*unsurtaustardep/MAX((d_top*unsurtaudet),1E-12) 2358 2380 qi=MAX(qi,0.)**(3./2.) 2359 2381 … … 2370 2392 ! Deposition term 2371 2393 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 2394 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 2395 BB=1./(rho(ind2p1)*Dv*qsith(ind2p1)) 2396 unsurtaustardep=C_cap*(Ni**0.66)*(qth(ind2p1)-qsith(ind2p1)) & 2397 /qsith(ind2p1)*4.*RPI/(AA+BB) & 2398 *(6.*rho(ind2p1)/rhoi/RPI/Gamma(4.))**(0.33) 2399 dep_term=rho(ind2p1)*fraca(ind2p1)*(qith(ind2p1)**0.33)*unsurtaustardep 2378 2400 2379 2401 ! Riming term 2380 2402 2381 2403 unsurtaurim=rho(ind2p1)*alpha*3./rhoi/2.*Ei*qlth(ind2p1)*((p0/pres(ind2p1))**0.4) 2382 rim_term=rho(ind2p1)* qith(ind2p1)*unsurtaurim2404 rim_term=rho(ind2p1)*fraca(ind2p1)*qith(ind2p1)*unsurtaurim 2383 2405 2384 2406 ! Precip term 2385 2407 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)) 2408 ! We assume that there is no solid precipitation outside thermals 2409 ! so the precipitation flux within thermals is equal to the precipitation flux 2410 ! at mesh-scale divided by thermals fraction 2411 2412 fp2=0. 2413 fp1=0. 2414 IF (fraca(ind2p1) .GT. 0.) THEN 2415 fp2=-qith(ind2p2)*rho(ind2p2)*vith(ind2p2)*fraca(ind2p2)! flux defined positive upward 2416 fp1=-qith(ind2p1)*rho(ind2p1)*vith(ind2p1)*fraca(ind2p1) 2417 ENDIF 2418 2419 precip_term=-1./deltazlev(ind2p1)*(fp2-fp1) 2389 2420 2390 2421 ! Calculation in a top-to-bottom loop 2391 2422 2392 2423 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) 2424 qi= 1./fm_therm(ind1,ind2p1)* & 2425 (deltazlev(ind2p1)*(-rim_term-dep_term-det_term-precip_term) + & 2426 fm_therm(ind1,ind2p2)*(qith(ind2p1))) 2427 END IF 2428 2429 ENDIF ! top thermals 2430 2431 qith(ind2)=MAX(0.,qi) 2403 2432 2404 2433 RETURN
Note: See TracChangeset
for help on using the changeset viewer.
