Index: /LMDZ6/trunk/libf/phylmd/lmdz_lscp_tools.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/lmdz_lscp_tools.F90	(revision 4665)
+++ /LMDZ6/trunk/libf/phylmd/lmdz_lscp_tools.F90	(revision 4665)
@@ -0,0 +1,456 @@
+MODULE lmdz_lscp_tools
+
+    IMPLICIT NONE
+
+CONTAINS
+
+!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+SUBROUTINE FALLICE_VELOCITY(klon,iwc,temp,rho,pres,ptconv,velo)
+
+    ! Ref:
+    ! Stubenrauch, C. J., Bonazzola, M.,
+    ! Protopapadaki, S. E., & Musat, I. (2019).
+    ! New cloud system metrics to assess bulk
+    ! ice cloud schemes in a GCM. Journal of
+    ! Advances in Modeling Earth Systems, 11,
+    ! 3212–3234. https://doi.org/10.1029/2019MS001642
+    
+    use lmdz_lscp_ini, only: iflag_vice, ffallv_con, ffallv_lsc
+    use lmdz_lscp_ini, only: cice_velo, dice_velo 
+
+    IMPLICIT NONE
+
+    INTEGER, INTENT(IN) :: klon
+    REAL, INTENT(IN), DIMENSION(klon) :: iwc       ! specific ice water content [kg/m3]
+    REAL, INTENT(IN), DIMENSION(klon) :: temp      ! temperature [K]
+    REAL, INTENT(IN), DIMENSION(klon) :: rho       ! dry air density [kg/m3]
+    REAL, INTENT(IN), DIMENSION(klon) :: pres      ! air pressure [Pa]
+    LOGICAL, INTENT(IN), DIMENSION(klon) :: ptconv    ! convective point  [-]
+
+    REAL, INTENT(OUT), DIMENSION(klon) :: velo    ! fallspeed velocity of crystals [m/s]
+
+
+    INTEGER i
+    REAL logvm,iwcg,tempc,phpa,fallv_tun
+    REAL m2ice, m2snow, vmice, vmsnow
+    REAL aice, bice, asnow, bsnow
+    
+
+    DO i=1,klon
+
+        IF (ptconv(i)) THEN
+            fallv_tun=ffallv_con
+        ELSE
+            fallv_tun=ffallv_lsc
+        ENDIF
+
+        tempc=temp(i)-273.15 ! celcius temp
+        iwcg=MAX(iwc(i)*1000.,1E-3) ! iwc in g/m3. We set a minimum value to prevent from division by 0
+        phpa=pres(i)/100.    ! pressure in hPa
+
+    IF (iflag_vice .EQ. 1) THEN
+        ! so-called 'empirical parameterization' in Stubenrauch et al. 2019
+        if (tempc .GE. -60.0) then
+
+            logvm= -0.0000414122*tempc*tempc*log(iwcg)-0.00538922*tempc*log(iwcg) &
+                    -0.0516344*log(iwcg)+0.00216078*tempc + 1.9714    
+            velo(i)=exp(logvm)
+        else
+            velo(i)=65.0*(iwcg**0.2)*(150./phpa)**0.15
+        endif
+        
+        velo(i)=fallv_tun*velo(i)/100.0 ! from cm/s to m/s
+
+    ELSE IF (iflag_vice .EQ. 2) THEN
+        ! so called  PSDM empirical coherent bulk ice scheme in Stubenrauch et al. 2019
+        aice=0.587
+        bice=2.45
+        asnow=0.0444
+        bsnow=2.1
+        
+        m2ice=((iwcg*0.001/aice)/(exp(13.6-bice*7.76+0.479*bice**2)* & 
+                exp((-0.0361+bice*0.0151+0.00149*bice**2)*tempc)))   &
+                **(1./(0.807+bice*0.00581+0.0457*bice**2))
+
+        vmice=100.*1042.4*exp(13.6-(bice+1)*7.76+0.479*(bice+1.)**2)*exp((-0.0361+ &
+                 (bice+1.)*0.0151+0.00149*(bice+1.)**2)*tempc) &
+                 *(m2ice**(0.807+(bice+1.)*0.00581+0.0457*(bice+1.)**2))/(iwcg*0.001/aice)
+
+       
+        vmice=vmice*((1000./phpa)**0.2)
+      
+        m2snow=((iwcg*0.001/asnow)/(exp(13.6-bsnow*7.76+0.479*bsnow**2)* &
+               exp((-0.0361+bsnow*0.0151+0.00149*bsnow**2)*tempc)))         &
+               **(1./(0.807+bsnow*0.00581+0.0457*bsnow**2))
+
+
+        vmsnow=100.*14.3*exp(13.6-(bsnow+.416)*7.76+0.479*(bsnow+.416)**2)&
+                  *exp((-0.0361+(bsnow+.416)*0.0151+0.00149*(bsnow+.416)**2)*tempc)&
+                  *(m2snow**(0.807+(bsnow+.416)*0.00581+0.0457*(bsnow+.416)**2))/(iwcg*0.001/asnow)
+
+        vmsnow=vmsnow*((1000./phpa)**0.35)
+        velo(i)=fallv_tun*min(vmsnow,vmice)/100. ! to m/s
+
+    ELSE
+        ! By default, fallspeed velocity of ice crystals according to Heymsfield & Donner 1990
+        velo(i) = fallv_tun*cice_velo*((iwcg/1000.)**dice_velo)
+    ENDIF
+    ENDDO
+
+END SUBROUTINE FALLICE_VELOCITY
+!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+SUBROUTINE ICEFRAC_LSCP(klon, temp, iflag_ice_thermo, distcltop, temp_cltop, icefrac, dicefracdT)
+!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+  
+  ! Compute the ice fraction 1-xliq (see e.g.
+  ! Doutriaux-Boucher & Quaas 2004, section 2.2.)
+  ! as a function of temperature
+  ! see also Fig 3 of Madeleine et al. 2020, JAMES
+!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+
+    USE print_control_mod, ONLY: lunout, prt_level
+    USE lmdz_lscp_ini, ONLY: t_glace_min, t_glace_max, exposant_glace, iflag_t_glace
+    USE lmdz_lscp_ini, ONLY : RTT, dist_liq, temp_nowater
+
+    IMPLICIT NONE
+
+
+    INTEGER, INTENT(IN)                 :: klon              ! number of horizontal grid points
+    REAL, INTENT(IN), DIMENSION(klon)   :: temp              ! temperature
+    REAL, INTENT(IN), DIMENSION(klon)   :: distcltop         ! distance to cloud top
+    REAL, INTENT(IN), DIMENSION(klon)   :: temp_cltop        ! temperature of cloud top
+    INTEGER, INTENT(IN)                 :: iflag_ice_thermo
+    REAL, INTENT(OUT), DIMENSION(klon)  :: icefrac
+    REAL, INTENT(OUT), DIMENSION(klon)  :: dicefracdT
+
+
+    INTEGER i
+    REAL    liqfrac_tmp, dicefrac_tmp
+    REAL    Dv, denomdep,beta,qsi,dqsidt
+    LOGICAL ice_thermo
+
+    CHARACTER (len = 20) :: modname = 'lscp_tools'
+    CHARACTER (len = 80) :: abort_message
+
+    IF ((iflag_t_glace.LT.2) .OR. (iflag_t_glace.GT.6)) THEN
+       abort_message = 'lscp cannot be used if iflag_t_glace<2 or >6'
+       CALL abort_physic(modname,abort_message,1)
+    ENDIF
+
+    IF (.NOT.((iflag_ice_thermo .EQ. 1).OR.(iflag_ice_thermo .GE. 3))) THEN
+       abort_message = 'lscp cannot be used without ice thermodynamics'
+       CALL abort_physic(modname,abort_message,1)
+    ENDIF
+
+
+    DO i=1,klon
+ 
+        ! old function with sole dependence upon temperature
+        IF (iflag_t_glace .EQ. 2) THEN
+            liqfrac_tmp = (temp(i)-t_glace_min) / (t_glace_max-t_glace_min)
+            liqfrac_tmp = MIN(MAX(liqfrac_tmp,0.0),1.0)
+            icefrac(i) = (1.0-liqfrac_tmp)**exposant_glace
+            IF (icefrac(i) .GT.0.) THEN
+                 dicefracdT(i)= exposant_glace * (icefrac(i)**(exposant_glace-1.)) &
+                           / (t_glace_min - t_glace_max)
+            ENDIF
+
+            IF ((icefrac(i).EQ.0).OR.(icefrac(i).EQ.1)) THEN
+                 dicefracdT(i)=0.
+            ENDIF
+
+        ENDIF
+
+        ! function of temperature used in CMIP6 physics
+        IF (iflag_t_glace .EQ. 3) THEN
+            liqfrac_tmp = (temp(i)-t_glace_min) / (t_glace_max-t_glace_min)
+            liqfrac_tmp = MIN(MAX(liqfrac_tmp,0.0),1.0)
+            icefrac(i) = 1.0-liqfrac_tmp**exposant_glace
+            IF ((icefrac(i) .GT.0.) .AND. (liqfrac_tmp .GT. 0.)) THEN
+                dicefracdT(i)= exposant_glace * ((liqfrac_tmp)**(exposant_glace-1.)) &
+                          / (t_glace_min - t_glace_max)
+            ELSE
+                dicefracdT(i)=0.
+            ENDIF
+        ENDIF
+
+        ! for iflag_t_glace .GE. 4, the liquid fraction depends upon temperature at cloud top
+        ! and then decreases with decreasing height 
+
+        !with linear function of temperature at cloud top
+        IF (iflag_t_glace .EQ. 4) THEN  
+                liqfrac_tmp = (temp(i)-t_glace_min) / (t_glace_max-t_glace_min)
+                liqfrac_tmp = MIN(MAX(liqfrac_tmp,0.0),1.0)
+                icefrac(i)    =  MAX(MIN(1.,1.0 - liqfrac_tmp*exp(-distcltop(i)/dist_liq)),0.)
+                dicefrac_tmp = - temp(i)/(t_glace_max-t_glace_min)
+                dicefracdT(i) = dicefrac_tmp*exp(-distcltop(i)/dist_liq)
+                IF ((liqfrac_tmp .LE.0) .OR. (liqfrac_tmp .GE. 1)) THEN
+                        dicefracdT(i) = 0.
+                ENDIF
+        ENDIF
+
+        ! with CMIP6 function of temperature at cloud top
+        IF (iflag_t_glace .EQ. 5) THEN 
+                liqfrac_tmp = (temp(i)-t_glace_min) / (t_glace_max-t_glace_min)
+                liqfrac_tmp =  MIN(MAX(liqfrac_tmp,0.0),1.0)
+                liqfrac_tmp = liqfrac_tmp**exposant_glace
+                icefrac(i)  =  MAX(MIN(1.,1.0 - liqfrac_tmp*exp(-distcltop(i)/dist_liq)),0.)
+                IF ((liqfrac_tmp .LE.0) .OR. (liqfrac_tmp .GE. 1)) THEN
+                        dicefracdT(i) = 0.
+                ELSE
+                        dicefracdT(i) = exposant_glace*((liqfrac_tmp)**(exposant_glace-1.))/(t_glace_min- t_glace_max) &
+                                        *exp(-distcltop(i)/dist_liq)
+                ENDIF
+        ENDIF
+
+        ! with modified function of temperature at cloud top 
+        ! to get largere values around 260 K, works well with t_glace_min = 241K
+        IF (iflag_t_glace .EQ. 6) THEN 
+                IF (temp(i) .GT. t_glace_max) THEN
+                        liqfrac_tmp = 1.
+                ELSE
+                        liqfrac_tmp = -((temp(i)-t_glace_max) / (t_glace_max-t_glace_min))**2+1.
+                ENDIF
+                liqfrac_tmp  = MIN(MAX(liqfrac_tmp,0.0),1.0)
+                icefrac(i)   = MAX(MIN(1.,1.0 - liqfrac_tmp*exp(-distcltop(i)/dist_liq)),0.)        
+                IF ((liqfrac_tmp .LE.0) .OR. (liqfrac_tmp .GE. 1)) THEN
+                        dicefracdT(i) = 0.
+                ELSE
+                        dicefracdT(i) = 2*((temp(i)-t_glace_max) / (t_glace_max-t_glace_min))/(t_glace_max-t_glace_min) &
+                                  *exp(-distcltop(i)/dist_liq)
+                ENDIF
+        ENDIF
+
+        ! if temperature of cloud top <-40°C, 
+        IF (iflag_t_glace .GE. 4) THEN
+                IF ((temp_cltop(i) .LE. temp_nowater) .AND. (temp(i) .LE. t_glace_max)) THEN 
+                        icefrac(i) = 1.
+                        dicefracdT(i) = 0.
+                ENDIF 
+        ENDIF
+
+
+     ENDDO ! klon
+ 
+     RETURN
+ 
+END SUBROUTINE ICEFRAC_LSCP
+!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+
+
+SUBROUTINE CALC_QSAT_ECMWF(klon,temp,qtot,pressure,tref,phase,flagth,qs,dqs)
+!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+    ! Calculate qsat following ECMWF method
+!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+
+    IMPLICIT NONE
+
+    include "YOMCST.h"
+    include "YOETHF.h"
+    include "FCTTRE.h"
+
+    INTEGER, INTENT(IN) :: klon  ! number of horizontal grid points
+    REAL, INTENT(IN), DIMENSION(klon) :: temp     ! temperature in K
+    REAL, INTENT(IN), DIMENSION(klon) :: qtot     ! total specific water in kg/kg
+    REAL, INTENT(IN), DIMENSION(klon) :: pressure ! pressure in Pa
+    REAL, INTENT(IN)                  :: tref     ! reference temperature in K
+    LOGICAL, INTENT(IN) :: flagth     ! flag for qsat calculation for thermals
+    INTEGER, INTENT(IN) :: phase 
+    ! phase: 0=depend on temperature sign (temp>tref -> liquid, temp<tref, solid)
+    !        1=liquid
+    !        2=solid
+
+    REAL, INTENT(OUT), DIMENSION(klon) :: qs      ! saturation specific humidity [kg/kg]
+    REAL, INTENT(OUT), DIMENSION(klon) :: dqs     ! derivation of saturation specific humidity wrt T
+
+    REAL delta, cor, cvm5
+    INTEGER i
+
+    DO i=1,klon
+
+    IF (phase .EQ. 1) THEN
+        delta=0.
+    ELSEIF (phase .EQ. 2) THEN
+        delta=1.
+    ELSE
+        delta=MAX(0.,SIGN(1.,tref-temp(i)))
+    ENDIF
+
+    IF (flagth) THEN
+    cvm5=R5LES*(1.-delta) + R5IES*delta
+    ELSE
+    cvm5 = R5LES*RLVTT*(1.-delta) + R5IES*RLSTT*delta
+    cvm5 = cvm5 /RCPD/(1.0+RVTMP2*(qtot(i)))
+    ENDIF
+
+    qs(i)= R2ES*FOEEW(temp(i),delta)/pressure(i)
+    qs(i)=MIN(0.5,qs(i))
+    cor=1./(1.-RETV*qs(i))
+    qs(i)=qs(i)*cor
+    dqs(i)= FOEDE(temp(i),delta,cvm5,qs(i),cor)
+
+    END DO
+
+END SUBROUTINE CALC_QSAT_ECMWF
+!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+
+!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+SUBROUTINE CALC_GAMMASAT(klon,temp,qtot,pressure,gammasat,dgammasatdt)
+
+!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+    ! programme that calculates the gammasat parameter that determines the
+    ! homogeneous condensation thresholds for cold (<0oC) clouds
+    ! condensation at q>gammasat*qsat
+    ! Etienne Vignon, March 2021
+!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+    use lmdz_lscp_ini, only: iflag_gammasat, t_glace_min, RTT
+
+    IMPLICIT NONE
+
+
+    INTEGER, INTENT(IN) :: klon                       ! number of horizontal grid points
+    REAL, INTENT(IN), DIMENSION(klon) :: temp         ! temperature in K
+    REAL, INTENT(IN), DIMENSION(klon) :: qtot         ! total specific water in kg/kg
+
+    REAL, INTENT(IN), DIMENSION(klon) :: pressure     ! pressure in Pa
+
+    REAL, INTENT(OUT), DIMENSION(klon) :: gammasat    ! coefficient to multiply qsat with to calculate saturation
+    REAL, INTENT(OUT), DIMENSION(klon) :: dgammasatdt ! derivative of gammasat wrt temperature
+
+    REAL, DIMENSION(klon) ::  qsi,qsl,dqsl,dqsi
+    REAL  fcirrus, fac
+    REAL, PARAMETER :: acirrus=2.349
+    REAL, PARAMETER :: bcirrus=259.0
+
+    INTEGER i
+    
+        CALL CALC_QSAT_ECMWF(klon,temp,qtot,pressure,RTT,1,.false.,qsl,dqsl)
+        CALL CALC_QSAT_ECMWF(klon,temp,qtot,pressure,RTT,2,.false.,qsi,dqsi)
+
+    DO i=1,klon
+
+        IF (temp(i) .GE. RTT) THEN
+            ! warm clouds: condensation at saturation wrt liquid
+            gammasat(i)=1.
+            dgammasatdt(i)=0.
+
+        ELSEIF ((temp(i) .LT. RTT) .AND. (temp(i) .GT. t_glace_min)) THEN
+            
+            IF (iflag_gammasat .GE. 2) THEN          
+               gammasat(i)=qsl(i)/qsi(i)
+               dgammasatdt(i)=(dqsl(i)*qsi(i)-dqsi(i)*qsl(i))/qsi(i)/qsi(i)
+            ELSE
+               gammasat(i)=1.
+               dgammasatdt(i)=0.
+            ENDIF
+
+        ELSE
+
+            IF (iflag_gammasat .GE.1) THEN
+            ! homogeneous freezing of aerosols, according to
+            ! Koop, 2000 and Karcher 2008, QJRMS
+            ! 'Cirrus regime'
+               fcirrus=acirrus-temp(i)/bcirrus
+               IF (fcirrus .LT. qsl(i)/qsi(i)) THEN
+                  gammasat(i)=qsl(i)/qsi(i)
+                  dgammasatdt(i)=(dqsl(i)*qsi(i)-dqsi(i)*qsl(i))/qsi(i)/qsi(i)
+               ELSE
+                  gammasat(i)=fcirrus
+                  dgammasatdt(i)=-1.0/bcirrus
+               ENDIF
+           
+            ELSE
+
+               gammasat(i)=1.
+               dgammasatdt(i)=0.
+
+            ENDIF
+
+        ENDIF
+   
+    END DO
+
+
+END SUBROUTINE CALC_GAMMASAT
+!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+
+!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+SUBROUTINE DISTANCE_TO_CLOUD_TOP(klon,klev,k,temp,pplay,paprs,rneb,distcltop1D,temp_cltop)
+!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+   
+   USE lmdz_lscp_ini, ONLY : rd,rg,tresh_cl
+
+   IMPLICIT NONE
+   
+   INTEGER, INTENT(IN) :: klon,klev                !number of horizontal and vertical grid points
+   INTEGER, INTENT(IN) :: k                        ! vertical index
+   REAL, INTENT(IN), DIMENSION(klon,klev) :: temp  ! temperature in K
+   REAL, INTENT(IN), DIMENSION(klon,klev) :: pplay ! pressure middle layer in Pa
+   REAL, INTENT(IN), DIMENSION(klon,klev+1) :: paprs ! pressure interfaces in Pa
+   REAL, INTENT(IN), DIMENSION(klon,klev) :: rneb  ! cloud fraction
+
+   REAL, INTENT(OUT), DIMENSION(klon) :: distcltop1D  ! distance from cloud top
+   REAL, INTENT(OUT), DIMENSION(klon) :: temp_cltop     ! temperature of cloud top
+   
+   REAL dzlay(klon,klev)
+   REAL zlay(klon,klev)
+   REAL dzinterf
+   INTEGER i,k_top, kvert
+   LOGICAL bool_cl
+
+
+   DO i=1,klon
+         ! Initialization height middle of first layer
+          dzlay(i,1) = Rd * temp(i,1) / rg * log(paprs(i,1)/paprs(i,2))
+          zlay(i,1) = dzlay(i,1)/2
+
+          DO kvert=2,klev
+                 IF (kvert.EQ.klev) THEN
+                       dzlay(i,kvert) = 2*(rd * temp(i,kvert) / rg * log(paprs(i,kvert)/pplay(i,kvert)))
+                 ELSE
+                       dzlay(i,kvert) = rd * temp(i,kvert) / rg * log(paprs(i,kvert)/paprs(i,kvert+1))
+                 ENDIF
+                       dzinterf       = rd * temp(i,kvert) / rg * log(pplay(i,kvert-1)/pplay(i,kvert))
+                       zlay(i,kvert)  = zlay(i,kvert-1) + dzinterf
+           ENDDO
+   ENDDO
+   
+   DO i=1,klon
+          k_top = k
+          IF (rneb(i,k) .LE. tresh_cl) THEN
+                 bool_cl = .FALSE.
+          ELSE
+                 bool_cl = .TRUE.
+          ENDIF
+
+          DO WHILE ((bool_cl) .AND. (k_top .LE. klev))
+          ! find cloud top
+                IF (rneb(i,k_top) .GT. tresh_cl) THEN
+                      k_top = k_top + 1
+                ELSE
+                      bool_cl = .FALSE.
+                      k_top   = k_top - 1
+                ENDIF
+          ENDDO
+          k_top=min(k_top,klev)
+
+          !dist to top is dist between current layer and layer of cloud top (from middle to middle) + dist middle to
+          !interf for layer of cloud top
+          distcltop1D(i) = zlay(i,k_top) - zlay(i,k) + dzlay(i,k_top)/2
+          temp_cltop(i)  = temp(i,k_top)
+   ENDDO ! klon
+
+END SUBROUTINE DISTANCE_TO_CLOUD_TOP
+!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+END MODULE lmdz_lscp_tools
+
+
Index: DZ6/trunk/libf/phylmd/lscp_tools_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/lscp_tools_mod.F90	(revision 4664)
+++ 	(revision )
@@ -1,456 +1,0 @@
-MODULE LSCP_TOOLS_MOD
-
-    IMPLICIT NONE
-
-CONTAINS
-
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-SUBROUTINE FALLICE_VELOCITY(klon,iwc,temp,rho,pres,ptconv,velo)
-
-    ! Ref:
-    ! Stubenrauch, C. J., Bonazzola, M.,
-    ! Protopapadaki, S. E., & Musat, I. (2019).
-    ! New cloud system metrics to assess bulk
-    ! ice cloud schemes in a GCM. Journal of
-    ! Advances in Modeling Earth Systems, 11,
-    ! 3212–3234. https://doi.org/10.1029/2019MS001642
-    
-    use lscp_ini_mod, only: iflag_vice, ffallv_con, ffallv_lsc
-    use lscp_ini_mod, only: cice_velo, dice_velo 
-
-    IMPLICIT NONE
-
-    INTEGER, INTENT(IN) :: klon
-    REAL, INTENT(IN), DIMENSION(klon) :: iwc       ! specific ice water content [kg/m3]
-    REAL, INTENT(IN), DIMENSION(klon) :: temp      ! temperature [K]
-    REAL, INTENT(IN), DIMENSION(klon) :: rho       ! dry air density [kg/m3]
-    REAL, INTENT(IN), DIMENSION(klon) :: pres      ! air pressure [Pa]
-    LOGICAL, INTENT(IN), DIMENSION(klon) :: ptconv    ! convective point  [-]
-
-    REAL, INTENT(OUT), DIMENSION(klon) :: velo    ! fallspeed velocity of crystals [m/s]
-
-
-    INTEGER i
-    REAL logvm,iwcg,tempc,phpa,fallv_tun
-    REAL m2ice, m2snow, vmice, vmsnow
-    REAL aice, bice, asnow, bsnow
-    
-
-    DO i=1,klon
-
-        IF (ptconv(i)) THEN
-            fallv_tun=ffallv_con
-        ELSE
-            fallv_tun=ffallv_lsc
-        ENDIF
-
-        tempc=temp(i)-273.15 ! celcius temp
-        iwcg=MAX(iwc(i)*1000.,1E-3) ! iwc in g/m3. We set a minimum value to prevent from division by 0
-        phpa=pres(i)/100.    ! pressure in hPa
-
-    IF (iflag_vice .EQ. 1) THEN
-        ! so-called 'empirical parameterization' in Stubenrauch et al. 2019
-        if (tempc .GE. -60.0) then
-
-            logvm= -0.0000414122*tempc*tempc*log(iwcg)-0.00538922*tempc*log(iwcg) &
-                    -0.0516344*log(iwcg)+0.00216078*tempc + 1.9714    
-            velo(i)=exp(logvm)
-        else
-            velo(i)=65.0*(iwcg**0.2)*(150./phpa)**0.15
-        endif
-        
-        velo(i)=fallv_tun*velo(i)/100.0 ! from cm/s to m/s
-
-    ELSE IF (iflag_vice .EQ. 2) THEN
-        ! so called  PSDM empirical coherent bulk ice scheme in Stubenrauch et al. 2019
-        aice=0.587
-        bice=2.45
-        asnow=0.0444
-        bsnow=2.1
-        
-        m2ice=((iwcg*0.001/aice)/(exp(13.6-bice*7.76+0.479*bice**2)* & 
-                exp((-0.0361+bice*0.0151+0.00149*bice**2)*tempc)))   &
-                **(1./(0.807+bice*0.00581+0.0457*bice**2))
-
-        vmice=100.*1042.4*exp(13.6-(bice+1)*7.76+0.479*(bice+1.)**2)*exp((-0.0361+ &
-                 (bice+1.)*0.0151+0.00149*(bice+1.)**2)*tempc) &
-                 *(m2ice**(0.807+(bice+1.)*0.00581+0.0457*(bice+1.)**2))/(iwcg*0.001/aice)
-
-       
-        vmice=vmice*((1000./phpa)**0.2)
-      
-        m2snow=((iwcg*0.001/asnow)/(exp(13.6-bsnow*7.76+0.479*bsnow**2)* &
-               exp((-0.0361+bsnow*0.0151+0.00149*bsnow**2)*tempc)))         &
-               **(1./(0.807+bsnow*0.00581+0.0457*bsnow**2))
-
-
-        vmsnow=100.*14.3*exp(13.6-(bsnow+.416)*7.76+0.479*(bsnow+.416)**2)&
-                  *exp((-0.0361+(bsnow+.416)*0.0151+0.00149*(bsnow+.416)**2)*tempc)&
-                  *(m2snow**(0.807+(bsnow+.416)*0.00581+0.0457*(bsnow+.416)**2))/(iwcg*0.001/asnow)
-
-        vmsnow=vmsnow*((1000./phpa)**0.35)
-        velo(i)=fallv_tun*min(vmsnow,vmice)/100. ! to m/s
-
-    ELSE
-        ! By default, fallspeed velocity of ice crystals according to Heymsfield & Donner 1990
-        velo(i) = fallv_tun*cice_velo*((iwcg/1000.)**dice_velo)
-    ENDIF
-    ENDDO
-
-END SUBROUTINE FALLICE_VELOCITY
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-SUBROUTINE ICEFRAC_LSCP(klon, temp, iflag_ice_thermo, distcltop, temp_cltop, icefrac, dicefracdT)
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-  
-  ! Compute the ice fraction 1-xliq (see e.g.
-  ! Doutriaux-Boucher & Quaas 2004, section 2.2.)
-  ! as a function of temperature
-  ! see also Fig 3 of Madeleine et al. 2020, JAMES
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-
-
-    USE print_control_mod, ONLY: lunout, prt_level
-    USE lscp_ini_mod, ONLY: t_glace_min, t_glace_max, exposant_glace, iflag_t_glace
-    USE lscp_ini_mod, ONLY : RTT, dist_liq, temp_nowater
-
-    IMPLICIT NONE
-
-
-    INTEGER, INTENT(IN)                 :: klon              ! number of horizontal grid points
-    REAL, INTENT(IN), DIMENSION(klon)   :: temp              ! temperature
-    REAL, INTENT(IN), DIMENSION(klon)   :: distcltop         ! distance to cloud top
-    REAL, INTENT(IN), DIMENSION(klon)   :: temp_cltop        ! temperature of cloud top
-    INTEGER, INTENT(IN)                 :: iflag_ice_thermo
-    REAL, INTENT(OUT), DIMENSION(klon)  :: icefrac
-    REAL, INTENT(OUT), DIMENSION(klon)  :: dicefracdT
-
-
-    INTEGER i
-    REAL    liqfrac_tmp, dicefrac_tmp
-    REAL    Dv, denomdep,beta,qsi,dqsidt
-    LOGICAL ice_thermo
-
-    CHARACTER (len = 20) :: modname = 'lscp_tools'
-    CHARACTER (len = 80) :: abort_message
-
-    IF ((iflag_t_glace.LT.2) .OR. (iflag_t_glace.GT.6)) THEN
-       abort_message = 'lscp cannot be used if iflag_t_glace<2 or >6'
-       CALL abort_physic(modname,abort_message,1)
-    ENDIF
-
-    IF (.NOT.((iflag_ice_thermo .EQ. 1).OR.(iflag_ice_thermo .GE. 3))) THEN
-       abort_message = 'lscp cannot be used without ice thermodynamics'
-       CALL abort_physic(modname,abort_message,1)
-    ENDIF
-
-
-    DO i=1,klon
- 
-        ! old function with sole dependence upon temperature
-        IF (iflag_t_glace .EQ. 2) THEN
-            liqfrac_tmp = (temp(i)-t_glace_min) / (t_glace_max-t_glace_min)
-            liqfrac_tmp = MIN(MAX(liqfrac_tmp,0.0),1.0)
-            icefrac(i) = (1.0-liqfrac_tmp)**exposant_glace
-            IF (icefrac(i) .GT.0.) THEN
-                 dicefracdT(i)= exposant_glace * (icefrac(i)**(exposant_glace-1.)) &
-                           / (t_glace_min - t_glace_max)
-            ENDIF
-
-            IF ((icefrac(i).EQ.0).OR.(icefrac(i).EQ.1)) THEN
-                 dicefracdT(i)=0.
-            ENDIF
-
-        ENDIF
-
-        ! function of temperature used in CMIP6 physics
-        IF (iflag_t_glace .EQ. 3) THEN
-            liqfrac_tmp = (temp(i)-t_glace_min) / (t_glace_max-t_glace_min)
-            liqfrac_tmp = MIN(MAX(liqfrac_tmp,0.0),1.0)
-            icefrac(i) = 1.0-liqfrac_tmp**exposant_glace
-            IF ((icefrac(i) .GT.0.) .AND. (liqfrac_tmp .GT. 0.)) THEN
-                dicefracdT(i)= exposant_glace * ((liqfrac_tmp)**(exposant_glace-1.)) &
-                          / (t_glace_min - t_glace_max)
-            ELSE
-                dicefracdT(i)=0.
-            ENDIF
-        ENDIF
-
-        ! for iflag_t_glace .GE. 4, the liquid fraction depends upon temperature at cloud top
-        ! and then decreases with decreasing height 
-
-        !with linear function of temperature at cloud top
-        IF (iflag_t_glace .EQ. 4) THEN  
-                liqfrac_tmp = (temp(i)-t_glace_min) / (t_glace_max-t_glace_min)
-                liqfrac_tmp = MIN(MAX(liqfrac_tmp,0.0),1.0)
-                icefrac(i)    =  MAX(MIN(1.,1.0 - liqfrac_tmp*exp(-distcltop(i)/dist_liq)),0.)
-                dicefrac_tmp = - temp(i)/(t_glace_max-t_glace_min)
-                dicefracdT(i) = dicefrac_tmp*exp(-distcltop(i)/dist_liq)
-                IF ((liqfrac_tmp .LE.0) .OR. (liqfrac_tmp .GE. 1)) THEN
-                        dicefracdT(i) = 0.
-                ENDIF
-        ENDIF
-
-        ! with CMIP6 function of temperature at cloud top
-        IF (iflag_t_glace .EQ. 5) THEN 
-                liqfrac_tmp = (temp(i)-t_glace_min) / (t_glace_max-t_glace_min)
-                liqfrac_tmp =  MIN(MAX(liqfrac_tmp,0.0),1.0)
-                liqfrac_tmp = liqfrac_tmp**exposant_glace
-                icefrac(i)  =  MAX(MIN(1.,1.0 - liqfrac_tmp*exp(-distcltop(i)/dist_liq)),0.)
-                IF ((liqfrac_tmp .LE.0) .OR. (liqfrac_tmp .GE. 1)) THEN
-                        dicefracdT(i) = 0.
-                ELSE
-                        dicefracdT(i) = exposant_glace*((liqfrac_tmp)**(exposant_glace-1.))/(t_glace_min- t_glace_max) &
-                                        *exp(-distcltop(i)/dist_liq)
-                ENDIF
-        ENDIF
-
-        ! with modified function of temperature at cloud top 
-        ! to get largere values around 260 K, works well with t_glace_min = 241K
-        IF (iflag_t_glace .EQ. 6) THEN 
-                IF (temp(i) .GT. t_glace_max) THEN
-                        liqfrac_tmp = 1.
-                ELSE
-                        liqfrac_tmp = -((temp(i)-t_glace_max) / (t_glace_max-t_glace_min))**2+1.
-                ENDIF
-                liqfrac_tmp  = MIN(MAX(liqfrac_tmp,0.0),1.0)
-                icefrac(i)   = MAX(MIN(1.,1.0 - liqfrac_tmp*exp(-distcltop(i)/dist_liq)),0.)        
-                IF ((liqfrac_tmp .LE.0) .OR. (liqfrac_tmp .GE. 1)) THEN
-                        dicefracdT(i) = 0.
-                ELSE
-                        dicefracdT(i) = 2*((temp(i)-t_glace_max) / (t_glace_max-t_glace_min))/(t_glace_max-t_glace_min) &
-                                  *exp(-distcltop(i)/dist_liq)
-                ENDIF
-        ENDIF
-
-        ! if temperature of cloud top <-40°C, 
-        IF (iflag_t_glace .GE. 4) THEN
-                IF ((temp_cltop(i) .LE. temp_nowater) .AND. (temp(i) .LE. t_glace_max)) THEN 
-                        icefrac(i) = 1.
-                        dicefracdT(i) = 0.
-                ENDIF 
-        ENDIF
-
-
-     ENDDO ! klon
- 
-     RETURN
- 
-END SUBROUTINE ICEFRAC_LSCP
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-
-
-
-SUBROUTINE CALC_QSAT_ECMWF(klon,temp,qtot,pressure,tref,phase,flagth,qs,dqs)
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-    ! Calculate qsat following ECMWF method
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-
-
-    IMPLICIT NONE
-
-    include "YOMCST.h"
-    include "YOETHF.h"
-    include "FCTTRE.h"
-
-    INTEGER, INTENT(IN) :: klon  ! number of horizontal grid points
-    REAL, INTENT(IN), DIMENSION(klon) :: temp     ! temperature in K
-    REAL, INTENT(IN), DIMENSION(klon) :: qtot     ! total specific water in kg/kg
-    REAL, INTENT(IN), DIMENSION(klon) :: pressure ! pressure in Pa
-    REAL, INTENT(IN)                  :: tref     ! reference temperature in K
-    LOGICAL, INTENT(IN) :: flagth     ! flag for qsat calculation for thermals
-    INTEGER, INTENT(IN) :: phase 
-    ! phase: 0=depend on temperature sign (temp>tref -> liquid, temp<tref, solid)
-    !        1=liquid
-    !        2=solid
-
-    REAL, INTENT(OUT), DIMENSION(klon) :: qs      ! saturation specific humidity [kg/kg]
-    REAL, INTENT(OUT), DIMENSION(klon) :: dqs     ! derivation of saturation specific humidity wrt T
-
-    REAL delta, cor, cvm5
-    INTEGER i
-
-    DO i=1,klon
-
-    IF (phase .EQ. 1) THEN
-        delta=0.
-    ELSEIF (phase .EQ. 2) THEN
-        delta=1.
-    ELSE
-        delta=MAX(0.,SIGN(1.,tref-temp(i)))
-    ENDIF
-
-    IF (flagth) THEN
-    cvm5=R5LES*(1.-delta) + R5IES*delta
-    ELSE
-    cvm5 = R5LES*RLVTT*(1.-delta) + R5IES*RLSTT*delta
-    cvm5 = cvm5 /RCPD/(1.0+RVTMP2*(qtot(i)))
-    ENDIF
-
-    qs(i)= R2ES*FOEEW(temp(i),delta)/pressure(i)
-    qs(i)=MIN(0.5,qs(i))
-    cor=1./(1.-RETV*qs(i))
-    qs(i)=qs(i)*cor
-    dqs(i)= FOEDE(temp(i),delta,cvm5,qs(i),cor)
-
-    END DO
-
-END SUBROUTINE CALC_QSAT_ECMWF
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-
-
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-SUBROUTINE CALC_GAMMASAT(klon,temp,qtot,pressure,gammasat,dgammasatdt)
-
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-    ! programme that calculates the gammasat parameter that determines the
-    ! homogeneous condensation thresholds for cold (<0oC) clouds
-    ! condensation at q>gammasat*qsat
-    ! Etienne Vignon, March 2021
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-
-    use lscp_ini_mod, only: iflag_gammasat, t_glace_min, RTT
-
-    IMPLICIT NONE
-
-
-    INTEGER, INTENT(IN) :: klon                       ! number of horizontal grid points
-    REAL, INTENT(IN), DIMENSION(klon) :: temp         ! temperature in K
-    REAL, INTENT(IN), DIMENSION(klon) :: qtot         ! total specific water in kg/kg
-
-    REAL, INTENT(IN), DIMENSION(klon) :: pressure     ! pressure in Pa
-
-    REAL, INTENT(OUT), DIMENSION(klon) :: gammasat    ! coefficient to multiply qsat with to calculate saturation
-    REAL, INTENT(OUT), DIMENSION(klon) :: dgammasatdt ! derivative of gammasat wrt temperature
-
-    REAL, DIMENSION(klon) ::  qsi,qsl,dqsl,dqsi
-    REAL  fcirrus, fac
-    REAL, PARAMETER :: acirrus=2.349
-    REAL, PARAMETER :: bcirrus=259.0
-
-    INTEGER i
-    
-        CALL CALC_QSAT_ECMWF(klon,temp,qtot,pressure,RTT,1,.false.,qsl,dqsl)
-        CALL CALC_QSAT_ECMWF(klon,temp,qtot,pressure,RTT,2,.false.,qsi,dqsi)
-
-    DO i=1,klon
-
-        IF (temp(i) .GE. RTT) THEN
-            ! warm clouds: condensation at saturation wrt liquid
-            gammasat(i)=1.
-            dgammasatdt(i)=0.
-
-        ELSEIF ((temp(i) .LT. RTT) .AND. (temp(i) .GT. t_glace_min)) THEN
-            
-            IF (iflag_gammasat .GE. 2) THEN          
-               gammasat(i)=qsl(i)/qsi(i)
-               dgammasatdt(i)=(dqsl(i)*qsi(i)-dqsi(i)*qsl(i))/qsi(i)/qsi(i)
-            ELSE
-               gammasat(i)=1.
-               dgammasatdt(i)=0.
-            ENDIF
-
-        ELSE
-
-            IF (iflag_gammasat .GE.1) THEN
-            ! homogeneous freezing of aerosols, according to
-            ! Koop, 2000 and Karcher 2008, QJRMS
-            ! 'Cirrus regime'
-               fcirrus=acirrus-temp(i)/bcirrus
-               IF (fcirrus .LT. qsl(i)/qsi(i)) THEN
-                  gammasat(i)=qsl(i)/qsi(i)
-                  dgammasatdt(i)=(dqsl(i)*qsi(i)-dqsi(i)*qsl(i))/qsi(i)/qsi(i)
-               ELSE
-                  gammasat(i)=fcirrus
-                  dgammasatdt(i)=-1.0/bcirrus
-               ENDIF
-           
-            ELSE
-
-               gammasat(i)=1.
-               dgammasatdt(i)=0.
-
-            ENDIF
-
-        ENDIF
-   
-    END DO
-
-
-END SUBROUTINE CALC_GAMMASAT
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-
-
-!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-SUBROUTINE DISTANCE_TO_CLOUD_TOP(klon,klev,k,temp,pplay,paprs,rneb,distcltop1D,temp_cltop)
-!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-   
-   USE lscp_ini_mod, ONLY : rd,rg,tresh_cl
-
-   IMPLICIT NONE
-   
-   INTEGER, INTENT(IN) :: klon,klev                !number of horizontal and vertical grid points
-   INTEGER, INTENT(IN) :: k                        ! vertical index
-   REAL, INTENT(IN), DIMENSION(klon,klev) :: temp  ! temperature in K
-   REAL, INTENT(IN), DIMENSION(klon,klev) :: pplay ! pressure middle layer in Pa
-   REAL, INTENT(IN), DIMENSION(klon,klev+1) :: paprs ! pressure interfaces in Pa
-   REAL, INTENT(IN), DIMENSION(klon,klev) :: rneb  ! cloud fraction
-
-   REAL, INTENT(OUT), DIMENSION(klon) :: distcltop1D  ! distance from cloud top
-   REAL, INTENT(OUT), DIMENSION(klon) :: temp_cltop     ! temperature of cloud top
-   
-   REAL dzlay(klon,klev)
-   REAL zlay(klon,klev)
-   REAL dzinterf
-   INTEGER i,k_top, kvert
-   LOGICAL bool_cl
-
-
-   DO i=1,klon
-         ! Initialization height middle of first layer
-          dzlay(i,1) = Rd * temp(i,1) / rg * log(paprs(i,1)/paprs(i,2))
-          zlay(i,1) = dzlay(i,1)/2
-
-          DO kvert=2,klev
-                 IF (kvert.EQ.klev) THEN
-                       dzlay(i,kvert) = 2*(rd * temp(i,kvert) / rg * log(paprs(i,kvert)/pplay(i,kvert)))
-                 ELSE
-                       dzlay(i,kvert) = rd * temp(i,kvert) / rg * log(paprs(i,kvert)/paprs(i,kvert+1))
-                 ENDIF
-                       dzinterf       = rd * temp(i,kvert) / rg * log(pplay(i,kvert-1)/pplay(i,kvert))
-                       zlay(i,kvert)  = zlay(i,kvert-1) + dzinterf
-           ENDDO
-   ENDDO
-   
-   DO i=1,klon
-          k_top = k
-          IF (rneb(i,k) .LE. tresh_cl) THEN
-                 bool_cl = .FALSE.
-          ELSE
-                 bool_cl = .TRUE.
-          ENDIF
-
-          DO WHILE ((bool_cl) .AND. (k_top .LE. klev))
-          ! find cloud top
-                IF (rneb(i,k_top) .GT. tresh_cl) THEN
-                      k_top = k_top + 1
-                ELSE
-                      bool_cl = .FALSE.
-                      k_top   = k_top - 1
-                ENDIF
-          ENDDO
-          k_top=min(k_top,klev)
-
-          !dist to top is dist between current layer and layer of cloud top (from middle to middle) + dist middle to
-          !interf for layer of cloud top
-          distcltop1D(i) = zlay(i,k_top) - zlay(i,k) + dzlay(i,k_top)/2
-          temp_cltop(i)  = temp(i,k_top)
-   ENDDO ! klon
-
-END SUBROUTINE DISTANCE_TO_CLOUD_TOP
-!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-
-END MODULE LSCP_TOOLS_MOD
-
-
