Changeset 4072


Ignore:
Timestamp:
Feb 1, 2022, 6:34:34 PM (2 years ago)
Author:
evignon
Message:

Mise à jour de la routine de nuages lscp
les principaux changements consistent en:

  • des corrections de bug (déclaration et division

par zero) pour que la routine compile avec les modifications d'OB

  • des restructurations pour que les fonctions de lscp_tools_mod

travaillent sur des tableaux et non des scalaires

  • une séparation des coefficients d'évaporation de la neige et de la pluie
Location:
LMDZ6/trunk/libf
Files:
8 edited

Legend:

Unmodified
Added
Removed

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

    r4010 r4072  
    230230
    231231INTEGER :: i,k,nsrf
    232 REAL :: ratiom, qsat2m, dqsatdT
    233 REAL, DIMENSION(klon) :: xsi0
     232REAL, DIMENSION(klon) :: xsi0, ratiom, qsat2m, dqsatdT
    234233REAL, DIMENSION (klon,klev) :: zlay
    235234
    … …  
    239238!-------------------------------------------
    240239
    241 DO i=1,klon
    242240   
    243     ratiom=0.
    244     xsi0(i)=0.
     241    ratiom(:)=0.
     242    xsi0(:)=0.
    245243   
    246244    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)
    250248    END DO
    251249   
    252     xsi0(i)=sqrt(xsi0(i))/(ratiom+1E-6)
    253 END DO
     250    xsi0(:)=sqrt(xsi0(:))/(ratiom(:)+1E-6)
    254251
    255252

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

    r3999 r4072  
    15391539
    15401540!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    1541 SUBROUTINE cloudth_mpc(klon,klev,ind2,mpc_bl_points,                     &
     1541SUBROUTINE cloudth_mpc(klon,klev,ind2,mpc_bl_points,                        &
    15421542&           temp,ztv,po,zqta,fraca,zpspsk,paprs,pplay,ztla,zthl,            &
    15431543&           ratqs,zqs,snowflux,qcloud,qincloud,icefrac,ctot,ctot_vol)
    … …  
    15541554      USE ioipsl_getin_p_mod, ONLY : getin_p
    15551555      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
     1556      USE lscp_tools_mod, ONLY: CALC_QSAT_ECMWF
     1557      USE phys_local_var_mod, ONLY : qlth, qith, qsith
    15581558
    15591559      IMPLICIT NONE
    … …  
    15731573
    15741574      INTEGER, INTENT(IN)                         :: klon,klev,ind2
    1575       INTEGER, DIMENSION(klon,klev), INTENT(INOUT)   ::  mpc_bl_points ! 1 where BL MPC, 0 otherwise
     1575     
    15761576
    15771577      REAL, DIMENSION(klon,klev), INTENT(IN)      ::  temp          ! Temperature [K]
    … …  
    15891589      REAL, DIMENSION(klon,klev+1), INTENT(IN)    ::  snowflux      ! snow flux at the interface of the layer [kg/m2/s]
    15901590
    1591 
     1591      INTEGER, DIMENSION(klon,klev), INTENT(INOUT) :: mpc_bl_points  ! grid points with convective (thermals) mixed phase clouds
     1592     
    15921593      REAL, DIMENSION(klon,klev), INTENT(OUT)      ::  ctot         ! Cloud fraction [0-1]
    15931594      REAL, DIMENSION(klon,klev), INTENT(OUT)      ::  ctot_vol     ! Volume cloud fraction [0-1]
    … …  
    16001601
    16011602      INTEGER itap,ind1,l,ig,iter,k
    1602       LOGICAL flag_topthermals
    1603 
    1604 
    1605       REAL zqsatth(klon,klev), zqsatenv(klon,klev)
     1603      INTEGER iflag_topthermals, niter
     1604
     1605
     1606      REAL zqsatth(klon), zqsatenv(klon), zqsatenvonly(klon)
    16061607      REAL sigma1(klon,klev)                                                         
    16071608      REAL sigma2(klon,klev)
    … …  
    16141615      REAL cenv_vol(klon,klev)
    16151616      REAL rneb(klon,klev)
    1616       REAL zqenv(klon)   
    1617 
     1617      REAL zqenv(klon), zqth(klon), zqenvonly(klon)
    16181618      REAL qsatmmussig1,qsatmmussig2,sqrtpi,sqrt2,sqrt2pi,pi
    16191619      REAL rdd,cppd,Lv
    … …  
    16241624      REAL xth1,xth2,xenv1,xenv2,deltasth, deltasenv
    16251625      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)
    16281630      REAL erf
    16291631      REAL zpdf_sig(klon),zpdf_k(klon),zpdf_delta(klon)
    … …  
    16321634      REAL zrho(klon,klev)
    16331635      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
    16361639      REAL senvi, senvl, qbase, sbase, qliqth, qiceth
    1637       REAL qimax, ttarget, stmp, cout, coutref
     1640      REAL qmax, ttarget, stmp, cout, coutref, fraci
    16381641      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)
    16421643
    16431644      ! Modifty the saturation deficit PDF in thermals
    … …  
    16451646      REAL,SAVE :: C_mpc
    16461647      !$OMP THREADPRIVATE(C_mpc)
     1648      REAL, SAVE    :: Ni,C_cap,Ei,d_top
     1649      !$OMP THREADPRIVATE(Ni,C_cap,Ei,d_top)
    16471650      ! Change the width of the PDF used for vertical subgrid scale heterogeneity
    16481651      ! (J Jouhaud, JL Dufresne, JB Madeleine)
    … …  
    17091712      sqrtpi=sqrt(pi)
    17101713      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.
    17111723
    17121724
    … …  
    17411753        WRITE(*,*) 'iflag_cloudth_vert_noratqs = ', iflag_cloudth_vert_noratqs
    17421754        ! Modifies the PDF in thermals when ice crystals are present
    1743         C_mpc=1.e2
     1755        C_mpc=1.e4
    17441756        CALL getin_p('C_mpc',C_mpc)
    17451757        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
    17461770
    17471771        firstcall=.FALSE.
    … …  
    17591783      DO ind1=1,klon
    17601784            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)  &
    17621786            .AND. (ztv(ind1,1).GT.ztv(ind1,2)) .AND.(fraca(ind1,ind2).GT.1.e-10)) THEN
    17631787                mpc_bl_points(ind1,ind2)=1
    … …  
    17721796!------------------------------------------------------------------------------- 
    17731797
    1774     DO ind1=1,klon
     1798! calculation of temperature, humidity and saturation specific humidity
     1799
     1800Tbefenv(:)=zthl(:,ind2)*zpspsk(:,ind2)
     1801Tbefth(:)=ztla(:,ind2)*zpspsk(:,ind2)
     1802Tbefenvonly(:)=temp(:,ind2)
     1803
     1804zqenv(:)=(po(:)-fraca(:,ind2)*zqta(:,ind2))/(1.-fraca(:,ind2)) !qt = a*qtth + (1-a)*qtenv
     1805zqth(:)=zqta(:,ind2)
     1806zqenvonly(:)=po(:)
     1807
     1808
     1809CALL CALC_QSAT_ECMWF(klon,Tbefenvonly,zqenvonly,paprs(:,ind2),RTT,0,.false.,zqsatenvonly,dqsatenv)
     1810
     1811CALL CALC_QSAT_ECMWF(klon,Tbefenv,zqenv,paprs(:,ind2),RTT,0,.false.,zqsatenv,dqsatenv)
     1812CALL CALC_QSAT_ECMWF(klon,Tbefenv,zqenv,paprs(:,ind2),RTT,1,.false.,qslenv,dqsatenv)
     1813
     1814CALL CALC_QSAT_ECMWF(klon,Tbefth,zqth,paprs(:,ind2),RTT,1,.false.,qslth,dqsatth)
     1815CALL CALC_QSAT_ECMWF(klon,Tbefth,zqth,paprs(:,ind2),RTT,2,.false.,qsith(:,ind2),dqsatth)
     1816CALL CALC_QSAT_ECMWF(klon,Tbefth,zqth,paprs(:,ind2),RTT,0,.false.,zqsatth,dqsatth)
     1817
     1818
     1819  DO ind1=1,klon
    17751820
    17761821
    … …  
    17801825! Environment:
    17811826
    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
    17891829               Lv=RLVTT
    17901830        ELSE
    … …  
    17931833       
    17941834
    1795         alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)     !qsl, p84
     1835        alenv=(0.622*Lv*zqsatenv(ind1))/(rdd*zthl(ind1,ind2)**2)     !qsl, p84
    17961836        aenv=1./(1.+(alenv*Lv/cppd))                                      !al, p84
    1797         senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))                          !s, p84
     1837        senv=aenv*(po(ind1)-zqsatenv(ind1))                          !s, p84
    17981838     
    17991839        ! For MPCs:
    18001840        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)     
    18031842        aenvl=1./(1.+(alenv*Lv/cppd))         
    1804         senvl=aenvl*(po(ind1)-qslenv)   
     1843        senvl=aenvl*(po(ind1)-qslenv(ind1))   
     1844        senvi=senv
    18051845        ENDIF
    18061846
    … …  
    18081848! Thermals:
    18091849
    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
    18151852            Lv=RLVTT
    18161853        ELSE
    … …  
    18191856
    18201857       
    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)       
    18221859        ath=1./(1.+(alth*Lv/cppd))                                                         
    1823         sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2))                     
     1860        sth=ath*(zqta(ind1,ind2)-zqsatth(ind1))                     
    18241861
    18251862       ! For MPCs:
    18261863        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))
    18331871        ENDIF     
    18341872
    … …  
    19231961                ctot(ind1,ind2)=0.
    19241962                ctot_vol(ind1,ind2)=0.
    1925                 qcloud(ind1)=zqsatenv(ind1,ind2)
     1963                qcloud(ind1)=zqsatenv(ind1)
    19261964                qincloud(ind1)=0.
    19271965            ELSE               
    … …  
    19411979            !...........
    19421980
     1981            qiceth=0.
    19431982            deltazlev_mpc=dz(ind1,:)
    19441983            temp_mpc=ztla(ind1,:)*zpspsk(ind1,:)
    … …  
    19461985            fraca_mpc=fraca(ind1,:)
    19471986            snowf_mpc=snowflux(ind1,:)
    1948             qth_mpc=zqta(ind1,:)
    1949             flag_topthermals=.FALSE.
     1987            iflag_topthermals=0
    19501988            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
    19521995            ENDIF
    19531996
    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)
    19712001            deltasth=athl*vert_alpha_th*sigma2s
    1972      
    1973             ! Liquid phase
    1974             !.............
    19752002            xth1=-(sthl+deltasth)/(sqrt(2.)*sigma2s)
    1976             xth2=-(sthl-deltasth)/(sqrt(2.)*sigma2s)           
     2003            xth2=-(sthl-deltasth)/(sqrt(2.)*sigma2s)
    19772004            exp_xth1 = exp(-1.*xth1**2)
    19782005            exp_xth2 = exp(-1.*xth2**2)
    … …  
    19842011            IntI1_CF=((sthl+deltasth)*(erf(xth2)-erf(xth1)))/(4*deltasth)
    19852012            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)           
    19912030            exp_xth1 = exp(-1.*xth1**2)
    19922031            exp_xth2 = exp(-1.*xth2**2)
    1993             IntJ=0.5*sthi*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2
     2032            IntJ=0.5*sthlc*(1-erf(xth2))+(sigma2sc/sqrt2pi)*exp_xth2
    19942033            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 
     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)
    20032042
    20042043            ! 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           
    20082045            qcth(ind1,ind2)=qlth(ind1,ind2)+qith(ind1,ind2)
    20092046
     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            !...............
    20102059            ! calculation of qbase which is the value of the water vapor within mixed phase clouds
    20112060            ! such that the total water in cloud = qbase+qliqth+qiceth
    … …  
    20152064           
    20162065            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
    20192069            pas=(maxi-mini)/niter
    20202070            stmp=mini
    … …  
    20312081                ENDIF
    20322082            ENDDO
    2033             qbase=MAX(0., sbase/athl+qslth)
     2083            qbase=MAX(0., sbase/athl+qslth(ind1))
    20342084
    20352085            ! surface cloud fraction in thermals
    … …  
    20532103              IntI1=(((sthl+deltasth-sbase)**2+(sigma2s)**2)/(8*deltasth))*(erf(xth2)-erf(xth1))
    20542104              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)
    20562106              IntI1_CF=((sthl-sbase+deltasth)*(erf(xth2)-erf(xth1)))/(4*deltasth)
    20572107              IntI3_CF=(sqrt2*sigma2s*(exp_xth1-exp_xth2))/(4*sqrtpi*deltasth)
    … …  
    20592109            ENDIF
    20602110              cth_vol(ind1,ind2)=MIN(MAX(0.,cth_vol(ind1,ind2)),1.)
     2111
     2112
    20612113
    20622114            ! Environment
    … …  
    21102162
    21112163           
    2112             ! Thermals + environment
     2164            ! Thermals + environment
     2165            ! =======================
    21132166            ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)
    21142167            qctot(ind1,ind2)=fraca(ind1,ind2)*qcth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qcenv(ind1,ind2)
    … …  
    21252178                ctot_vol(ind1,ind2)=0.
    21262179                qincloud(ind1)=0.
    2127                 qcloud(ind1)=zqsatenv(ind1,ind2)
     2180                qcloud(ind1)=zqsatenv(ind1)
    21282181            ELSE               
    21292182                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))
    21312184                qincloud(ind1)=MAX(fraca(ind1,ind2)*(qcth(ind1,ind2)/cth(ind1,ind2)) &
    21322185                            +(1.-1.*fraca(ind1,ind2))*(qcenv(ind1,ind2)/cenv(ind1,ind2)),0.)
    … …  
    21392192
    21402193     
    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
    21482197                Lv=RLVTT
    21492198        ELSE
    … …  
    21532202
    21542203        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)
    21562205        aenv=1./(1.+(alenv*Lv/cppd))
    2157         senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
     2206        senv=aenv*(po(ind1)-zqsatenvonly(ind1))
    21582207        sth=0.
    21592208     
    2160         sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
     2209        sigma1s=ratqs(ind1,ind2)*zqenvonly(ind1)
    21612210        sigma2s=0.
    21622211
    … …  
    21692218        IF (ctot(ind1,ind2).LT.1.e-3) THEN
    21702219          ctot(ind1,ind2)=0.
    2171           qcloud(ind1)=zqsatenv(ind1,ind2)
     2220          qcloud(ind1)=zqsatenvonly(ind1)
    21722221          qincloud(ind1)=0.
    21732222        ELSE
    2174           qcloud(ind1)=qctot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
     2223          qcloud(ind1)=qctot(ind1,ind2)/ctot(ind1,ind2)+zqsatenvonly(ind1)
    21752224          qincloud(ind1)=MAX(qctot(ind1,ind2)/ctot(ind1,ind2),0.)
    21762225        ENDIF
    … …  
    21942243
    21952244!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    2196 SUBROUTINE ICE_MPC_BL_CLOUDS(ind1,ind2,klev,flag_topthermals,temp,pres,qth,qlth,qith,deltazlev,snowf,fraca,qi)
     2245SUBROUTINE ICE_MPC_BL_CLOUDS(ind1,ind2,klev,Ni,Ei,C_cap,d_top,iflag_topthermals,temp,pres,qth,qsith,qlth,deltazlev,snowf,fraca,qith)
    21972246
    21982247! parameterization of ice for boundary
    … …  
    22602309    IMPLICIT none
    22612310
    2262 
    22632311    INCLUDE "YOMCST.h"
    22642312    INCLUDE "nuage.h"
    22652313
    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
    22682320    REAL,  DIMENSION(klev), INTENT(IN) :: temp       ! temperature [K] within thermals
    22692321    REAL,  DIMENSION(klev), INTENT(IN) :: pres       ! pressure [Pa]
    22702322    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]
    22712324    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]
    22732325    REAL,  DIMENSION(klev), INTENT(IN) :: deltazlev  ! layer thickness [m]
    22742326    REAL,  DIMENSION(klev+1), INTENT(IN) :: snowf      ! snow flux at the upper inferface
    22752327    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]
    22842329
    22852330
    … …  
    22872332    REAL rho(klev)
    22882333    REAL unsurtaudet, unsurtaustardep, unsurtaurim
    2289     REAL qsl, qsi, dqs, AA, BB, Ka, Dv, rhoi
    2290     REAL  p0, t0
     2334    REAL qi, AA, BB, Ka, Dv, rhoi
     2335    REAL p0, t0, fp1, fp2
    22912336    REAL alpha, flux_term
    22922337    REAL det_term, precip_term, rim_term, dep_term
    22932338
    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 
    23162339
    23172340    ind2p1=ind2+1
    23182341    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     ! ==================
    23262342
    23272343    rho=pres/temp/RD  ! air density kg/m3
    … …  
    23342350
    23352351
    2336     IF (flag_topthermals) THEN ! uppermost thermals level, solve a third order polynomial with Cardan's method
     2352    IF (iflag_topthermals .GT. 0) THEN ! uppermost thermals levels
    23372353
    23382354    Dv=0.0001*0.211*(p0/pres(ind2))*((temp(ind2)/T0)**1.94) ! water vapor diffusivity in air, SI
    23392355
    23402356    ! Detrainment term:
     2357
    23412358    unsurtaudet=detr_therm(ind1,ind2)/rho(ind2)/deltazlev(ind2)
    23422359 
    2343     ! vertical flux
    2344    
    2345     flux_term=d_top*fm_therm(ind1,ind2)/deltazlev(ind2)
    23462360
    23472361    ! 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)
    23502362    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)
    23532366
    23542367    ! Riming term  neglected at this level
    23552368    !unsurtaurim=rho(ind2)*alpha*3./rhoi/2.*Ei*qlth(ind2)*((p0/pres(ind2))**0.4)
    23562369
    2357     qi=rho(ind2)*unsurtaustardep/MAX((rho(ind2)*unsurtaudet-flux_term),1E-12)
     2370    qi=fraca(ind2)*unsurtaustardep/MAX((d_top*unsurtaudet),1E-12)
    23582371    qi=MAX(qi,0.)**(3./2.)
    23592372
    … …  
    23702383    ! Deposition term
    23712384
    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)
    23742385    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
     2386    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
    23782389 
    23792390    ! Riming term
    23802391
    23812392    unsurtaurim=rho(ind2p1)*alpha*3./rhoi/2.*Ei*qlth(ind2p1)*((p0/pres(ind2p1))**0.4)
    2382     rim_term=rho(ind2p1)*qith(ind2p1)*unsurtaurim
     2393    rim_term=rho(ind2p1)*fraca(ind2p1)*qith(ind2p1)*unsurtaurim
    23832394
    23842395    ! Precip term
    23852396
    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)
    23892409
    23902410    ! Calculation in a top-to-bottom loop
    23912411
    23922412    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)
    24032421
    24042422    RETURN

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

    r4062 r4072  
    175175    INTEGER,SAVE :: iflag_clw_omp
    176176    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
    178178    LOGICAL,SAVE :: reevap_ice_omp
    179179    INTEGER,SAVE :: iflag_pdf_omp
    … …  
    11311131    CALL getin('coef_eva',coef_eva_omp)
    11321132    !
     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    !
    11331141    !Config Key  = reevap_ice
    11341142    !Config Desc = 
    … …  
    24922500    ffallv_con = ffallv_con_omp
    24932501    coef_eva = coef_eva_omp
     2502    coef_eva_i = coef_eva_i_omp
    24942503    reevap_ice = reevap_ice_omp
    24952504    iflag_pdf = iflag_pdf_omp
    … …  
    29142923    WRITE(lunout,*) ' ffallv_con = ', ffallv_con
    29152924    WRITE(lunout,*) ' coef_eva = ', coef_eva
     2925    WRITE(lunout,*) ' coef_eva_i = ', coef_eva_i
    29162926    WRITE(lunout,*) ' reevap_ice = ', reevap_ice
    29172927    WRITE(lunout,*) ' iflag_pdf = ', iflag_pdf

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

    r2415 r4072  
    1111      REAL cld_tau_lsc,cld_tau_con
    1212      REAL ffallv_lsc,ffallv_con
    13       REAL coef_eva
     13      REAL coef_eva,coef_eva_i
    1414      LOGICAL reevap_ice
    1515      INTEGER iflag_pdf
    … …  
    2525     &     ,ffallv_con                                                  &
    2626     &     ,coef_eva                                                    &
     27     &     ,coef_eva_i                                                  &
    2728     &     ,reevap_ice                                                  &
    2829     &     ,iflag_fisrtilp_qsat                                         &

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

    r4062 r4072  
    66
    77!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    8 SUBROUTINE LSCP(dtime,missing_val,                      & 
     8SUBROUTINE LSCP(dtime,missing_val,                      &
    99     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,             &
    1111     radliq, radicefrac, rain, snow,                    &
    1212     pfrac_impa, pfrac_nucl, pfrac_1nucl,               &
    … …  
    2525!     
    2626! This code is a new version of the fisrtilp.F90 routine, which is itself a
    27 ! fusion of 'first' (superrsaturation physics, P. LeVan K. Laval)
     27! merge of 'first' (superrsaturation physics, P. LeVan K. Laval)
    2828! and 'ilp' (il pleut, L. Li)
    2929!
    … …  
    126126  INTEGER,                         INTENT(IN)   :: iflag_ice_thermo! flag to activate the ice thermodynamics
    127127                                                                   ! 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
    129129
    130130  LOGICAL, DIMENSION(klon,klev),   INTENT(IN)   :: ptconv          ! grid points where deep convection scheme is active
    … …  
    138138  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: zpspsk          ! exner potential (p/100000)**(R/cp)
    139139  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: ztla            ! liquid temperature within thermals [K]
    140   REAL, DIMENSION(klon,klev),      INTENT(INOUT)   :: zthl            ! liquid potential temperature [K]
     140  REAL, DIMENSION(klon,klev),      INTENT(INOUT)   :: zthl         ! liquid potential temperature [K]
    141141
    142142  ! INPUT/OUTPUT variables
    … …  
    145145  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: ratqs            ! function of pressure that sets the large-scale
    146146                                                                    ! cloud PDF (sigma=ratqs*qt)
     147 
    147148
    148149  ! Input sursaturation en glace
    … …  
    204205  PARAMETER (ztfondue=278.15)
    205206
    206   REAL, SAVE    :: rain_int_min=0.001                               ! Minimum local rain intensity [mm/s] before the decrease in  associates precipitation fraction
     207  REAL, SAVE    :: rain_int_min=0.001                               ! Minimum local rain intensity [mm/s] before the decrease in associated precipitation fraction
    207208  !$OMP THREADPRIVATE(rain_int_min)
    208209
    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
    211216  ! LOCAL VARIABLES:
    212217  !----------------
    213218
    214219
    215   REAL qsl, qsi
     220  REAL qsl(klon), qsi(klon)
    216221  REAL zct, zcl
    217222  REAL ctot(klon,klev)
    218223  REAL ctot_vol(klon,klev)
    219   INTEGER mpc_bl_points(klon,klev)
    220224  REAL zqs(klon), zdqs(klon), zdelta, zcor, zcvm5 
    221225  REAL zdqsdT_raw(klon)
    … …  
    230234  REAL zifl(klon), zifln(klon), zqev0,zqevi, zqevti
    231235  REAL zoliq(klon), zcond(klon), zq(klon), zqn(klon)
    232   REAL zoliqp(klon), zoliqi(klon)
     236  REAL zoliql(klon), zoliqi(klon)
    233237  REAL zt(klon)
    234   REAL zdz(klon),zrho(klon),ztot, zrhol(klon)
     238  REAL zdz(klon),zrho(klon),iwc(klon)
    235239  REAL zchau,zfroi,zfice(klon),zneb(klon),znebprecip(klon)
    236   REAL zmelt,zpluie,zice
     240  REAL zmelt,zrain,zsnow,zprecip
    237241  REAL dzfice(klon)
    238   REAL zsolid,qtot,dqsl,dqsi
     242  REAL zsolid
     243  REAL qtot(klon), qzero(klon)
     244  REAL dqsl(klon),dqsi(klon)
    239245  REAL smallestreal
    240246  !  Variables for Bergeron process
    … …  
    262268  REAL tot_zneb(klon), tot_znebn(klon), d_tot_zneb(klon)                                             
    263269  REAL d_znebprecip_clr_cld(klon), d_znebprecip_cld_clr(klon)
    264   REAL velo(klon)
     270  REAL velo(klon,klev), vr
    265271  REAL qlmpc, qimpc, rnebmpc
    266   REAL radliqi(klon,klev)
     272  REAL radliqi(klon,klev), radliql(klon,klev)
    267273
    268274  INTEGER i, k, n, kk
    269275  INTEGER n_i(klon), iter
    270276  INTEGER ncoreczq 
     277  INTEGER mpc_bl_points(klon,klev)
    271278  INTEGER,SAVE :: itap=0
    272279  !$OMP THREADPRIVATE(itap)
    … …  
    334341    CALL getin_p('seuil_neb',seuil_neb)
    335342    CALL getin_p('rain_int_min',rain_int_min)
     343    CALL getin_p('ok_radliq_precip',ok_radliq_precip)
    336344    WRITE(lunout,*) 'lscp, ninter:', ninter
    337345    WRITE(lunout,*) 'lscp, iflag_evap_prec:', iflag_evap_prec
    338346    WRITE(lunout,*) 'lscp, seuil_neb:', seuil_neb
    339347    WRITE(lunout,*) 'lscp, rain_int_min:', rain_int_min
     348    WRITE(lunout,*) 'lscp, ok_radliq_precip:', ok_radliq_precip
    340349
    341350    ! check for precipitation sub-time steps
    … …  
    397406tot_znebn(:) = 0.0                                                                               
    398407d_tot_zneb(:) = 0.0   
     408qzero(:)=0.0
    399409
    400410!--ice sursaturation
    … …  
    468478    ! --------------------------------------------------------------------
    469479    ! A part of the precipitation coming from above is evaporated/sublimated.
    470     !
    471480    ! --------------------------------------------------------------------
    472481
    … …  
    474483    IF (iflag_evap_prec.GE.1) THEN
    475484
     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(:))
    476492
    477493        DO i = 1, klon
    … …  
    480496            IF (zrfl(i)+zifl(i).GT.0.) THEN
    481497                   
    482                 CALL CALC_QSAT_ECMWF(zt(i),0.,pplay(i,k),RTT,0,.false.,zqs(i),zdqs(i))
    483 
    484498            ! LTP: we only account for precipitation evaporation in the clear-sky (iflag_evap_prec=4).
    485499                IF (iflag_evap_prec.EQ.4) THEN                                                         
    … …  
    503517                ENDIF   
    504518
    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
    511520                ! Evaporation of liquid precipitation coming from above
    512521                ! dP/dz=beta*(1-q/qsat)*sqrt(P)
    … …  
    515524
    516525                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)) &
    518527                    *SQRT(zrfl(i)/max(1.e-4,znebprecip(i)))        &
    519528                    *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
    520529                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)) &
    522531                    *SQRT(zrfl(i)/max(1.e-8,znebprecipclr(i))) &
    523532                    *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
    524533                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)) &
    526535                    *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
    527536                ENDIF
    … …  
    531540                *RG*dtime/(paprs(i,k)-paprs(i,k+1))
    532541         
    533                 ! qsat wrt ice according to thermcep
    534                 CALL CALC_QSAT_ECMWF(zt(i),0.,pplay(i,k),RTT,2,.false.,qsi,dqsi)
    535 
    536542                ! sublimation of the solid precipitation coming from above
    537543                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)) &
    539545                    *SQRT(zifl(i)/max(1.e-4,znebprecip(i))) &
    540546                    *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
    541547                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)) &
    543549                    *SQRT(zifl(i)/max(1.e-8,znebprecipclr(i))) &
    544550                    *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
    545551                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)) &
    547553                    *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
    548554                ENDIF
    … …  
    566572
    567573
    568                 ! EV: agrees with JL's comments below so correct and comment
    569                 ! JLD: I do not understand the lines below. We distribute the liquid
    570                 ! and solid parts of the precipitation even if the layer does not get
    571                 ! saturated. To my opinion, we should all replace with:
    572                 ! zqev=zqevt
    573                 ! zqevi=zqevti
    574                 !    IF (zqevt+zqevti.GT.0.) THEN
    575                 !        zqev=MIN(zqev0*zqevt/(zqevt+zqevti),zqevt)
    576                 !        zqevi=MIN(zqev0*zqevti/(zqevt+zqevti),zqevti)
    577                 !    ELSE
    578                 !        zqev=0.
    579                 !        zqevi=0.
    580                 !    ENDIF
    581                 ! ENDIF
    582 
    583574                ! New solid and liquid precipitation fluxes
    584575                zrfln(i) = Max(0.,zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1))  &
    … …  
    659650    ! Calculation of qsat, L/Cp*dqsat/dT and ncoreczq counter
    660651    !-------------------------------------------------------
    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
    663656            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))
    666657            zdqsdT_raw(i) = zdqs(i)*RCPD*(1.0+RVTMP2*zq(i)) / (RLVTT*(1.-zdelta) + RLSTT*zdelta)
    667          
    668658            IF (zq(i) .LT. 1.e-15) THEN
    669659                ncoreczq=ncoreczq+1
    … …  
    671661            ENDIF
    672662
    673         ENDDO
     663     ENDDO
    674664   
    675665
    … …  
    777767                DO i=1,klon
    778768
    779                     ! convergence = .true. since convergence is not satisfied
     769                    ! convergence = .true. until when convergence is not satisfied
    780770                    convergence(i)=ABS(DT(i)).GT.DDT0
    781771                   
    … …  
    796786
    797787                        ! 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
    804806
    805807                        zpdf_sig(i)=ratqs(i,k)*zq(i)
    … …  
    826828                          ENDIF
    827829
    828                           rnebss(i,k)=0.0   !--ajout 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)
    829831                          fcontrN(i,k)=0.0  !--idem
    830832                          fcontrP(i,k)=0.0  !--idem
    831833                          qss(i,k)=0.0      !--idem
    832 
     834                       
    833835                        ELSE
     836               
    834837                        !------------------------------------
    835                         ! SURSATURATION EN GLACE
     838                        ! ICE SUPERSATURATION
    836839                        !------------------------------------
    837840
    838841                        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),                                 & 
    841844                             Tcontr(i,k), qcontr(i,k), qcontr2(i,k), fcontrN(i,k), fcontrP(i,k) )
    842845
    … …  
    850853
    851854
    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
    860861               
    861862                        IF (zfice(i).LT.1) THEN
    … …  
    896897
    897898        ! 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(:))
    899900
    900901        DO i=1, klon
    … …  
    10021003
    10031004    DO i = 1, klon
    1004         IF (rneb(i,k).GT.0.0) THEN
    10051005            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)
    10101012            radliq(i,k) = zoliq(i)/REAL(ninter+1)
    10111013            radicefrac(i,k) = zfice(i)
    10121014            radliqi(i,k)=zoliq(i)*zfice(i)/REAL(ninter+1)
    1013         ENDIF
     1015            radliql(i,k)=zoliq(i)*(1.-zfice(i))/REAL(ninter+1)
    10141016    ENDDO
    10151017
    … …  
    10241026        DO i=1,klon
    10251027            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
    10271029            ENDIF
    10281030        ENDDO
    10291031
    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))
    10311033
    10321034        DO i = 1, klon
    … …  
    10341036            IF (rneb(i,k).GT.0.0) THEN
    10351037
    1036                 ! Initialization of zpluie, zice and ztot:
    1037                 zpluie=0.
    1038                 zice=0.
    1039                 ztot=0.
     1038                ! Initialization of zrain, zsnow and zprecip:
     1039                zrain=0.
     1040                zsnow=0.
     1041                zprecip=0.
    10401042
    10411043                IF (zneb(i).EQ.seuil_neb) THEN
    1042                     ztot = 0.0
    1043                     zice = 0.0
    1044                     zpluie= 0.0
     1044                    zprecip = 0.0
     1045                    zsnow = 0.0
     1046                    zrain= 0.0
    10451047                ELSE
    10461048                    ! water quantity to remove: zchau (Sundqvist, 1978)
    … …  
    10491051                        zcl=cld_lc_con
    10501052                        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)
    10521054                    ELSE
    10531055                        zcl=cld_lc_lsc
    10541056                        zct=1./cld_tau_lsc
    10551057                        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)
    10571059                    ENDIF
    10581060
    … …  
    10611063                    ! surface fraction (which is larger and artificially
    10621064                    ! reduces the in-cloud water).
     1065
    10631066                    IF ((iflag_cloudth_vert.GE.3).AND.(iflag_rain_incloud_vol.EQ.1)) THEN
    10641067                        zchau = zct   *dtime/REAL(ninter) * zoliq(i) &
    … …  
    10691072                    ENDIF
    10701073
    1071                     zpluie = MIN(MAX(zchau,0.0),zoliq(i)*(1.-zfice(i)))
    1072                     zice   = MIN(MAX(zfroi,0.0),zoliq(i)*zfice(i))
    1073                     ztot   = 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)
    10741077
    10751078                ENDIF
    10761079           
    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)
    10841083
    10851084                radliq(i,k) = radliq(i,k) + zoliq(i)/REAL(ninter+1)
     1085                radliql(i,k) = radliql(i,k) + zoliql(i)/REAL(ninter+1)
    10861086                radliqi(i,k) = radliqi(i,k) + zoliqi(i)/REAL(ninter+1)
    10871087           
    1088             ENDIF
     1088            ENDIF ! rneb >0
    10891089
    10901090        ENDDO  ! i = 1,klon
    … …  
    10931093
    10941094   
    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
    10961106    DO i=1,klon
    1097         IF (radliq(i,k) .GT. 0) THEN
     1107        IF (radliq(i,k) .GT. 0.) THEN
    10981108            radicefrac(i,k)=MIN(MAX(radliqi(i,k)/radliq(i,k),0.),1.)
    10991109        ENDIF
    … …  
    11011111
    11021112
    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
    11141114
    11151115    DO i = 1, klon
    11161116
    11171117            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
    11401152           
    1141                 d_ql(i,k) = (1-zfice(i))*zoliq(i)
    1142                 d_qi(i,k) = zfice(i)*zoliq(i)
    1143        
    11441153                IF (iflag_evap_prec.EQ.4) THEN                                                           
    11451154                    zrflcld(i) = zrflcld(i)+zqprecl(i) &                                                   
    … …  
    11651174    IF (iflag_evap_prec.EQ.4) THEN 
    11661175
    1167         DO i=1, klon                                       
     1176        DO i=1,klon                                       
    11681177               
    11691178            IF ((zrflclr(i) + ziflclr(i)) .GT. 0. ) THEN     
    … …  
    11921201    ! Outputs:
    11931202    ! Precipitation fluxes at layer interfaces
    1194     ! and temperature and vapor tendencies
     1203    ! and temperature and water species tendencies
    11951204    DO i = 1, klon
    11961205        psfl(i,k)=zifl(i)
    11971206        prfl(i,k)=zrfl(i)
     1207        d_ql(i,k) = (1-zfice(i))*zoliq(i)
     1208        d_qi(i,k) = zfice(i)*zoliq(i)
    11981209        d_q(i,k) = zq(i) - q(i,k)
    11991210        d_t(i,k) = zt(i) - t(i,k)
    … …  
    12641275    ENDDO
    12651276     
    1266     !--save some variables for ice sursaturation
     1277    !--save some variables for ice supersaturation
    12671278    !
    12681279    DO i = 1, klon
    1269         ! pour la mémoire
     1280        ! for memory
    12701281        rneb_seri(i,k) = rneb(i,k)
    12711282
    1272         ! pour les diagnostics
     1283        ! for diagnostics
    12731284        rnebclr(i,k) = 1.0 - rneb(i,k) - rnebss(i,k)
    12741285
    12751286        qvc(i,k) = zqs(i) * rneb(i,k)
    1276         qclr(i,k) = MAX(1.e-10,zq(i) - qvc(i,k) - qss(i,k))  !--ajout OB a cause de cas pathologiques avec lognormale=F
     1287        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
    12771288        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
     1294END DO
    12861295
    12871296!======================================================================

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

    r4059 r4072  
    4747
    4848        tempc=temp(i)-273.15 ! celcius temp
    49         iwcg=iwc(i)*1000.    ! iwc in g/m3
     49        iwcg=MAX(iwc(i)*1000.,1E-3) ! iwc in g/m3. We set a minimum value to prevent from division by 0
    5050        phpa=pres(i)/100.    ! pressure in hPa
    5151
    … …  
    7676                **(1./(0.807+bice*0.00581+0.0457*bice**2))
    7777
    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)
    8181
    8282       
    … …  
    9393
    9494        vmsnow=vmsnow*((1000./phpa)**0.35)
    95 
    9695        velo(i)=fallv_tun*min(vmsnow,vmice)/100. ! to m/s
    9796        dvel=0.2
    98         cvel=velo(i)/((iwc(i)*rho(i))**dvel)
     97        cvel=velo(i)/((iwcg/1000.*rho(i))**dvel)
    9998
    10099    ELSE
    101100        ! 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)
    103102        dvel=0.16
    104103        cvel=fallv_tun*3.29/2.0*(rho(i)**0.16)
    105104    ENDIF
    106 
    107105    ENDDO
    108106
    … …  
    202200    ENDDO
    203201
     202
    204203    RETURN
    205204
    … …  
    209208
    210209
    211 SUBROUTINE CALC_QSAT_ECMWF(temp,qtot,pressure,tref,phase,flagth,qs,dqs)
     210SUBROUTINE CALC_QSAT_ECMWF(klon,temp,qtot,pressure,tref,phase,flagth,qs,dqs)
    212211!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    213212    ! Calculate qsat following ECMWF method
    214213!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     214
    215215
    216216    IMPLICIT NONE
    … …  
    220220    include "FCTTRE.h"
    221221
    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
    226227    LOGICAL, INTENT(IN) :: flagth     ! flag for qsat calculation for thermals
    227 
    228228    INTEGER, INTENT(IN) :: phase
    229229    ! phase: 0=depend on temperature sign (temp>tref -> liquid, temp<tref, solid)
    … …  
    231231    !        2=solid
    232232
    233     REAL, INTENT(OUT) :: qs      ! saturation specific humidity [kg/kg]
    234     REAL, INTENT(OUT) :: dqs     ! derivation of saturation specific humidity wrt T
     233    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
    235235
    236236    REAL delta, cor, cvm5
    237    
     237    INTEGER i
     238
     239    DO i=1,klon
     240
    238241    IF (phase .EQ. 1) THEN
    239242        delta=0.
    … …  
    241244        delta=1.
    242245    ELSE
    243         delta=MAX(0.,SIGN(1.,tref-temp))
     246        delta=MAX(0.,SIGN(1.,tref-temp(i)))
    244247    ENDIF
    245248
    … …  
    248251    ELSE
    249252    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
    258263
    259264END SUBROUTINE CALC_QSAT_ECMWF
    … …  
    262267
    263268!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    264 SUBROUTINE CALC_GAMMASAT(temp,qtot,pressure,gammasat,dgammasatdt)
     269SUBROUTINE CALC_GAMMASAT(klon,temp,qtot,pressure,gammasat,dgammasatdt)
    265270
    266271!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    … …  
    279284    include "nuage.h"
    280285
    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
    291297    REAL, PARAMETER :: acirrus=2.349
    292298    REAL, PARAMETER :: bcirrus=259.0
    293299
    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
    299308            ! 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)) THEN
     309            gammasat(i)=1.
     310            dgammasatdt(i)=0.
     311
     312        ELSEIF ((temp(i) .LT. RTT) .AND. (temp(i) .GT. t_glace_min)) THEN
    304313           
    305314            IF (iflag_gammasat .GE. 2) THEN         
    306                gammasat=qsl/qsi
    307                dgammasatdt=(dqsl*qsi-dqsi*qsl)/qsi/qsi
     315               gammasat(i)=qsl(i)/qsi(i)
     316               dgammasatdt(i)=(dqsl(i)*qsi(i)-dqsi(i)*qsl(i))/qsi(i)/qsi(i)
    308317            ELSE
    309                gammasat=1.
    310                dgammasatdt=0.
     318               gammasat(i)=1.
     319               dgammasatdt(i)=0.
    311320            ENDIF
    312321
    … …  
    317326            ! Koop, 2000 and Karcher 2008, QJRMS
    318327            ! 'Cirrus regime'
    319                fcirrus=acirrus-temp/bcirrus
    320                IF (fcirrus .LT. qsl/qsi) THEN
    321                   gammasat=qsl/qsi
    322                   dgammasatdt=(dqsl*qsi-dqsi*qsl)/qsi/qsi
     328               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)
    323332               ELSE
    324                   gammasat=fcirrus
    325                   dgammasatdt=-1.0/bcirrus
     333                  gammasat(i)=fcirrus
     334                  dgammasatdt(i)=-1.0/bcirrus
    326335               ENDIF
    327336           
    328337            ELSE
    329338
    330                gammasat=1.
    331                dgammasatdt=0.
     339               gammasat(i)=1.
     340               dgammasatdt(i)=0.
    332341
    333342            ENDIF
    … …  
    335344        ENDIF
    336345   
     346    END DO
     347
     348
    337349END SUBROUTINE CALC_GAMMASAT
    338350

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

    r4059 r4072  
    438438      REAL,ALLOCATABLE,SAVE,DIMENSION(:,:) :: cldemi, cldfra, cldtau, fiwc, fl, re, flwc
    439439!$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)
    442442      REAL,ALLOCATABLE,SAVE,DIMENSION(:,:) :: ref_liq, ref_ice, theta, zphi
    443443!$OMP THREADPRIVATE(ref_liq, ref_ice, theta, zphi)
    … …  
    458458      REAL,ALLOCATABLE,SAVE,DIMENSION(:,:,:) :: phiwriteSTD, qwriteSTD, twriteSTD, rhwriteSTD
    459459!$OMP THREADPRIVATE(phiwriteSTD, qwriteSTD, twriteSTD, rhwriteSTD)
     460
    460461
    461462! ug et d'autres encore:
    … …  
    777778      ALLOCATE(rain_lsc(klon))
    778779      ALLOCATE(rain_num(klon))
    779       ALLOCATE(qlth(klon,klev), qith(klon,klev))
     780      ALLOCATE(qlth(klon,klev), qith(klon,klev), qsith(klon,klev))
    780781      !
    781782      ALLOCATE(sens_x(klon), sens_w(klon))
    … …  
    823824      ALLOCATE(ref_liq_pi(klon, klev), ref_ice_pi(klon, klev))
    824825      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
    825827      ALLOCATE(pmfd(klon, klev), pmfu(klon, klev))
    826828
    … …  
    10871089      DEALLOCATE(rain_lsc)
    10881090      DEALLOCATE(rain_num)
    1089       DEALLOCATE(qlth, qith)
     1091      DEALLOCATE(qlth, qith, qsith)
    10901092!
    10911093      DEALLOCATE(sens_x, sens_w)

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

    r4065 r4072  
    559559      REAL,ALLOCATABLE,SAVE,DIMENSION(:,:) :: cldemi, cldfra, cldtau, fiwc, fl, re, flwc
    560560!$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)
    563563      REAL,ALLOCATABLE,SAVE,DIMENSION(:,:) :: ref_liq, ref_ice, theta, zphi
    564564!$OMP THREADPRIVATE(ref_liq, ref_ice, theta, zphi)
    … …  
    943943      ALLOCATE(rain_lsc(klon))
    944944      ALLOCATE(rain_num(klon))
    945       ALLOCATE(qlth(klon,klev), qith(klon,klev))
     945      ALLOCATE(qlth(klon,klev), qith(klon,klev), qsith(klon,klev))
    946946      !
    947947#ifdef ISO
    … …  
    10041004      ALLOCATE(ref_liq_pi(klon, klev), ref_ice_pi(klon, klev))
    10051005      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
    10061007      ALLOCATE(pmfd(klon, klev), pmfu(klon, klev))
    10071008
    … …  
    13241325      DEALLOCATE(rain_lsc)
    13251326      DEALLOCATE(rain_num)
    1326       DEALLOCATE(qlth, qith)
     1327      DEALLOCATE(qlth, qith, qsith)
    13271328!
    13281329      DEALLOCATE(sens_x, sens_w)
Note: See TracChangeset for help on using the changeset viewer.