Changeset 4072 for LMDZ6/trunk/libf
- Timestamp:
- Feb 1, 2022, 6:34:34 PM (2 years ago)
- Location:
- LMDZ6/trunk/libf
- Files:
-
- 8 edited
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/trunk/libf/phylmd/calcratqs_multi_mod.F90
r4010 r4072 230 230 231 231 INTEGER :: i,k,nsrf 232 REAL :: ratiom, qsat2m, dqsatdT 233 REAL, DIMENSION(klon) :: xsi0 232 REAL, DIMENSION(klon) :: xsi0, ratiom, qsat2m, dqsatdT 234 233 REAL, DIMENSION (klon,klev) :: zlay 235 234 … … 239 238 !------------------------------------------- 240 239 241 DO i=1,klon242 240 243 ratiom =0.244 xsi0( i)=0.241 ratiom(:)=0. 242 xsi0(:)=0. 245 243 246 244 DO nsrf=1,nbsrf 247 CALL CALC_QSAT_ECMWF( t2m(i,nsrf),q2m(i,nsrf),paprs(i,1),RTT,0,.false.,qsat2m,dqsatdT)248 ratiom =ratiom+pctsrf(i,nsrf)*(q2m(i,nsrf)/qsat2m)249 xsi0( i)=xsi0(i)+pctsrf(i,nsrf)*((q2m(i,nsrf)/qsat2m-ratiom)**2)245 CALL CALC_QSAT_ECMWF(klon,t2m(:,nsrf),q2m(:,nsrf),paprs(:,1),RTT,0,.false.,qsat2m,dqsatdT) 246 ratiom(:)=ratiom(:)+pctsrf(:,nsrf)*(q2m(:,nsrf)/qsat2m(:)) 247 xsi0(:)=xsi0(:)+pctsrf(:,nsrf)*((q2m(:,nsrf)/qsat2m(:)-ratiom(:))**2) 250 248 END DO 251 249 252 xsi0(i)=sqrt(xsi0(i))/(ratiom+1E-6) 253 END DO 250 xsi0(:)=sqrt(xsi0(:))/(ratiom(:)+1E-6) 254 251 255 252 -
LMDZ6/trunk/libf/phylmd/cloudth_mod.F90
r3999 r4072 1539 1539 1540 1540 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1541 SUBROUTINE cloudth_mpc(klon,klev,ind2,mpc_bl_points, &1541 SUBROUTINE cloudth_mpc(klon,klev,ind2,mpc_bl_points, & 1542 1542 & temp,ztv,po,zqta,fraca,zpspsk,paprs,pplay,ztla,zthl, & 1543 1543 & ratqs,zqs,snowflux,qcloud,qincloud,icefrac,ctot,ctot_vol) … … 1554 1554 USE ioipsl_getin_p_mod, ONLY : getin_p 1555 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_LSCP1557 USE phys_local_var_mod, ONLY : qlth, qith 1556 USE lscp_tools_mod, ONLY: CALC_QSAT_ECMWF 1557 USE phys_local_var_mod, ONLY : qlth, qith, qsith 1558 1558 1559 1559 IMPLICIT NONE … … 1573 1573 1574 1574 INTEGER, INTENT(IN) :: klon,klev,ind2 1575 INTEGER, DIMENSION(klon,klev), INTENT(INOUT) :: mpc_bl_points ! 1 where BL MPC, 0 otherwise1575 1576 1576 1577 1577 REAL, DIMENSION(klon,klev), INTENT(IN) :: temp ! Temperature [K] … … 1589 1589 REAL, DIMENSION(klon,klev+1), INTENT(IN) :: snowflux ! snow flux at the interface of the layer [kg/m2/s] 1590 1590 1591 1591 INTEGER, DIMENSION(klon,klev), INTENT(INOUT) :: mpc_bl_points ! grid points with convective (thermals) mixed phase clouds 1592 1592 1593 REAL, DIMENSION(klon,klev), INTENT(OUT) :: ctot ! Cloud fraction [0-1] 1593 1594 REAL, DIMENSION(klon,klev), INTENT(OUT) :: ctot_vol ! Volume cloud fraction [0-1] … … 1600 1601 1601 1602 INTEGER itap,ind1,l,ig,iter,k 1602 LOGICAL flag_topthermals1603 1604 1605 REAL zqsatth(klon ,klev), zqsatenv(klon,klev)1603 INTEGER iflag_topthermals, niter 1604 1605 1606 REAL zqsatth(klon), zqsatenv(klon), zqsatenvonly(klon) 1606 1607 REAL sigma1(klon,klev) 1607 1608 REAL sigma2(klon,klev) … … 1614 1615 REAL cenv_vol(klon,klev) 1615 1616 REAL rneb(klon,klev) 1616 REAL zqenv(klon) 1617 1617 REAL zqenv(klon), zqth(klon), zqenvonly(klon) 1618 1618 REAL qsatmmussig1,qsatmmussig2,sqrtpi,sqrt2,sqrt2pi,pi 1619 1619 REAL rdd,cppd,Lv … … 1624 1624 REAL xth1,xth2,xenv1,xenv2,deltasth, deltasenv 1625 1625 REAL IntJ,IntI1,IntI2,IntI3,IntJ_CF,IntI1_CF,IntI3_CF,coeffqlenv,coeffqlth 1626 REAL Tbef,zdelta,qsatbef,zcor 1627 REAL qlbef,dqsatdt 1626 REAL zdelta,qsatbef,zcor 1627 REAL Tbefth(klon), Tbefenv(klon), Tbefenvonly(klon) 1628 REAL qlbef 1629 REAL dqsatenv(klon), dqsatth(klon), dqsatenvonly(klon) 1628 1630 REAL erf 1629 1631 REAL zpdf_sig(klon),zpdf_k(klon),zpdf_delta(klon) … … 1632 1634 REAL zrho(klon,klev) 1633 1635 REAL dz(klon,klev) 1634 REAL qslth, qsith, qslenv, alenvl, aenvl 1635 REAL sthi, sthl, althl, athl 1636 REAL qslenv(klon), qslth(klon) 1637 REAL alenvl, aenvl 1638 REAL sthi, sthl, sthil, althl, athl, althi, athi, sthlc, deltasthc, sigma2sc 1636 1639 REAL senvi, senvl, qbase, sbase, qliqth, qiceth 1637 REAL q imax, ttarget, stmp, cout, coutref1640 REAL qmax, ttarget, stmp, cout, coutref, fraci 1638 1641 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 REAL deltazlev_mpc(klev), temp_mpc(klev), pres_mpc(klev), fraca_mpc(klev+1), snowf_mpc(klev+1) 1642 1643 1643 1644 ! Modifty the saturation deficit PDF in thermals … … 1645 1646 REAL,SAVE :: C_mpc 1646 1647 !$OMP THREADPRIVATE(C_mpc) 1648 REAL, SAVE :: Ni,C_cap,Ei,d_top 1649 !$OMP THREADPRIVATE(Ni,C_cap,Ei,d_top) 1647 1650 ! Change the width of the PDF used for vertical subgrid scale heterogeneity 1648 1651 ! (J Jouhaud, JL Dufresne, JB Madeleine) … … 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 1804 zqenv(:)=(po(:)-fraca(:,ind2)*zqta(:,ind2))/(1.-fraca(:,ind2)) !qt = a*qtth + (1-a)*qtenv 1805 zqth(:)=zqta(:,ind2) 1806 zqenvonly(:)=po(:) 1807 1808 1809 CALL CALC_QSAT_ECMWF(klon,Tbefenvonly,zqenvonly,paprs(:,ind2),RTT,0,.false.,zqsatenvonly,dqsatenv) 1810 1811 CALL CALC_QSAT_ECMWF(klon,Tbefenv,zqenv,paprs(:,ind2),RTT,0,.false.,zqsatenv,dqsatenv) 1812 CALL CALC_QSAT_ECMWF(klon,Tbefenv,zqenv,paprs(:,ind2),RTT,1,.false.,qslenv,dqsatenv) 1813 1814 CALL CALC_QSAT_ECMWF(klon,Tbefth,zqth,paprs(:,ind2),RTT,1,.false.,qslth,dqsatth) 1815 CALL CALC_QSAT_ECMWF(klon,Tbefth,zqth,paprs(:,ind2),RTT,2,.false.,qsith(:,ind2),dqsatth) 1816 CALL CALC_QSAT_ECMWF(klon,Tbefth,zqth,paprs(:,ind2),RTT,0,.false.,zqsatth,dqsatth) 1817 1818 1819 DO ind1=1,klon 1775 1820 1776 1821 … … 1780 1825 ! Environment: 1781 1826 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 1827 1828 IF (Tbefenv(ind1) .GE. RTT) THEN 1789 1829 Lv=RLVTT 1790 1830 ELSE … … 1793 1833 1794 1834 1795 alenv=(0.622*Lv*zqsatenv(ind1 ,ind2))/(rdd*zthl(ind1,ind2)**2) !qsl, p841835 alenv=(0.622*Lv*zqsatenv(ind1))/(rdd*zthl(ind1,ind2)**2) !qsl, p84 1796 1836 aenv=1./(1.+(alenv*Lv/cppd)) !al, p84 1797 senv=aenv*(po(ind1)-zqsatenv(ind1 ,ind2)) !s, p841837 senv=aenv*(po(ind1)-zqsatenv(ind1)) !s, p84 1798 1838 1799 1839 ! For MPCs: 1800 1840 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) 1841 alenvl=(0.622*RLVTT*qslenv(ind1))/(rdd*zthl(ind1,ind2)**2) 1803 1842 aenvl=1./(1.+(alenv*Lv/cppd)) 1804 senvl=aenvl*(po(ind1)-qslenv) 1843 senvl=aenvl*(po(ind1)-qslenv(ind1)) 1844 senvi=senv 1805 1845 ENDIF 1806 1846 … … 1808 1848 ! Thermals: 1809 1849 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 1850 1851 IF (Tbefth(ind1) .GE. RTT) THEN 1815 1852 Lv=RLVTT 1816 1853 ELSE … … 1819 1856 1820 1857 1821 alth=(0.622*Lv*zqsatth(ind1 ,ind2))/(rdd*ztla(ind1,ind2)**2)1858 alth=(0.622*Lv*zqsatth(ind1))/(rdd*ztla(ind1,ind2)**2) 1822 1859 ath=1./(1.+(alth*Lv/cppd)) 1823 sth=ath*(zqta(ind1,ind2)-zqsatth(ind1 ,ind2))1860 sth=ath*(zqta(ind1,ind2)-zqsatth(ind1)) 1824 1861 1825 1862 ! For MPCs: 1826 1863 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) 1864 althi=alth 1865 althl=(0.622*RLVTT*qslth(ind1))/(rdd*ztla(ind1,ind2)**2) 1866 athl=1./(1.+(alth*RLVTT/cppd)) 1867 athi=alth 1868 sthl=athl*(zqta(ind1,ind2)-qslth(ind1)) 1869 sthi=athi*(zqta(ind1,ind2)-qsith(ind1,ind2)) 1870 sthil=athi*(zqta(ind1,ind2)-qslth(ind1)) 1833 1871 ENDIF 1834 1872 … … 1923 1961 ctot(ind1,ind2)=0. 1924 1962 ctot_vol(ind1,ind2)=0. 1925 qcloud(ind1)=zqsatenv(ind1 ,ind2)1963 qcloud(ind1)=zqsatenv(ind1) 1926 1964 qincloud(ind1)=0. 1927 1965 ELSE … … 1941 1979 !........... 1942 1980 1981 qiceth=0. 1943 1982 deltazlev_mpc=dz(ind1,:) 1944 1983 temp_mpc=ztla(ind1,:)*zpspsk(ind1,:) … … 1946 1985 fraca_mpc=fraca(ind1,:) 1947 1986 snowf_mpc=snowflux(ind1,:) 1948 qth_mpc=zqta(ind1,:) 1949 flag_topthermals=.FALSE. 1987 iflag_topthermals=0 1950 1988 IF ((mpc_bl_points(ind1,ind2) .EQ. 1) .AND. (mpc_bl_points(ind1,ind2+1) .EQ. 0)) THEN 1951 flag_topthermals = .TRUE. 1989 iflag_topthermals = 1 1990 ELSE IF ((mpc_bl_points(ind1,ind2) .EQ. 1) .AND. (mpc_bl_points(ind1,ind2+1) .EQ. 1) & 1991 .AND. (mpc_bl_points(ind1,ind2+2) .EQ. 0) ) THEN 1992 iflag_topthermals = 2 1993 ELSE 1994 iflag_topthermals = 0 1952 1995 ENDIF 1953 1996 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) 1997 CALL ICE_MPC_BL_CLOUDS(ind1,ind2,klev,Ni,Ei,C_cap,d_top,iflag_topthermals,temp_mpc,pres_mpc,zqta(ind1,:),qsith(ind1,:),qlth(ind1,:),deltazlev_mpc,snowf_mpc,fraca_mpc,qith(ind1,:)) 1998 1999 ! qmax calculation 2000 sigma2s=(sigma2s_factor*((MAX((sthl-senvl),0.)**2)**0.5)/((fraca(ind1,ind2)+0.02)**sigma2s_power))+0.002*zqta(ind1,ind2) 1971 2001 deltasth=athl*vert_alpha_th*sigma2s 1972 1973 ! Liquid phase1974 !.............1975 2002 xth1=-(sthl+deltasth)/(sqrt(2.)*sigma2s) 1976 xth2=-(sthl-deltasth)/(sqrt(2.)*sigma2s) 2003 xth2=-(sthl-deltasth)/(sqrt(2.)*sigma2s) 1977 2004 exp_xth1 = exp(-1.*xth1**2) 1978 2005 exp_xth2 = exp(-1.*xth2**2) … … 1984 2011 IntI1_CF=((sthl+deltasth)*(erf(xth2)-erf(xth1)))/(4*deltasth) 1985 2012 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) 2013 qmax=MAX(IntJ+IntI1+IntI2+IntI3,0.) 2014 2015 2016 ! Liquid phase 2017 !................ 2018 ! We account for the effect of ice crystals in thermals on sthl 2019 ! and on the width of the distribution 2020 2021 sthlc=sthl*1./(1.+C_mpc*qith(ind1,ind2)) & 2022 + (1.-1./(1.+C_mpc*qith(ind1,ind2))) * athl*(qsith(ind1,ind2)-qslth(ind1)) 2023 2024 sigma2sc=(sigma2s_factor*((MAX((sthlc-senvl),0.)**2)**0.5)/((fraca(ind1,ind2)+0.02)**sigma2s_power))+0.002*zqta(ind1,ind2) 2025 deltasthc=athl*vert_alpha_th*sigma2sc 2026 2027 2028 xth1=-(sthlc+deltasthc)/(sqrt(2.)*sigma2sc) 2029 xth2=-(sthlc-deltasthc)/(sqrt(2.)*sigma2sc) 1991 2030 exp_xth1 = exp(-1.*xth1**2) 1992 2031 exp_xth2 = exp(-1.*xth2**2) 1993 IntJ=0.5*sth i*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth22032 IntJ=0.5*sthlc*(1-erf(xth2))+(sigma2sc/sqrt2pi)*exp_xth2 1994 2033 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 2034 IntI1=(((sthlc+deltasthc)**2+(sigma2sc)**2)/(8*deltasthc))*(erf(xth2)-erf(xth1)) 2035 IntI2=(sigma2sc**2/(4*deltasthc*sqrtpi))*(xth1*exp_xth1-xth2*exp_xth2) 2036 IntI3=((sqrt2*sigma2sc*(sthlc+deltasthc))/(4*sqrtpi*deltasthc))*(exp_xth1-exp_xth2) 2037 IntI1_CF=((sthlc+deltasthc)*(erf(xth2)-erf(xth1)))/(4*deltasthc) 2038 IntI3_CF=(sqrt2*sigma2sc*(exp_xth1-exp_xth2))/(4*sqrtpi*deltasthc) 2039 qliqth=IntJ+IntI1+IntI2+IntI3 2040 2041 qlth(ind1,ind2)=MAX(0.,qliqth) 2003 2042 2004 2043 ! 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)) 2044 2008 2045 qcth(ind1,ind2)=qlth(ind1,ind2)+qith(ind1,ind2) 2009 2046 2047 2048 ! consistency with subgrid distribution 2049 2050 IF ((qcth(ind1,ind2) .GT. qmax) .AND. (qcth(ind1,ind2) .GT. 0)) THEN 2051 fraci=qith(ind1,ind2)/qcth(ind1,ind2) 2052 qcth(ind1,ind2)=qmax 2053 qith(ind1,ind2)=fraci*qmax 2054 qlth(ind1,ind2)=(1.-fraci)*qmax 2055 ENDIF 2056 2057 ! Cloud Fraction 2058 !............... 2010 2059 ! calculation of qbase which is the value of the water vapor within mixed phase clouds 2011 2060 ! such that the total water in cloud = qbase+qliqth+qiceth … … 2015 2064 2016 2065 ttarget=qcth(ind1,ind2) 2017 mini=athl*(qsith-qslth) 2018 maxi=0. 2066 mini= athl*(qsith(ind1,ind2)-qslth(ind1)) 2067 maxi= 0. !athl*(3.*sqrt(sigma2s)) 2068 niter=20 2019 2069 pas=(maxi-mini)/niter 2020 2070 stmp=mini … … 2031 2081 ENDIF 2032 2082 ENDDO 2033 qbase=MAX(0., sbase/athl+qslth )2083 qbase=MAX(0., sbase/athl+qslth(ind1)) 2034 2084 2035 2085 ! surface cloud fraction in thermals … … 2053 2103 IntI1=(((sthl+deltasth-sbase)**2+(sigma2s)**2)/(8*deltasth))*(erf(xth2)-erf(xth1)) 2054 2104 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)2105 IntI3=((sqrt2*sigma2s*(sthl+deltasth))/(4*sqrtpi*deltasth))*(exp_xth1-exp_xth2) 2056 2106 IntI1_CF=((sthl-sbase+deltasth)*(erf(xth2)-erf(xth1)))/(4*deltasth) 2057 2107 IntI3_CF=(sqrt2*sigma2s*(exp_xth1-exp_xth2))/(4*sqrtpi*deltasth) … … 2059 2109 ENDIF 2060 2110 cth_vol(ind1,ind2)=MIN(MAX(0.,cth_vol(ind1,ind2)),1.) 2111 2112 2061 2113 2062 2114 ! Environment … … 2110 2162 2111 2163 2112 ! Thermals + environment 2164 ! Thermals + environment 2165 ! ======================= 2113 2166 ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2) 2114 2167 qctot(ind1,ind2)=fraca(ind1,ind2)*qcth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qcenv(ind1,ind2) … … 2125 2178 ctot_vol(ind1,ind2)=0. 2126 2179 qincloud(ind1)=0. 2127 qcloud(ind1)=zqsatenv(ind1 ,ind2)2180 qcloud(ind1)=zqsatenv(ind1) 2128 2181 ELSE 2129 2182 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 )2183 +(1.-1.*fraca(ind1,ind2))*(qcenv(ind1,ind2)/cenv(ind1,ind2)+qslenv(ind1)) 2131 2184 qincloud(ind1)=MAX(fraca(ind1,ind2)*(qcth(ind1,ind2)/cth(ind1,ind2)) & 2132 2185 +(1.-1.*fraca(ind1,ind2))*(qcenv(ind1,ind2)/cenv(ind1,ind2)),0.) … … 2139 2192 2140 2193 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 2194 2195 2196 IF (Tbefenvonly(ind1) .GE. RTT) THEN 2148 2197 Lv=RLVTT 2149 2198 ELSE … … 2153 2202 2154 2203 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)2204 alenv=(0.622*Lv*zqsatenvonly(ind1))/(rdd*zthl(ind1,ind2)**2) 2156 2205 aenv=1./(1.+(alenv*Lv/cppd)) 2157 senv=aenv*(po(ind1)-zqsatenv (ind1,ind2))2206 senv=aenv*(po(ind1)-zqsatenvonly(ind1)) 2158 2207 sth=0. 2159 2208 2160 sigma1s=ratqs(ind1,ind2)*zqenv (ind1)2209 sigma1s=ratqs(ind1,ind2)*zqenvonly(ind1) 2161 2210 sigma2s=0. 2162 2211 … … 2169 2218 IF (ctot(ind1,ind2).LT.1.e-3) THEN 2170 2219 ctot(ind1,ind2)=0. 2171 qcloud(ind1)=zqsatenv (ind1,ind2)2220 qcloud(ind1)=zqsatenvonly(ind1) 2172 2221 qincloud(ind1)=0. 2173 2222 ELSE 2174 qcloud(ind1)=qctot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv (ind1,ind2)2223 qcloud(ind1)=qctot(ind1,ind2)/ctot(ind1,ind2)+zqsatenvonly(ind1) 2175 2224 qincloud(ind1)=MAX(qctot(ind1,ind2)/ctot(ind1,ind2),0.) 2176 2225 ENDIF … … 2194 2243 2195 2244 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 2196 SUBROUTINE ICE_MPC_BL_CLOUDS(ind1,ind2,klev, flag_topthermals,temp,pres,qth,qlth,qith,deltazlev,snowf,fraca,qi)2245 SUBROUTINE ICE_MPC_BL_CLOUDS(ind1,ind2,klev,Ni,Ei,C_cap,d_top,iflag_topthermals,temp,pres,qth,qsith,qlth,deltazlev,snowf,fraca,qith) 2197 2246 2198 2247 ! parameterization of ice for boundary … … 2260 2309 IMPLICIT none 2261 2310 2262 2263 2311 INCLUDE "YOMCST.h" 2264 2312 INCLUDE "nuage.h" 2265 2313 2266 INTEGER, INTENT(IN) :: ind1,ind2, klev ! horizontal and vertical indices and dimensions 2267 LOGICAL, INTENT(IN) :: flag_topthermals ! uppermost layer of thermals ? 2314 INTEGER, INTENT(IN) :: ind1,ind2, klev ! horizontal and vertical indices and dimensions 2315 INTEGER, INTENT(IN) :: iflag_topthermals ! uppermost layer of thermals ? 2316 REAL, INTENT(IN) :: Ni ! ice number concentration [m-3] 2317 REAL, INTENT(IN) :: Ei ! ice-droplet collision efficiency 2318 REAL, INTENT(IN) :: C_cap ! ice crystal capacitance 2319 REAL, INTENT(IN) :: d_top ! cloud-top ice crystal mixing parameter 2268 2320 REAL, DIMENSION(klev), INTENT(IN) :: temp ! temperature [K] within thermals 2269 2321 REAL, DIMENSION(klev), INTENT(IN) :: pres ! pressure [Pa] 2270 2322 REAL, DIMENSION(klev), INTENT(IN) :: qth ! mean specific water content in thermals [kg/kg] 2323 REAL, DIMENSION(klev), INTENT(IN) :: qsith ! saturation specific humidity wrt ice in thermals [kg/kg] 2271 2324 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 2325 REAL, DIMENSION(klev), INTENT(IN) :: deltazlev ! layer thickness [m] 2274 2326 REAL, DIMENSION(klev+1), INTENT(IN) :: snowf ! snow flux at the upper inferface 2275 2327 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) 2328 REAL, DIMENSION(klev), INTENT(INOUT) :: qith ! condensed ice water , thermals [kg/kg] 2284 2329 2285 2330 … … 2287 2332 REAL rho(klev) 2288 2333 REAL unsurtaudet, unsurtaustardep, unsurtaurim 2289 REAL q sl, qsi, dqs, AA, BB, Ka, Dv, rhoi2290 REAL p0, t02334 REAL qi, AA, BB, Ka, Dv, rhoi 2335 REAL p0, t0, fp1, fp2 2291 2336 REAL alpha, flux_term 2292 2337 REAL det_term, precip_term, rim_term, dep_term 2293 2338 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 2339 2317 2340 ind2p1=ind2+1 2318 2341 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 2342 2327 2343 rho=pres/temp/RD ! air density kg/m3 … … 2334 2350 2335 2351 2336 IF ( flag_topthermals) THEN ! uppermost thermals level, solve a third order polynomial with Cardan's method2352 IF (iflag_topthermals .GT. 0) THEN ! uppermost thermals levels 2337 2353 2338 2354 Dv=0.0001*0.211*(p0/pres(ind2))*((temp(ind2)/T0)**1.94) ! water vapor diffusivity in air, SI 2339 2355 2340 2356 ! Detrainment term: 2357 2341 2358 unsurtaudet=detr_therm(ind1,ind2)/rho(ind2)/deltazlev(ind2) 2342 2359 2343 ! vertical flux2344 2345 flux_term=d_top*fm_therm(ind1,ind2)/deltazlev(ind2)2346 2360 2347 2361 ! 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 2362 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) 2363 BB=1./(rho(ind2)*Dv*qsith(ind2)) 2364 unsurtaustardep=C_cap*(Ni**0.66)*(qth(ind2)-qsith(ind2))/qsith(ind2) & 2365 *4.*RPI/(AA+BB)*(6.*rho(ind2)/rhoi/RPI/Gamma(4.))**(0.33) 2353 2366 2354 2367 ! Riming term neglected at this level 2355 2368 !unsurtaurim=rho(ind2)*alpha*3./rhoi/2.*Ei*qlth(ind2)*((p0/pres(ind2))**0.4) 2356 2369 2357 qi= rho(ind2)*unsurtaustardep/MAX((rho(ind2)*unsurtaudet-flux_term),1E-12)2370 qi=fraca(ind2)*unsurtaustardep/MAX((d_top*unsurtaudet),1E-12) 2358 2371 qi=MAX(qi,0.)**(3./2.) 2359 2372 … … 2370 2383 ! Deposition term 2371 2384 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 2385 AA=RLSTT/Ka/temp(ind2p1)*(RLSTT/RV/temp(ind2p1)-1.) 2375 BB=1./(rho(ind2p1)*Dv*qsi )2376 unsurtaustardep=C_cap*(Ni**0.66)*(q sl-qsi)/qsi*4.*RPI/(AA+BB)*(6.*rho(ind2p1)/rhoi/RPI/Gamma(4.))**(0.33)2377 dep_term=rho(ind2p1)* (qith(ind2p1)**0.33)*unsurtaustardep2386 BB=1./(rho(ind2p1)*Dv*qsith(ind2p1)) 2387 unsurtaustardep=C_cap*(Ni**0.66)*(qth(ind2p1)-qsith(ind2p1))/qsith(ind2p1)*4.*RPI/(AA+BB)*(6.*rho(ind2p1)/rhoi/RPI/Gamma(4.))**(0.33) 2388 dep_term=rho(ind2p1)*fraca(ind2p1)*(qith(ind2p1)**0.33)*unsurtaustardep 2378 2389 2379 2390 ! Riming term 2380 2391 2381 2392 unsurtaurim=rho(ind2p1)*alpha*3./rhoi/2.*Ei*qlth(ind2p1)*((p0/pres(ind2p1))**0.4) 2382 rim_term=rho(ind2p1)* qith(ind2p1)*unsurtaurim2393 rim_term=rho(ind2p1)*fraca(ind2p1)*qith(ind2p1)*unsurtaurim 2383 2394 2384 2395 ! Precip term 2385 2396 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)) 2397 ! We assume that there is no solid precipitation outside thermals 2398 ! so the precipitation flux within thermals is equal to the precipitation flux 2399 ! at mesh-scale divided by thermals fraction 2400 2401 fp2=0. 2402 fp1=0. 2403 IF (fraca(ind2p1) .GT. 0.) THEN 2404 fp2=-snowf(ind2p2)/fraca(ind2p1) ! flux defined positive upward 2405 fp1=-snowf(ind2p1)/fraca(ind2p1) 2406 ENDIF 2407 2408 precip_term=-1./deltazlev(ind2p1)*(fp2-fp1) 2389 2409 2390 2410 ! Calculation in a top-to-bottom loop 2391 2411 2392 2412 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) 2413 qi= 1./fm_therm(ind1,ind2p1)* & 2414 (deltazlev(ind2p1)*(-rim_term-dep_term-det_term-precip_term) + & 2415 fm_therm(ind1,ind2p2)*(qith(ind2p1))) 2416 END IF 2417 2418 ENDIF ! top thermals 2419 2420 qith(ind2)=MAX(0.,qi) 2403 2421 2404 2422 RETURN -
LMDZ6/trunk/libf/phylmd/conf_phys_m.F90
r4062 r4072 175 175 INTEGER,SAVE :: iflag_clw_omp 176 176 REAL,SAVE :: cld_lc_lsc_omp,cld_lc_con_omp,cld_tau_lsc_omp,cld_tau_con_omp 177 REAL,SAVE :: ffallv_lsc_omp, ffallv_con_omp,coef_eva_omp 177 REAL,SAVE :: ffallv_lsc_omp, ffallv_con_omp,coef_eva_omp,coef_eva_i_omp 178 178 LOGICAL,SAVE :: reevap_ice_omp 179 179 INTEGER,SAVE :: iflag_pdf_omp … … 1131 1131 CALL getin('coef_eva',coef_eva_omp) 1132 1132 ! 1133 !Config Key = coef_eva_i 1134 !Config Desc = 1135 !Config Def = 2.e-5 1136 !Config Help = 1137 ! 1138 coef_eva_i_omp = coef_eva_omp 1139 CALL getin('coef_eva_i',coef_eva_i_omp) 1140 ! 1133 1141 !Config Key = reevap_ice 1134 1142 !Config Desc = … … 2492 2500 ffallv_con = ffallv_con_omp 2493 2501 coef_eva = coef_eva_omp 2502 coef_eva_i = coef_eva_i_omp 2494 2503 reevap_ice = reevap_ice_omp 2495 2504 iflag_pdf = iflag_pdf_omp … … 2914 2923 WRITE(lunout,*) ' ffallv_con = ', ffallv_con 2915 2924 WRITE(lunout,*) ' coef_eva = ', coef_eva 2925 WRITE(lunout,*) ' coef_eva_i = ', coef_eva_i 2916 2926 WRITE(lunout,*) ' reevap_ice = ', reevap_ice 2917 2927 WRITE(lunout,*) ' iflag_pdf = ', iflag_pdf -
LMDZ6/trunk/libf/phylmd/fisrtilp.h
r2415 r4072 11 11 REAL cld_tau_lsc,cld_tau_con 12 12 REAL ffallv_lsc,ffallv_con 13 REAL coef_eva 13 REAL coef_eva,coef_eva_i 14 14 LOGICAL reevap_ice 15 15 INTEGER iflag_pdf … … 25 25 & ,ffallv_con & 26 26 & ,coef_eva & 27 & ,coef_eva_i & 27 28 & ,reevap_ice & 28 29 & ,iflag_fisrtilp_qsat & -
LMDZ6/trunk/libf/phylmd/lscp_mod.F90
r4062 r4072 6 6 7 7 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 8 SUBROUTINE LSCP(dtime,missing_val, & 8 SUBROUTINE LSCP(dtime,missing_val, & 9 9 paprs,pplay,t,q,ptconv,ratqs, & 10 d_t, d_q, d_ql, d_qi, rneb, rneb_seri, & 10 d_t, d_q, d_ql, d_qi, rneb, rneb_seri, & 11 11 radliq, radicefrac, rain, snow, & 12 12 pfrac_impa, pfrac_nucl, pfrac_1nucl, & … … 25 25 ! 26 26 ! This code is a new version of the fisrtilp.F90 routine, which is itself a 27 ! fusionof 'first' (superrsaturation physics, P. LeVan K. Laval)27 ! merge of 'first' (superrsaturation physics, P. LeVan K. Laval) 28 28 ! and 'ilp' (il pleut, L. Li) 29 29 ! … … 126 126 INTEGER, INTENT(IN) :: iflag_ice_thermo! flag to activate the ice thermodynamics 127 127 ! CR: if iflag_ice_thermo=2, only convection is active 128 LOGICAL, INTENT(IN) :: ok_ice_sursat ! flag to determine if ice sursaturation is activated 128 LOGICAL, INTENT(IN) :: ok_ice_sursat ! flag to determine if ice sursaturation is activated 129 129 130 130 LOGICAL, DIMENSION(klon,klev), INTENT(IN) :: ptconv ! grid points where deep convection scheme is active … … 138 138 REAL, DIMENSION(klon,klev), INTENT(IN) :: zpspsk ! exner potential (p/100000)**(R/cp) 139 139 REAL, DIMENSION(klon,klev), INTENT(IN) :: ztla ! liquid temperature within thermals [K] 140 REAL, DIMENSION(klon,klev), INTENT(INOUT) :: zthl 140 REAL, DIMENSION(klon,klev), INTENT(INOUT) :: zthl ! liquid potential temperature [K] 141 141 142 142 ! INPUT/OUTPUT variables … … 145 145 REAL, DIMENSION(klon,klev), INTENT(INOUT):: ratqs ! function of pressure that sets the large-scale 146 146 ! cloud PDF (sigma=ratqs*qt) 147 147 148 148 149 ! Input sursaturation en glace … … 204 205 PARAMETER (ztfondue=278.15) 205 206 206 REAL, SAVE :: rain_int_min=0.001 ! Minimum local rain intensity [mm/s] before the decrease in associatesprecipitation fraction207 REAL, SAVE :: rain_int_min=0.001 ! Minimum local rain intensity [mm/s] before the decrease in associated precipitation fraction 207 208 !$OMP THREADPRIVATE(rain_int_min) 208 209 209 210 210 LOGICAL, SAVE :: ok_radliq_precip=.false. ! Inclusion of all precipitation in radliq (cloud water seen by radiation) 211 !$OMP THREADPRIVATE(ok_radliq_precip) 212 213 214 215 211 216 ! LOCAL VARIABLES: 212 217 !---------------- 213 218 214 219 215 REAL qsl , qsi220 REAL qsl(klon), qsi(klon) 216 221 REAL zct, zcl 217 222 REAL ctot(klon,klev) 218 223 REAL ctot_vol(klon,klev) 219 INTEGER mpc_bl_points(klon,klev)220 224 REAL zqs(klon), zdqs(klon), zdelta, zcor, zcvm5 221 225 REAL zdqsdT_raw(klon) … … 230 234 REAL zifl(klon), zifln(klon), zqev0,zqevi, zqevti 231 235 REAL zoliq(klon), zcond(klon), zq(klon), zqn(klon) 232 REAL zoliq p(klon), zoliqi(klon)236 REAL zoliql(klon), zoliqi(klon) 233 237 REAL zt(klon) 234 REAL zdz(klon),zrho(klon), ztot, zrhol(klon)238 REAL zdz(klon),zrho(klon),iwc(klon) 235 239 REAL zchau,zfroi,zfice(klon),zneb(klon),znebprecip(klon) 236 REAL zmelt,z pluie,zice240 REAL zmelt,zrain,zsnow,zprecip 237 241 REAL dzfice(klon) 238 REAL zsolid,qtot,dqsl,dqsi 242 REAL zsolid 243 REAL qtot(klon), qzero(klon) 244 REAL dqsl(klon),dqsi(klon) 239 245 REAL smallestreal 240 246 ! Variables for Bergeron process … … 262 268 REAL tot_zneb(klon), tot_znebn(klon), d_tot_zneb(klon) 263 269 REAL d_znebprecip_clr_cld(klon), d_znebprecip_cld_clr(klon) 264 REAL velo(klon )270 REAL velo(klon,klev), vr 265 271 REAL qlmpc, qimpc, rnebmpc 266 REAL radliqi(klon,klev) 272 REAL radliqi(klon,klev), radliql(klon,klev) 267 273 268 274 INTEGER i, k, n, kk 269 275 INTEGER n_i(klon), iter 270 276 INTEGER ncoreczq 277 INTEGER mpc_bl_points(klon,klev) 271 278 INTEGER,SAVE :: itap=0 272 279 !$OMP THREADPRIVATE(itap) … … 334 341 CALL getin_p('seuil_neb',seuil_neb) 335 342 CALL getin_p('rain_int_min',rain_int_min) 343 CALL getin_p('ok_radliq_precip',ok_radliq_precip) 336 344 WRITE(lunout,*) 'lscp, ninter:', ninter 337 345 WRITE(lunout,*) 'lscp, iflag_evap_prec:', iflag_evap_prec 338 346 WRITE(lunout,*) 'lscp, seuil_neb:', seuil_neb 339 347 WRITE(lunout,*) 'lscp, rain_int_min:', rain_int_min 348 WRITE(lunout,*) 'lscp, ok_radliq_precip:', ok_radliq_precip 340 349 341 350 ! check for precipitation sub-time steps … … 397 406 tot_znebn(:) = 0.0 398 407 d_tot_zneb(:) = 0.0 408 qzero(:)=0.0 399 409 400 410 !--ice sursaturation … … 468 478 ! -------------------------------------------------------------------- 469 479 ! A part of the precipitation coming from above is evaporated/sublimated. 470 !471 480 ! -------------------------------------------------------------------- 472 481 … … 474 483 IF (iflag_evap_prec.GE.1) THEN 475 484 485 ! Calculation of saturation specific humidity 486 ! depending on temperature: 487 CALL CALC_QSAT_ECMWF(klon,zt(:),qzero(:),pplay(:,k),RTT,0,.false.,zqs(:),zdqs(:)) 488 ! wrt liquid water 489 CALL CALC_QSAT_ECMWF(klon,zt(:),qzero(:),pplay(:,k),RTT,1,.false.,qsl(:),dqsl(:)) 490 ! wrt ice 491 CALL CALC_QSAT_ECMWF(klon,zt(:),qzero(:),pplay(:,k),RTT,2,.false.,qsi(:),dqsi(:)) 476 492 477 493 DO i = 1, klon … … 480 496 IF (zrfl(i)+zifl(i).GT.0.) THEN 481 497 482 CALL CALC_QSAT_ECMWF(zt(i),0.,pplay(i,k),RTT,0,.false.,zqs(i),zdqs(i))483 484 498 ! LTP: we only account for precipitation evaporation in the clear-sky (iflag_evap_prec=4). 485 499 IF (iflag_evap_prec.EQ.4) THEN … … 503 517 ENDIF 504 518 505 ! A. JAM 506 ! We consider separately qsatice and qstal 507 ! qsat wrt liquid phase according to thermcep 508 CALL CALC_QSAT_ECMWF(zt(i),0.,pplay(i,k),RTT,1,.false.,qsl,dqsl) 509 510 519 511 520 ! Evaporation of liquid precipitation coming from above 512 521 ! dP/dz=beta*(1-q/qsat)*sqrt(P) … … 515 524 516 525 IF (iflag_evap_prec.EQ.3) THEN 517 zqevt = znebprecip(i)*coef_eva*(1.0-zq(i)/qsl ) &526 zqevt = znebprecip(i)*coef_eva*(1.0-zq(i)/qsl(i)) & 518 527 *SQRT(zrfl(i)/max(1.e-4,znebprecip(i))) & 519 528 *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG 520 529 ELSE IF (iflag_evap_prec.EQ.4) THEN 521 zqevt = znebprecipclr(i)*coef_eva*(1.0-zq(i)/qsl ) &530 zqevt = znebprecipclr(i)*coef_eva*(1.0-zq(i)/qsl(i)) & 522 531 *SQRT(zrfl(i)/max(1.e-8,znebprecipclr(i))) & 523 532 *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG 524 533 ELSE 525 zqevt = 1.*coef_eva*(1.0-zq(i)/qsl )*SQRT(zrfl(i)) &534 zqevt = 1.*coef_eva*(1.0-zq(i)/qsl(i))*SQRT(zrfl(i)) & 526 535 *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG 527 536 ENDIF … … 531 540 *RG*dtime/(paprs(i,k)-paprs(i,k+1)) 532 541 533 ! qsat wrt ice according to thermcep534 CALL CALC_QSAT_ECMWF(zt(i),0.,pplay(i,k),RTT,2,.false.,qsi,dqsi)535 536 542 ! sublimation of the solid precipitation coming from above 537 543 IF (iflag_evap_prec.EQ.3) THEN 538 zqevti = znebprecip(i)*coef_eva *(1.0-zq(i)/qsi) &544 zqevti = znebprecip(i)*coef_eva_i*(1.0-zq(i)/qsi(i)) & 539 545 *SQRT(zifl(i)/max(1.e-4,znebprecip(i))) & 540 546 *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG 541 547 ELSE IF (iflag_evap_prec.EQ.4) THEN 542 zqevti = znebprecipclr(i)*coef_eva *(1.0-zq(i)/qsi) &548 zqevti = znebprecipclr(i)*coef_eva_i*(1.0-zq(i)/qsi(i)) & 543 549 *SQRT(zifl(i)/max(1.e-8,znebprecipclr(i))) & 544 550 *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG 545 551 ELSE 546 zqevti = 1.*coef_eva *(1.0-zq(i)/qsi)*SQRT(zifl(i)) &552 zqevti = 1.*coef_eva_i*(1.0-zq(i)/qsi(i))*SQRT(zifl(i)) & 547 553 *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG 548 554 ENDIF … … 566 572 567 573 568 ! EV: agrees with JL's comments below so correct and comment569 ! JLD: I do not understand the lines below. We distribute the liquid570 ! and solid parts of the precipitation even if the layer does not get571 ! saturated. To my opinion, we should all replace with:572 ! zqev=zqevt573 ! zqevi=zqevti574 ! IF (zqevt+zqevti.GT.0.) THEN575 ! zqev=MIN(zqev0*zqevt/(zqevt+zqevti),zqevt)576 ! zqevi=MIN(zqev0*zqevti/(zqevt+zqevti),zqevti)577 ! ELSE578 ! zqev=0.579 ! zqevi=0.580 ! ENDIF581 ! ENDIF582 583 574 ! New solid and liquid precipitation fluxes 584 575 zrfln(i) = Max(0.,zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1)) & … … 659 650 ! Calculation of qsat, L/Cp*dqsat/dT and ncoreczq counter 660 651 !------------------------------------------------------- 661 662 DO i = 1, klon 652 653 qtot(:)=zq(:)+zmqc(:) 654 CALL CALC_QSAT_ECMWF(klon,zt(:),qtot(:),pplay(:,k),RTT,0,.false.,zqs(:),zdqs(:)) 655 DO i = 1, klon 663 656 zdelta = MAX(0.,SIGN(1.,RTT-zt(i))) 664 qtot=zq(i)+zmqc(i)665 CALL CALC_QSAT_ECMWF(zt(i),qtot,pplay(i,k),RTT,0,.false.,zqs(i),zdqs(i))666 657 zdqsdT_raw(i) = zdqs(i)*RCPD*(1.0+RVTMP2*zq(i)) / (RLVTT*(1.-zdelta) + RLSTT*zdelta) 667 668 658 IF (zq(i) .LT. 1.e-15) THEN 669 659 ncoreczq=ncoreczq+1 … … 671 661 ENDIF 672 662 673 663 ENDDO 674 664 675 665 … … 777 767 DO i=1,klon 778 768 779 ! convergence = .true. sinceconvergence is not satisfied769 ! convergence = .true. until when convergence is not satisfied 780 770 convergence(i)=ABS(DT(i)).GT.DDT0 781 771 … … 796 786 797 787 ! Rneb, qzn and zcond for lognormal PDFs 798 qtot=zq(i)+zmqc(i) 799 CALL CALC_QSAT_ECMWF(Tbef(i),qtot,pplay(i,k),RTT,0,.false.,zqs(i),zdqs(i)) 800 CALL CALC_GAMMASAT(Tbef(i),qtot,pplay(i,k),gammasat(i),dgammasatdt(i)) 801 802 ! saturation may occur at a humidity different from qsat (gamma qsat), so gamma correction for dqs 803 zdqs(i) = gammasat(i)*zdqs(i)+zqs(i)*dgammasatdt(i) 788 qtot(i)=zq(i)+zmqc(i) 789 790 ENDIF 791 792 ENDDO 793 794 ! Calculation of saturation specific humidity and ce fraction 795 CALL CALC_QSAT_ECMWF(klon,Tbef(:),qtot(:),pplay(:,k),RTT,0,.false.,zqs(:),zdqs(:)) 796 CALL CALC_GAMMASAT(klon,Tbef(:),qtot(:),pplay(:,k),gammasat(:),dgammasatdt(:)) 797 ! saturation may occur at a humidity different from qsat (gamma qsat), so gamma correction for dqs 798 zdqs(:) = gammasat(:)*zdqs(:)+zqs(:)*dgammasatdt(:) 799 CALL ICEFRAC_LSCP(klon, zt(:),pplay(:,k)/paprs(:,1),zfice(:),dzfice(:)) 800 801 802 803 DO i=1,klon 804 805 IF ((convergence(i) .OR. (n_i(i) .EQ. 0)) .AND. lognormale(i)) THEN 804 806 805 807 zpdf_sig(i)=ratqs(i,k)*zq(i) … … 826 828 ENDIF 827 829 828 rnebss(i,k)=0.0 !--a jout OB (necessaire car boucle de convergence sur le temps)830 rnebss(i,k)=0.0 !--added by OB (needed because of the convergence loop on time) 829 831 fcontrN(i,k)=0.0 !--idem 830 832 fcontrP(i,k)=0.0 !--idem 831 833 qss(i,k)=0.0 !--idem 832 834 833 835 ELSE 836 834 837 !------------------------------------ 835 ! SURSATURATION EN GLACE838 ! ICE SUPERSATURATION 836 839 !------------------------------------ 837 840 838 841 CALL ice_sursat(pplay(i,k), paprs(i,k)-paprs(i,k+1), dtime, i, k, t(i,k), zq(i), & 839 gamma_ss(i,k), zqs(i), Tbef(i), rneb_seri(i,k), ratqs(i,k), & 840 rneb(i,k), zqn(i), rnebss(i,k), qss(i,k), & 842 gamma_ss(i,k), zqs(i), Tbef(i), rneb_seri(i,k), ratqs(i,k), & 843 rneb(i,k), zqn(i), rnebss(i,k), qss(i,k), & 841 844 Tcontr(i,k), qcontr(i,k), qcontr2(i,k), fcontrN(i,k), fcontrP(i,k) ) 842 845 … … 850 853 851 854 852 ! P2.2.2> Approximative calculation of temperature variation DT 853 ! due to condensation. 854 ! Calculated variables: 855 ! dT : temperature change due to condensation 856 ! --------------------------------------------------------------- 857 858 ! EV: calculation of icefrac in one sole function 859 CALL icefrac_lscp(klon, zt(:),pplay(:,k)/paprs(:,1),zfice(:),dzfice(:)) 855 ! P2.2.2> Approximative calculation of temperature variation DT 856 ! due to condensation. 857 ! Calculated variables: 858 ! dT : temperature change due to condensation 859 !--------------------------------------------------------------- 860 860 861 861 862 IF (zfice(i).LT.1) THEN … … 896 897 897 898 ! Partition function in stratiform clouds (will be overwritten in boundary-layer MPCs) 898 CALL icefrac_lscp(klon,zt(:),pplay(:,k)/paprs(:,1),zfice(:), dzfice(:))899 CALL ICEFRAC_LSCP(klon,zt(:),pplay(:,k)/paprs(:,1),zfice(:), dzfice(:)) 899 900 900 901 DO i=1, klon … … 1002 1003 1003 1004 DO i = 1, klon 1004 IF (rneb(i,k).GT.0.0) THEN1005 1005 zoliq(i) = zcond(i) 1006 zrho(i) = pplay(i,k) / zt(i) / RD 1007 zdz(i) = (paprs(i,k)-paprs(i,k+1)) / (zrho(i)*RG) 1008 1009 zneb(i) = MAX(rneb(i,k), seuil_neb) 1006 zoliqi(i)= zoliq(i)*zfice(i) 1007 zoliql(i)= zoliq(i)*(1.-zfice(i)) 1008 zrho(i) = pplay(i,k) / zt(i) / RD 1009 zdz(i) = (paprs(i,k)-paprs(i,k+1)) / (zrho(i)*RG) 1010 iwc(i) = 0. 1011 zneb(i) = MAX(rneb(i,k),seuil_neb) 1010 1012 radliq(i,k) = zoliq(i)/REAL(ninter+1) 1011 1013 radicefrac(i,k) = zfice(i) 1012 1014 radliqi(i,k)=zoliq(i)*zfice(i)/REAL(ninter+1) 1013 ENDIF1015 radliql(i,k)=zoliq(i)*(1.-zfice(i))/REAL(ninter+1) 1014 1016 ENDDO 1015 1017 … … 1024 1026 DO i=1,klon 1025 1027 IF (rneb(i,k).GT.0.0) THEN 1026 zrhol(i) = zrho(i) * zoliq(i) / zneb(i)1028 iwc(i) = zrho(i) * zoliqi(i) / zneb(i) ! in-cloud ice condensate content 1027 1029 ENDIF 1028 1030 ENDDO 1029 1031 1030 CALL FALLICE_VELOCITY(klon, zrhol(:),zt(:),zrho(:),pplay(:,k),ptconv(:,k),velo(:))1032 CALL FALLICE_VELOCITY(klon,iwc(:),zt(:),zrho(:),pplay(:,k),ptconv(:,k),velo(:,k)) 1031 1033 1032 1034 DO i = 1, klon … … 1034 1036 IF (rneb(i,k).GT.0.0) THEN 1035 1037 1036 ! Initialization of z pluie, zice and ztot:1037 z pluie=0.1038 z ice=0.1039 z tot=0.1038 ! Initialization of zrain, zsnow and zprecip: 1039 zrain=0. 1040 zsnow=0. 1041 zprecip=0. 1040 1042 1041 1043 IF (zneb(i).EQ.seuil_neb) THEN 1042 z tot= 0.01043 z ice= 0.01044 z pluie= 0.01044 zprecip = 0.0 1045 zsnow = 0.0 1046 zrain= 0.0 1045 1047 ELSE 1046 1048 ! water quantity to remove: zchau (Sundqvist, 1978) … … 1049 1051 zcl=cld_lc_con 1050 1052 zct=1./cld_tau_con 1051 zfroi = dtime/REAL(ninter)/zdz(i)*zoliq(i)*velo(i )*zfice(i)1053 zfroi = dtime/REAL(ninter)/zdz(i)*zoliq(i)*velo(i,k)*zfice(i) 1052 1054 ELSE 1053 1055 zcl=cld_lc_lsc 1054 1056 zct=1./cld_tau_lsc 1055 1057 zfroi = dtime/REAL(ninter)/zdz(i)*zoliq(i) & ! dqice/dt=1/rho*d(rho*wice*qice)/dz 1056 *velo(i ) * zfice(i)1058 *velo(i,k) * zfice(i) 1057 1059 ENDIF 1058 1060 … … 1061 1063 ! surface fraction (which is larger and artificially 1062 1064 ! reduces the in-cloud water). 1065 1063 1066 IF ((iflag_cloudth_vert.GE.3).AND.(iflag_rain_incloud_vol.EQ.1)) THEN 1064 1067 zchau = zct *dtime/REAL(ninter) * zoliq(i) & … … 1069 1072 ENDIF 1070 1073 1071 z pluie= MIN(MAX(zchau,0.0),zoliq(i)*(1.-zfice(i)))1072 z ice= MIN(MAX(zfroi,0.0),zoliq(i)*zfice(i))1073 z tot = MAX(zpluie + zice,0.0)1074 zrain = MIN(MAX(zchau,0.0),zoliq(i)*(1.-zfice(i))) 1075 zsnow = MIN(MAX(zfroi,0.0),zoliq(i)*zfice(i)) 1076 zprecip = MAX(zrain + zsnow,0.0) 1074 1077 1075 1078 ENDIF 1076 1079 1077 ztot = MIN(ztot,zoliq(i)) 1078 zice = MIN(zice,ztot) 1079 zpluie = MIN(zpluie,ztot) 1080 1081 zoliqp(i) = MAX(zoliq(i)*(1.-zfice(i))-1.*zpluie , 0.0) 1082 zoliqi(i) = MAX(zoliq(i)*zfice(i)-1.*zice , 0.0) 1083 zoliq(i) = MAX(zoliq(i)-ztot , 0.0) 1080 zoliql(i) = MAX(zoliq(i)*(1.-zfice(i))-1.*zrain , 0.0) 1081 zoliqi(i) = MAX(zoliq(i)*zfice(i)-1.*zsnow , 0.0) 1082 zoliq(i) = MAX(zoliq(i)-zprecip , 0.0) 1084 1083 1085 1084 radliq(i,k) = radliq(i,k) + zoliq(i)/REAL(ninter+1) 1085 radliql(i,k) = radliql(i,k) + zoliql(i)/REAL(ninter+1) 1086 1086 radliqi(i,k) = radliqi(i,k) + zoliqi(i)/REAL(ninter+1) 1087 1087 1088 ENDIF 1088 ENDIF ! rneb >0 1089 1089 1090 1090 ENDDO ! i = 1,klon … … 1093 1093 1094 1094 1095 ! Caculate the percentage of ice in "radliq" so cloud+precip seen by the radiation scheme 1095 ! Include the contribution to non-evaporated precipitation to radliq 1096 IF ((ok_radliq_precip) .AND. (k .LT. klev)) THEN 1097 ! for rain drops, we assume a fallspeed of 5m/s 1098 vr=5.0 1099 radliql(:,k)=radliql(:,k)+zrfl(:)/vr/zrho(:) 1100 ! for ice crystals, we take the fallspeed from level above 1101 radliqi(:,k)=radliqi(:,k)+zifl(:)/zrho(:)/velo(:,k+1) 1102 radliq(:,k)=radliql(:,k)+radliqi(:,k) 1103 ENDIF 1104 1105 ! caculate the percentage of ice in "radliq" so cloud+precip seen by the radiation scheme 1096 1106 DO i=1,klon 1097 IF (radliq(i,k) .GT. 0 ) THEN1107 IF (radliq(i,k) .GT. 0.) THEN 1098 1108 radicefrac(i,k)=MIN(MAX(radliqi(i,k)/radliq(i,k),0.),1.) 1099 1109 ENDIF … … 1101 1111 1102 1112 1103 ! CR&JYG: We account for the Wegener-Findeisen-Bergeron process in the precipitation flux: 1104 ! If T<0C, liquid precip are converted into ice, which leads to an increase in 1105 ! temperature DeltaT. The effect of DeltaT on condensates and precipitation is roughly 1106 ! taken into account through a linearization of the equations and by approximating 1107 ! the liquid precipitation process with a "threshold" process. We assume that 1108 ! condensates are not modified during this operation. Liquid precipitation is 1109 ! removed (in the limit DeltaT<273.15-T). Solid precipitation increases. 1110 ! Water vapor increases as well 1111 ! Note that compared to fisrtilp, we always assume iflag_bergeron=2 1112 1113 1113 ! Precipitation flux calculation 1114 1114 1115 1115 DO i = 1, klon 1116 1116 1117 1117 IF (rneb(i,k) .GT. 0.0) THEN 1118 1119 zqpreci(i)=(zcond(i)-zoliq(i))*zfice(i) 1120 zqprecl(i)=(zcond(i)-zoliq(i))*(1.-zfice(i)) 1121 zcp=RCPD*(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i))) 1122 coef1 = rneb(i,k)*RLSTT/zcp*zdqsdT_raw(i) 1123 ! Computation of DT if all the liquid precip freezes 1124 DeltaT = RLMLT*zqprecl(i) / (zcp*(1.+coef1)) 1125 ! T should not exceed the freezing point 1126 ! that is Delta > RTT-zt(i) 1127 DeltaT = max( min( RTT-zt(i), DeltaT) , 0. ) 1128 zt(i) = zt(i) + DeltaT 1129 ! water vaporization due to temp. increase 1130 Deltaq = rneb(i,k)*zdqsdT_raw(i)*DeltaT 1131 ! we add this vaporized water to the total vapor and we remove it from the precip 1132 zq(i) = zq(i) + Deltaq 1133 ! The three "max" lines herebelow protect from rounding errors 1134 zcond(i) = max( zcond(i)- Deltaq, 0. ) 1135 ! liquid precipitation freezes oe evaporates 1136 Deltaqprecl = -zcp/RLMLT*(1.+coef1)*DeltaT 1137 zqprecl(i) = max( zqprecl(i) + Deltaqprecl, 0. ) 1138 ! iced water budget 1139 zqpreci(i) = max (zqpreci(i) - Deltaqprecl - Deltaq, 0.) 1118 1119 1120 ! CR&JYG: We account for the Wegener-Findeisen-Bergeron process in the precipitation flux: 1121 ! If T<0C, liquid precip are converted into ice, which leads to an increase in 1122 ! temperature DeltaT. The effect of DeltaT on condensates and precipitation is roughly 1123 ! taken into account through a linearization of the equations and by approximating 1124 ! the liquid precipitation process with a "threshold" process. We assume that 1125 ! condensates are not modified during this operation. Liquid precipitation is 1126 ! removed (in the limit DeltaT<273.15-T). Solid precipitation increases. 1127 ! Water vapor increases as well 1128 ! Note that compared to fisrtilp, we always assume iflag_bergeron=2 1129 1130 zqpreci(i)=(zcond(i)-zoliq(i))*zfice(i) 1131 zqprecl(i)=(zcond(i)-zoliq(i))*(1.-zfice(i)) 1132 zcp=RCPD*(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i))) 1133 coef1 = rneb(i,k)*RLSTT/zcp*zdqsdT_raw(i) 1134 ! Computation of DT if all the liquid precip freezes 1135 DeltaT = RLMLT*zqprecl(i) / (zcp*(1.+coef1)) 1136 ! T should not exceed the freezing point 1137 ! that is Delta > RTT-zt(i) 1138 DeltaT = max( min( RTT-zt(i), DeltaT) , 0. ) 1139 zt(i) = zt(i) + DeltaT 1140 ! water vaporization due to temp. increase 1141 Deltaq = rneb(i,k)*zdqsdT_raw(i)*DeltaT 1142 ! we add this vaporized water to the total vapor and we remove it from the precip 1143 zq(i) = zq(i) + Deltaq 1144 ! The three "max" lines herebelow protect from rounding errors 1145 zcond(i) = max( zcond(i)- Deltaq, 0. ) 1146 ! liquid precipitation converted to ice precip 1147 Deltaqprecl = -zcp/RLMLT*(1.+coef1)*DeltaT 1148 zqprecl(i) = max( zqprecl(i) + Deltaqprecl, 0. ) 1149 ! iced water budget 1150 zqpreci(i) = max (zqpreci(i) - Deltaqprecl - Deltaq, 0.) 1151 1140 1152 1141 d_ql(i,k) = (1-zfice(i))*zoliq(i)1142 d_qi(i,k) = zfice(i)*zoliq(i)1143 1144 1153 IF (iflag_evap_prec.EQ.4) THEN 1145 1154 zrflcld(i) = zrflcld(i)+zqprecl(i) & … … 1165 1174 IF (iflag_evap_prec.EQ.4) THEN 1166 1175 1167 DO i=1, 1176 DO i=1,klon 1168 1177 1169 1178 IF ((zrflclr(i) + ziflclr(i)) .GT. 0. ) THEN … … 1192 1201 ! Outputs: 1193 1202 ! Precipitation fluxes at layer interfaces 1194 ! and temperature and vaportendencies1203 ! and temperature and water species tendencies 1195 1204 DO i = 1, klon 1196 1205 psfl(i,k)=zifl(i) 1197 1206 prfl(i,k)=zrfl(i) 1207 d_ql(i,k) = (1-zfice(i))*zoliq(i) 1208 d_qi(i,k) = zfice(i)*zoliq(i) 1198 1209 d_q(i,k) = zq(i) - q(i,k) 1199 1210 d_t(i,k) = zt(i) - t(i,k) … … 1264 1275 ENDDO 1265 1276 1266 !--save some variables for ice su rsaturation1277 !--save some variables for ice supersaturation 1267 1278 ! 1268 1279 DO i = 1, klon 1269 ! pour la mémoire1280 ! for memory 1270 1281 rneb_seri(i,k) = rneb(i,k) 1271 1282 1272 ! pour lesdiagnostics1283 ! for diagnostics 1273 1284 rnebclr(i,k) = 1.0 - rneb(i,k) - rnebss(i,k) 1274 1285 1275 1286 qvc(i,k) = zqs(i) * rneb(i,k) 1276 qclr(i,k) = MAX(1.e-10,zq(i) - qvc(i,k) - qss(i,k)) !--a jout OB a cause de cas pathologiques avec lognormale=F1287 qclr(i,k) = MAX(1.e-10,zq(i) - qvc(i,k) - qss(i,k)) !--added by OB because of pathologic cases with the lognormal 1277 1288 qcld(i,k) = qvc(i,k) + zcond(i) 1278 1279 !q_sat 1280 CALL CALC_QSAT_ECMWF(Tbef(i),0.,pplay(i,k),RTT,1,.false.,zqsatl(i,k),zdqs(i)) 1281 CALL CALC_QSAT_ECMWF(Tbef(i),0.,pplay(i,k),RTT,2,.false.,zqsats(i,k),zdqs(i)) 1282 1283 ENDDO 1284 1285 ENDDO 1289 END DO 1290 !q_sat 1291 CALL CALC_QSAT_ECMWF(klon,Tbef(:),qzero(:),pplay(:,k),RTT,1,.false.,zqsatl(:,k),zdqs(:)) 1292 CALL CALC_QSAT_ECMWF(klon,Tbef(:),qzero(:),pplay(:,k),RTT,2,.false.,zqsats(:,k),zdqs(:)) 1293 1294 END DO 1286 1295 1287 1296 !====================================================================== -
LMDZ6/trunk/libf/phylmd/lscp_tools_mod.F90
r4059 r4072 47 47 48 48 tempc=temp(i)-273.15 ! celcius temp 49 iwcg= iwc(i)*1000. ! iwc in g/m349 iwcg=MAX(iwc(i)*1000.,1E-3) ! iwc in g/m3. We set a minimum value to prevent from division by 0 50 50 phpa=pres(i)/100. ! pressure in hPa 51 51 … … 76 76 **(1./(0.807+bice*0.00581+0.0457*bice**2)) 77 77 78 vmice=100.*1042.4*exp(13.6-(bice+1)*7.76+0.479*(bice+1.)**2)*exp((-0.0361+ &79 (bice+1.)*0.0151+0.00149*(bice+1.)**2)*tempc) &80 *(m2ice**(0.807+(bice+1.)*0.00581+0.0457*(bice+1.)**2))/(iwcg*0.001/aice)78 vmice=100.*1042.4*exp(13.6-(bice+1)*7.76+0.479*(bice+1.)**2)*exp((-0.0361+ & 79 (bice+1.)*0.0151+0.00149*(bice+1.)**2)*tempc) & 80 *(m2ice**(0.807+(bice+1.)*0.00581+0.0457*(bice+1.)**2))/(iwcg*0.001/aice) 81 81 82 82 … … 93 93 94 94 vmsnow=vmsnow*((1000./phpa)**0.35) 95 96 95 velo(i)=fallv_tun*min(vmsnow,vmice)/100. ! to m/s 97 96 dvel=0.2 98 cvel=velo(i)/((iwc (i)*rho(i))**dvel)97 cvel=velo(i)/((iwcg/1000.*rho(i))**dvel) 99 98 100 99 ELSE 101 100 ! By default, fallspeed velocity of ice crystals according to Heymsfield & Donner 1990 102 velo(i) = fallv_tun*3.29/2.0 * ((iwc(i))**0.16)101 velo(i) = fallv_tun*3.29/2.0*((iwcg/1000.)**0.16) 103 102 dvel=0.16 104 103 cvel=fallv_tun*3.29/2.0*(rho(i)**0.16) 105 104 ENDIF 106 107 105 ENDDO 108 106 … … 202 200 ENDDO 203 201 202 204 203 RETURN 205 204 … … 209 208 210 209 211 SUBROUTINE CALC_QSAT_ECMWF( temp,qtot,pressure,tref,phase,flagth,qs,dqs)210 SUBROUTINE CALC_QSAT_ECMWF(klon,temp,qtot,pressure,tref,phase,flagth,qs,dqs) 212 211 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 213 212 ! Calculate qsat following ECMWF method 214 213 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 214 215 215 216 216 IMPLICIT NONE … … 220 220 include "FCTTRE.h" 221 221 222 REAL, INTENT(IN) :: temp ! temperature in K 223 REAL, INTENT(IN) :: qtot ! total specific water in kg/kg 224 REAL, INTENT(IN) :: pressure ! pressure in Pa 225 REAL, INTENT(IN) :: tref ! reference temperature in K 222 INTEGER, INTENT(IN) :: klon ! number of horizontal grid points 223 REAL, INTENT(IN), DIMENSION(klon) :: temp ! temperature in K 224 REAL, INTENT(IN), DIMENSION(klon) :: qtot ! total specific water in kg/kg 225 REAL, INTENT(IN), DIMENSION(klon) :: pressure ! pressure in Pa 226 REAL, INTENT(IN) :: tref ! reference temperature in K 226 227 LOGICAL, INTENT(IN) :: flagth ! flag for qsat calculation for thermals 227 228 228 INTEGER, INTENT(IN) :: phase 229 229 ! phase: 0=depend on temperature sign (temp>tref -> liquid, temp<tref, solid) … … 231 231 ! 2=solid 232 232 233 REAL, INTENT(OUT) :: qs ! saturation specific humidity [kg/kg]234 REAL, INTENT(OUT) :: dqs ! derivation of saturation specific humidity wrt T233 REAL, INTENT(OUT), DIMENSION(klon) :: qs ! saturation specific humidity [kg/kg] 234 REAL, INTENT(OUT), DIMENSION(klon) :: dqs ! derivation of saturation specific humidity wrt T 235 235 236 236 REAL delta, cor, cvm5 237 237 INTEGER i 238 239 DO i=1,klon 240 238 241 IF (phase .EQ. 1) THEN 239 242 delta=0. … … 241 244 delta=1. 242 245 ELSE 243 delta=MAX(0.,SIGN(1.,tref-temp ))246 delta=MAX(0.,SIGN(1.,tref-temp(i))) 244 247 ENDIF 245 248 … … 248 251 ELSE 249 252 cvm5 = R5LES*RLVTT*(1.-delta) + R5IES*RLSTT*delta 250 cvm5 = cvm5 /RCPD/(1.0+RVTMP2*(qtot)) 251 ENDIF 252 253 qs= R2ES*FOEEW(temp,delta)/pressure 254 qs=MIN(0.5,qs) 255 cor=1./(1.-RETV*qs) 256 qs=qs*cor 257 dqs= FOEDE(temp,delta,cvm5,qs,cor) 253 cvm5 = cvm5 /RCPD/(1.0+RVTMP2*(qtot(i))) 254 ENDIF 255 256 qs(i)= R2ES*FOEEW(temp(i),delta)/pressure(i) 257 qs(i)=MIN(0.5,qs(i)) 258 cor=1./(1.-RETV*qs(i)) 259 qs(i)=qs(i)*cor 260 dqs(i)= FOEDE(temp(i),delta,cvm5,qs(i),cor) 261 262 END DO 258 263 259 264 END SUBROUTINE CALC_QSAT_ECMWF … … 262 267 263 268 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 264 SUBROUTINE CALC_GAMMASAT( temp,qtot,pressure,gammasat,dgammasatdt)269 SUBROUTINE CALC_GAMMASAT(klon,temp,qtot,pressure,gammasat,dgammasatdt) 265 270 266 271 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ … … 279 284 include "nuage.h" 280 285 281 282 REAL, INTENT(IN) :: temp ! temperature in K 283 REAL, INTENT(IN) :: qtot ! total specific water in kg/kg 284 285 REAL, INTENT(IN) :: pressure ! pressure in Pa 286 287 REAL, INTENT(OUT) :: gammasat ! coefficient to multiply qsat with to calculate saturation 288 REAL, INTENT(OUT) :: dgammasatdt ! derivative of gammasat wrt temperature 289 290 REAL qsi,qsl,fac,dqsl,dqsi,fcirrus 286 INTEGER, INTENT(IN) :: klon ! number of horizontal grid points 287 REAL, INTENT(IN), DIMENSION(klon) :: temp ! temperature in K 288 REAL, INTENT(IN), DIMENSION(klon) :: qtot ! total specific water in kg/kg 289 290 REAL, INTENT(IN), DIMENSION(klon) :: pressure ! pressure in Pa 291 292 REAL, INTENT(OUT), DIMENSION(klon) :: gammasat ! coefficient to multiply qsat with to calculate saturation 293 REAL, INTENT(OUT), DIMENSION(klon) :: dgammasatdt ! derivative of gammasat wrt temperature 294 295 REAL, DIMENSION(klon) :: qsi,qsl,dqsl,dqsi 296 REAL fcirrus, fac 291 297 REAL, PARAMETER :: acirrus=2.349 292 298 REAL, PARAMETER :: bcirrus=259.0 293 299 294 295 CALL CALC_QSAT_ECMWF(temp,qtot,pressure,RTT,1,.false.,qsl,dqsl) 296 CALL CALC_QSAT_ECMWF(temp,qtot,pressure,RTT,2,.false.,qsi,dqsi) 297 298 IF (temp .GE. RTT) THEN 300 INTEGER i 301 302 CALL CALC_QSAT_ECMWF(klon,temp,qtot,pressure,RTT,1,.false.,qsl,dqsl) 303 CALL CALC_QSAT_ECMWF(klon,temp,qtot,pressure,RTT,2,.false.,qsi,dqsi) 304 305 DO i=1,klon 306 307 IF (temp(i) .GE. RTT) THEN 299 308 ! warm clouds: condensation at saturation wrt liquid 300 gammasat =1.301 dgammasatdt =0.302 303 ELSEIF ((temp .LT. RTT) .AND. (temp.GT. t_glace_min)) THEN309 gammasat(i)=1. 310 dgammasatdt(i)=0. 311 312 ELSEIF ((temp(i) .LT. RTT) .AND. (temp(i) .GT. t_glace_min)) THEN 304 313 305 314 IF (iflag_gammasat .GE. 2) THEN 306 gammasat =qsl/qsi307 dgammasatdt =(dqsl*qsi-dqsi*qsl)/qsi/qsi315 gammasat(i)=qsl(i)/qsi(i) 316 dgammasatdt(i)=(dqsl(i)*qsi(i)-dqsi(i)*qsl(i))/qsi(i)/qsi(i) 308 317 ELSE 309 gammasat =1.310 dgammasatdt =0.318 gammasat(i)=1. 319 dgammasatdt(i)=0. 311 320 ENDIF 312 321 … … 317 326 ! Koop, 2000 and Karcher 2008, QJRMS 318 327 ! 'Cirrus regime' 319 fcirrus=acirrus-temp /bcirrus320 IF (fcirrus .LT. qsl /qsi) THEN321 gammasat =qsl/qsi322 dgammasatdt =(dqsl*qsi-dqsi*qsl)/qsi/qsi328 fcirrus=acirrus-temp(i)/bcirrus 329 IF (fcirrus .LT. qsl(i)/qsi(i)) THEN 330 gammasat(i)=qsl(i)/qsi(i) 331 dgammasatdt(i)=(dqsl(i)*qsi(i)-dqsi(i)*qsl(i))/qsi(i)/qsi(i) 323 332 ELSE 324 gammasat =fcirrus325 dgammasatdt =-1.0/bcirrus333 gammasat(i)=fcirrus 334 dgammasatdt(i)=-1.0/bcirrus 326 335 ENDIF 327 336 328 337 ELSE 329 338 330 gammasat =1.331 dgammasatdt =0.339 gammasat(i)=1. 340 dgammasatdt(i)=0. 332 341 333 342 ENDIF … … 335 344 ENDIF 336 345 346 END DO 347 348 337 349 END SUBROUTINE CALC_GAMMASAT 338 350 -
LMDZ6/trunk/libf/phylmd/phys_local_var_mod.F90
r4059 r4072 438 438 REAL,ALLOCATABLE,SAVE,DIMENSION(:,:) :: cldemi, cldfra, cldtau, fiwc, fl, re, flwc 439 439 !$OMP THREADPRIVATE(cldemi, cldfra, cldtau, fiwc, fl, re, flwc) 440 REAL,ALLOCATABLE,SAVE,DIMENSION(:,:) :: qlth, qith 441 !$OMP THREADPRIVATE(qlth, qith )440 REAL,ALLOCATABLE,SAVE,DIMENSION(:,:) :: qlth, qith, qsith 441 !$OMP THREADPRIVATE(qlth, qith, qsith) 442 442 REAL,ALLOCATABLE,SAVE,DIMENSION(:,:) :: ref_liq, ref_ice, theta, zphi 443 443 !$OMP THREADPRIVATE(ref_liq, ref_ice, theta, zphi) … … 458 458 REAL,ALLOCATABLE,SAVE,DIMENSION(:,:,:) :: phiwriteSTD, qwriteSTD, twriteSTD, rhwriteSTD 459 459 !$OMP THREADPRIVATE(phiwriteSTD, qwriteSTD, twriteSTD, rhwriteSTD) 460 460 461 461 462 ! ug et d'autres encore: … … 777 778 ALLOCATE(rain_lsc(klon)) 778 779 ALLOCATE(rain_num(klon)) 779 ALLOCATE(qlth(klon,klev), qith(klon,klev) )780 ALLOCATE(qlth(klon,klev), qith(klon,klev), qsith(klon,klev)) 780 781 ! 781 782 ALLOCATE(sens_x(klon), sens_w(klon)) … … 823 824 ALLOCATE(ref_liq_pi(klon, klev), ref_ice_pi(klon, klev)) 824 825 ALLOCATE(zphi(klon, klev), zx_rh(klon, klev), zx_rhl(klon,klev), zx_rhi(klon,klev)) 826 zx_rhl(:,:)=0.; zx_rhi(:,:)=0. ! because not always defined 825 827 ALLOCATE(pmfd(klon, klev), pmfu(klon, klev)) 826 828 … … 1087 1089 DEALLOCATE(rain_lsc) 1088 1090 DEALLOCATE(rain_num) 1089 DEALLOCATE(qlth, qith )1091 DEALLOCATE(qlth, qith, qsith) 1090 1092 ! 1091 1093 DEALLOCATE(sens_x, sens_w) -
LMDZ6/trunk/libf/phylmdiso/phys_local_var_mod.F90
r4065 r4072 559 559 REAL,ALLOCATABLE,SAVE,DIMENSION(:,:) :: cldemi, cldfra, cldtau, fiwc, fl, re, flwc 560 560 !$OMP THREADPRIVATE(cldemi, cldfra, cldtau, fiwc, fl, re, flwc) 561 REAL,ALLOCATABLE,SAVE,DIMENSION(:,:) :: qlth, qith 562 !$OMP THREADPRIVATE(qlth, qith )561 REAL,ALLOCATABLE,SAVE,DIMENSION(:,:) :: qlth, qith, qsith 562 !$OMP THREADPRIVATE(qlth, qith, qsith) 563 563 REAL,ALLOCATABLE,SAVE,DIMENSION(:,:) :: ref_liq, ref_ice, theta, zphi 564 564 !$OMP THREADPRIVATE(ref_liq, ref_ice, theta, zphi) … … 943 943 ALLOCATE(rain_lsc(klon)) 944 944 ALLOCATE(rain_num(klon)) 945 ALLOCATE(qlth(klon,klev), qith(klon,klev) )945 ALLOCATE(qlth(klon,klev), qith(klon,klev), qsith(klon,klev)) 946 946 ! 947 947 #ifdef ISO … … 1004 1004 ALLOCATE(ref_liq_pi(klon, klev), ref_ice_pi(klon, klev)) 1005 1005 ALLOCATE(zphi(klon, klev), zx_rh(klon, klev), zx_rhl(klon,klev), zx_rhi(klon,klev)) 1006 zx_rhl(:,:)=0.; zx_rhi(:,:)=0. ! because not always defined 1006 1007 ALLOCATE(pmfd(klon, klev), pmfu(klon, klev)) 1007 1008 … … 1324 1325 DEALLOCATE(rain_lsc) 1325 1326 DEALLOCATE(rain_num) 1326 DEALLOCATE(qlth, qith )1327 DEALLOCATE(qlth, qith, qsith) 1327 1328 ! 1328 1329 DEALLOCATE(sens_x, sens_w)
Note: See TracChangeset
for help on using the changeset viewer.