From 602a524216b58b8a365a1b8d2f2431576bc5e12a Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Beno=C3=AEt=20Vi=C3=A9?= <benoit.vie@meteo.fr>
Date: Tue, 29 Mar 2022 11:57:45 +0200
Subject: [PATCH] New key LSNOW_T / Add J.Wurtz in routines update list / Add
 LSNOW_T option in lidar, radar etc

---
 src/MNH/aircraft_balloon_evol.f90        |  22 +++-
 src/MNH/default_desfmn.f90               |   3 +-
 src/MNH/ini_ice_c1r3.f90                 |  41 +++----
 src/MNH/ini_lima_cold_mixed.f90          | 143 +++++++++++------------
 src/MNH/ini_param_elec.f90               |   1 +
 src/MNH/ini_rain_ice.f90                 |   1 +
 src/MNH/ini_rain_ice_elec.f90            |   1 +
 src/MNH/init_aerosol_concentration.f90   |   6 +-
 src/MNH/lidar.f90                        |  25 +++-
 src/MNH/lima_cold_slow_processes.f90     |  55 ++++-----
 src/MNH/lima_conversion_melting_snow.f90 |  13 +--
 src/MNH/lima_droplets_riming_snow.f90    |   7 +-
 src/MNH/lima_graupel.f90                 |  20 ++--
 src/MNH/lima_ice_aggregation_snow.f90    |   7 +-
 src/MNH/lima_ice_snow_deposition.f90     |  11 +-
 src/MNH/lima_mixed.f90                   |  28 ++---
 src/MNH/lima_mixed_fast_processes.f90    |  70 ++---------
 src/MNH/lima_rain_accr_snow.f90          |  47 ++------
 src/MNH/lima_sedimentation.f90           |   5 +-
 src/MNH/lima_snow_deposition.f90         |  11 +-
 src/MNH/lima_tendencies.f90              |  28 +++--
 src/MNH/modd_param_lima.f90              |   1 +
 src/MNH/modd_param_lima_cold.f90         |   5 +-
 src/MNH/modd_param_lima_mixed.f90        |   1 +
 src/MNH/modn_param_lima.f90              |   2 +-
 src/MNH/profilern.f90                    |  22 +++-
 src/MNH/radar_rain_ice.f90               |  26 ++++-
 src/MNH/radar_scattering.f90             |  41 ++++---
 src/MNH/rrcolss.f90                      |   9 +-
 src/MNH/rscolrg.f90                      |  11 +-
 src/MNH/rzcolx.f90                       |  13 ++-
 src/MNH/write_lfifm1_for_diag.f90        |  16 +--
 32 files changed, 333 insertions(+), 359 deletions(-)

diff --git a/src/MNH/aircraft_balloon_evol.f90 b/src/MNH/aircraft_balloon_evol.f90
index 9642bcc4e..43cab9ac9 100644
--- a/src/MNH/aircraft_balloon_evol.f90
+++ b/src/MNH/aircraft_balloon_evol.f90
@@ -144,10 +144,11 @@ USE MODD_NSV,              ONLY : NSV_LIMA_NI,NSV_LIMA_NR,NSV_LIMA_NC
 USE MODD_PARAMETERS
 USE MODD_PARAM_LIMA,       ONLY: XALPHAR_L=>XALPHAR,XNUR_L=>XNUR,XALPHAS_L=>XALPHAS,XNUS_L=>XNUS,&
                                  XALPHAG_L=>XALPHAG,XNUG_L=>XNUG, XALPHAI_L=>XALPHAI,XNUI_L=>XNUI,&
-                                 XRTMIN_L=>XRTMIN,XALPHAC_L=>XALPHAC,XNUC_L=>XNUC
+                                 XRTMIN_L=>XRTMIN,XALPHAC_L=>XALPHAC,XNUC_L=>XNUC,LSNOW_T_L=>LSNOW_T
 USE MODD_PARAM_LIMA_COLD,  ONLY: XDI_L=>XDI,XLBEXI_L=>XLBEXI,XLBI_L=>XLBI,XAI_L=>XAI,XBI_L=>XBI,XC_I_L=>XC_I,&
                                  XLBEXS_L=>XLBEXS,XLBS_L=>XLBS,XCCS_L=>XCCS,&
-                                 XAS_L=>XAS,XBS_L=>XBS,XCXS_L=>XCXS
+                                 XAS_L=>XAS,XBS_L=>XBS,XCXS_L=>XCXS,        &
+                                 XLBDAS_MAX,XLBDAS_MIN
 USE MODD_PARAM_LIMA_MIXED, ONLY: XDG_L=>XDG,XLBEXG_L=>XLBEXG,XLBG_L=>XLBG,XCCG_L=>XCCG,&
                                  XAG_L=>XAG,XBG_L=>XBG,XCXG_L=>XCXG,XCG_L=>XCG
 USE MODD_PARAM_LIMA_WARM,  ONLY: XLBEXR_L=>XLBEXR,XLBR_L=>XLBR,XBR_L=>XBR,XAR_L=>XAR,&
@@ -318,6 +319,7 @@ REAL                           :: ZAETOT,ZAETMP,ZREFLOC,ZQSCA,ZQBACK,ZQEXT ! tem
 REAL,DIMENSION(:),ALLOCATABLE  :: ZAELOC,ZZMZ ! temporary arrays
 INTEGER                        :: JPTS_GAULAG=7 ! number of points for Gauss-Laguerre quadrature
 REAL                           :: ZLBDA   ! slope distribution parameter
+REAL                           :: ZN   ! number concentration
 REAL                           :: ZFRAC_ICE  ! ice water fraction
 REAL                           :: ZDELTA_EQUIV ! mass-equivalent Gauss-Laguerre point
 REAL                           :: ZFW ! liquid fraction
@@ -1177,7 +1179,17 @@ IF ( TPFLYER%FLY) THEN
                     ZLB=( ZA*ZCC*MOMG(ZALPHA,ZNU,ZB) )**(-ZLBEX)
                   ENDIF
               END SELECT
-              ZLBDA=ZLB*(ZRHODREFZ(JK)*ZRZ(JK,JLOOP))**ZLBEX
+              IF (JLOOP.EQ.5 .AND. CCLOUD=='LIMA' .AND. LSNOW_T_L) THEN
+                 IF (ZTEMPZ(JK)>-10.) THEN
+                    ZLBDA = MAX(MIN(XLBDAS_MAX, 10**(14.554-0.0423*(ZTEMPZ(JK)+273.15))),XLBDAS_MIN)
+                 ELSE
+                    ZLBDA = MAX(MIN(XLBDAS_MAX, 10**(6.226-0.0106*(ZTEMPZ(JK)+273.15))),XLBDAS_MIN)
+                 END IF
+                 ZN=XLBS_L*ZRHODREFZ(JK)*ZRZ(JK,JLOOP)*ZLBDA**ZB
+              ELSE
+                 ZLBDA=ZLB*(ZRHODREFZ(JK)*ZRZ(JK,JLOOP))**ZLBEX
+                 ZN=ZCC*ZLBDA**ZCX
+              END IF
               ZREFLOC=0.
               ZAETMP=0.
               DO JJ=1,JPTS_GAULAG ! Gauss-Laguerre quadrature
@@ -1209,8 +1221,8 @@ IF ( TPFLYER%FLY) THEN
                 ZREFLOC=ZREFLOC+ZQBACK*ZX(JJ)**(ZNU-1)*ZDELTA_EQUIV**2*ZW(JJ)
                 ZAETMP =ZAETMP +ZQEXT *ZX(JJ)**(ZNU-1)*ZDELTA_EQUIV**2*ZW(JJ)
               END DO
-              ZREFLOC=ZREFLOC*(XLAM_CRAD/XPI)**4*ZCC*ZLBDA**ZCX/(4.*GAMMA(ZNU)*.93)
-              ZAETMP=ZAETMP  *           XPI    *ZCC*ZLBDA**ZCX/(4.*GAMMA(ZNU))
+              ZREFLOC=ZREFLOC*(XLAM_CRAD/XPI)**4*ZN/(4.*GAMMA(ZNU)*.93)
+              ZAETMP=ZAETMP  *           XPI    *ZN/(4.*GAMMA(ZNU))
               TPFLYER%CRARE(IN,JK)=TPFLYER%CRARE(IN,JK)+ZREFLOC
               ZAELOC(JK)=ZAELOC(JK)+ZAETMP
             END IF
diff --git a/src/MNH/default_desfmn.f90 b/src/MNH/default_desfmn.f90
index c656a19af..76ccae576 100644
--- a/src/MNH/default_desfmn.f90
+++ b/src/MNH/default_desfmn.f90
@@ -276,7 +276,7 @@ USE MODD_ALLSTATION_n
 USE MODD_PARAM_LIMA, ONLY : LCOLD, LNUCL, LSEDI, LHHONI, LSNOW, LHAIL, LMEYERS,       &
                             NMOD_IFN, XIFN_CONC, LIFN_HOM, CIFN_SPECIES,              &
                             CINT_MIXING, NMOD_IMM, NIND_SPECIE,                       &
-                            CPRISTINE_ICE_LIMA, CHEVRIMED_ICE_LIMA,                   &
+                            LSNOW_T, CPRISTINE_ICE_LIMA, CHEVRIMED_ICE_LIMA,          &
                             XFACTNUC_DEP, XFACTNUC_CON,                               &
                             OWARM=>LWARM, LACTI, ORAIN=>LRAIN, OSEDC=>LSEDC,          &
                             OACTIT=>LACTIT, LBOUND, LSPRO, LADJ,                      &
@@ -1024,6 +1024,7 @@ IF (KMI == 1) THEN
   LSEDI  = .TRUE.
   LSNOW  = .TRUE.
   LHAIL  = .FALSE.
+  LSNOW_T = .TRUE.
   CPRISTINE_ICE_LIMA = 'PLAT'
   CHEVRIMED_ICE_LIMA = 'GRAU'
   XFACTNUC_DEP = 1.0  
diff --git a/src/MNH/ini_ice_c1r3.f90 b/src/MNH/ini_ice_c1r3.f90
index 771069f6b..883e56a1c 100644
--- a/src/MNH/ini_ice_c1r3.f90
+++ b/src/MNH/ini_ice_c1r3.f90
@@ -90,6 +90,7 @@ END MODULE MODI_INI_ICE_C1R3
 !!  Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
 !  P. Wautelet 10/04/2019: replace ABORT and STOP calls by Print_msg
 !  P. Wautelet 26/04/2019: replace non-standard FLOAT function by REAL function
+!  J. Wurtz       03/2022: new snow characteristics
 !
 !-------------------------------------------------------------------------------
 !
@@ -724,18 +725,18 @@ IF( (KACCLBDAS/=NACCLBDAS) .OR. (KACCLBDAR/=NACCLBDAR) .OR. (KND/=IND) .OR. &
     (PACCLBDAS_MAX/=XACCLBDAS_MAX) .OR. (PACCLBDAR_MAX/=XACCLBDAR_MAX) .OR. &
     (PACCLBDAS_MIN/=XACCLBDAS_MIN) .OR. (PACCLBDAR_MIN/=XACCLBDAR_MIN) .OR. &
     (PFDINFTY/=ZFDINFTY)                                               ) THEN
-!  CALL RRCOLSS ( IND, XALPHAS, XNUS, XALPHAR, XNUR,                          &
-!                 ZESR, XBR, XCS, XDS, XCR, XDR,                              &
-!                 XACCLBDAS_MAX, XACCLBDAR_MAX, XACCLBDAS_MIN, XACCLBDAR_MIN, &
-!                 ZFDINFTY, XKER_RACCSS, XAG, XBS, XAS                        )
-!  CALL RZCOLX  ( IND, XALPHAS, XNUS, XALPHAR, XNUR,                          &
-!                 ZESR, XBR, XCS, XDS, XCR, XDR,                              &
-!                 XACCLBDAS_MAX, XACCLBDAR_MAX, XACCLBDAS_MIN, XACCLBDAR_MIN, &
-!                 ZFDINFTY, XKER_RACCS                                        )
-!  CALL RSCOLRG ( IND, XALPHAS, XNUS, XALPHAR, XNUR,                          &
-!                 ZESR, XBS, XCS, XDS, XCR, XDR,                              &
-!                 XACCLBDAS_MAX, XACCLBDAR_MAX, XACCLBDAS_MIN, XACCLBDAR_MIN, &
-!                 ZFDINFTY, XKER_SACCRG,XAG, XBS, XAS                         )
+  CALL RRCOLSS ( IND, XALPHAS, XNUS, XALPHAR, XNUR,                          &
+                 ZESR, XBR, XCS, XDS, 0., XCR, XDR,                          &
+                 XACCLBDAS_MAX, XACCLBDAR_MAX, XACCLBDAS_MIN, XACCLBDAR_MIN, &
+                 ZFDINFTY, XKER_RACCSS, XAG, XBS, XAS                        )
+  CALL RZCOLX  ( IND, XALPHAS, XNUS, XALPHAR, XNUR,                          &
+                 ZESR, XBR, XCS, XDS, 0., XCR, XDR, 0.,                      &
+                 XACCLBDAS_MAX, XACCLBDAR_MAX, XACCLBDAS_MIN, XACCLBDAR_MIN, &
+                 ZFDINFTY, XKER_RACCS                                        )
+  CALL RSCOLRG ( IND, XALPHAS, XNUS, XALPHAR, XNUR,                          &
+                 ZESR, XBS, XCS, XDS, 0., XCR, XDR,                          &
+                 XACCLBDAS_MAX, XACCLBDAR_MAX, XACCLBDAS_MIN, XACCLBDAR_MIN, &
+                 ZFDINFTY, XKER_SACCRG,XAG, XBS, XAS                         )
   WRITE(UNIT=ILUOUT0,FMT='("*****************************************")')
   WRITE(UNIT=ILUOUT0,FMT='("**** UPDATE NEW SET OF RACSS KERNELS ****")')
   WRITE(UNIT=ILUOUT0,FMT='("**** UPDATE NEW SET OF RACS  KERNELS ****")')
@@ -937,10 +938,10 @@ IF( (KDRYLBDAG/=NDRYLBDAG) .OR. (KDRYLBDAS/=NDRYLBDAS) .OR. (KND/=IND) .OR. &
     (PDRYLBDAG_MAX/=XDRYLBDAG_MAX) .OR. (PDRYLBDAS_MAX/=XDRYLBDAS_MAX) .OR. &
     (PDRYLBDAG_MIN/=XDRYLBDAG_MIN) .OR. (PDRYLBDAS_MIN/=XDRYLBDAS_MIN) .OR. &
     (PFDINFTY/=ZFDINFTY)                                               ) THEN
-!  CALL RZCOLX ( IND, XALPHAG, XNUG, XALPHAS, XNUS,                          &
-!                ZEGS, XBS, XCG, XDG, XCS, XDS,                              &
-!                XDRYLBDAG_MAX, XDRYLBDAS_MAX, XDRYLBDAG_MIN, XDRYLBDAS_MIN, &
-!                ZFDINFTY, XKER_SDRYG                                        )
+  CALL RZCOLX ( IND, XALPHAG, XNUG, XALPHAS, XNUS,                          &
+                ZEGS, XBS, XCG, XDG, 0., XCS, XDS, 0.,                      &
+                XDRYLBDAG_MAX, XDRYLBDAS_MAX, XDRYLBDAG_MIN, XDRYLBDAS_MIN, &
+                ZFDINFTY, XKER_SDRYG                                        )
   WRITE(UNIT=ILUOUT0,FMT='("*****************************************")')
   WRITE(UNIT=ILUOUT0,FMT='("**** UPDATE NEW SET OF SDRYG KERNELS ****")')
   WRITE(UNIT=ILUOUT0,FMT='("*****************************************")')
@@ -1003,10 +1004,10 @@ IF( (KDRYLBDAG/=NDRYLBDAG) .OR. (KDRYLBDAR/=NDRYLBDAR) .OR. (KND/=IND) .OR. &
     (PDRYLBDAG_MAX/=XDRYLBDAG_MAX) .OR. (PDRYLBDAR_MAX/=XDRYLBDAR_MAX) .OR. &
     (PDRYLBDAG_MIN/=XDRYLBDAG_MIN) .OR. (PDRYLBDAR_MIN/=XDRYLBDAR_MIN) .OR. &
     (PFDINFTY/=ZFDINFTY)                                               ) THEN
-!  CALL RZCOLX ( IND, XALPHAG, XNUG, XALPHAR, XNUR,                          &
-!                ZEGR, XBR, XCG, XDG, XCR, XDR,                              &
-!                XDRYLBDAG_MAX, XDRYLBDAR_MAX, XDRYLBDAG_MIN, XDRYLBDAR_MIN, &
-!                ZFDINFTY, XKER_RDRYG                                        )
+  CALL RZCOLX ( IND, XALPHAG, XNUG, XALPHAR, XNUR,                          &
+                ZEGR, XBR, XCG, XDG, 0., XCR, XDR, 0.,                      &
+                XDRYLBDAG_MAX, XDRYLBDAR_MAX, XDRYLBDAG_MIN, XDRYLBDAR_MIN, &
+                ZFDINFTY, XKER_RDRYG                                        )
   WRITE(UNIT=ILUOUT0,FMT='("*****************************************")')
   WRITE(UNIT=ILUOUT0,FMT='("**** UPDATE NEW SET OF RDRYG KERNELS ****")')
   WRITE(UNIT=ILUOUT0,FMT='("*****************************************")')
diff --git a/src/MNH/ini_lima_cold_mixed.f90 b/src/MNH/ini_lima_cold_mixed.f90
index 993dee565..ea9b388c3 100644
--- a/src/MNH/ini_lima_cold_mixed.f90
+++ b/src/MNH/ini_lima_cold_mixed.f90
@@ -41,6 +41,7 @@ END MODULE MODI_INI_LIMA_COLD_MIXED
 !!  Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
 !  P. Wautelet 10/04/2019: replace ABORT and STOP calls by Print_msg
 !  P. Wautelet 26/04/2019: replace non-standard FLOAT function by REAL function
+!  J. Wurtz       03/2022: new snow characteristics
 !
 !-------------------------------------------------------------------------------
 !
@@ -172,25 +173,25 @@ XF1IS = 0.28
 XAS = 0.02
 XBS = 1.9
 
-!OPER
-!XCS = 5.      ! Wurtz
-!XDS = 0.27    ! Wurtz
-!XFVELOS = 0.  ! Wurtz
-
-!Cas GAMMAGEN
-XCS = 11.52     ! Wurtz
-XDS = 0.39      ! Wurtz
-XFVELOS =0.097  ! Wurtz
-
+IF (LSNOW_T) THEN
+!Cas Gamma generalisee
+XCS = 11.52
+XDS = 0.39
+XFVELOS =0.097
 !Cas MP
-!XCS = 13.2        ! Wurtz
-!XDS = 0.423       ! Wurtz        
-!XFVELOS = 25.14   ! Wurtz
+!XCS = 13.2
+!XDS = 0.423       
+!XFVELOS = 25.14
+ELSE
+XCS = 5.
+XDS = 0.27
+XFVELOS = 0.
+END IF
 
-!
-!XCCS = 5.0 ! Wurtz unused
-!XCXS = 1.0 ! Wurtz unused
-!
+IF (.NOT. LSNOW_T) THEN
+   XCCS = 5.0
+   XCXS = 1.0
+END IF
 
 XF0S = 0.86
 XF1S = 0.28
@@ -249,18 +250,17 @@ XC1H = 1./2.
 XALPHAI = 3.0  ! Gamma law for the ice crystal volume
 XNUI    = 3.0  ! Gamma law with little dispersion
 !
-!Wurtz
-!Cas MP
-!XALPHAS = 1.0  ! Exponential law
-!XNUS    = 1.0  ! Exponential law
-
-!Wurtz
+IF (LSNOW_T) THEN
 !Cas GAMMAGEN
-XALPHAS = .214   ! Generalized gamma law    !Wurtz
-XNUS    = 43.7   ! Generalized gamma law    !Wurtz
-XTRANS_MP_GAMMAS = SQRT( ( GAMMA(XNUS + 2./XALPHAS)*GAMMA(XNUS + 4./XALPHAS) ) / &
-                         ( 8.* GAMMA(XNUS + 1./XALPHAS)*GAMMA(XNUS + 3./XALPHAS) ) )    ! Wurtz
-
+   XALPHAS = .214   ! Generalized gamma law
+   XNUS    = 43.7   ! Generalized gamma law
+   XTRANS_MP_GAMMAS = SQRT( ( GAMMA(XNUS + 2./XALPHAS)*GAMMA(XNUS + 4./XALPHAS) ) / &
+                            ( 8.* GAMMA(XNUS + 1./XALPHAS)*GAMMA(XNUS + 3./XALPHAS) ) )
+ELSE
+   XALPHAS = 1.0  ! Exponential law
+   XNUS    = 1.0  ! Exponential law
+   XTRANS_MP_GAMMAS = 1.
+END IF
 !
 XALPHAG = 1.0  ! Exponential law
 XNUG    = 1.0  ! Exponential law
@@ -273,10 +273,13 @@ XNUH    = 8.0  ! Gamma law with little dispersion
 XLBEXI = 1.0/XBI
 XLBI   = XAI*MOMG(XALPHAI,XNUI,XBI)
 !
-!!XLBEXS = 1.0/(XCXS-XBS) ! Wurtz unused
-XLBEXS = 0. !1.0/(XCXS-XBS) ! Wurtz unused but used for XLB
-!XLBS   = ( XAS*XCCS*MOMG(XALPHAS,XNUS,XBS) )**(-XLBEXS)
-XLBS   = 1.0/(XAS*MOMG(XALPHAS,XNUS,XBS)) ! Wurtz
+IF (LSNOW_T) THEN
+   XLBEXS = 0. ! Not used
+   XLBS   = 1.0/(XAS*MOMG(XALPHAS,XNUS,XBS))
+ELSE
+   XLBEXS = 1.0/(XCXS-XBS)
+   XLBS   = ( XAS*XCCS*MOMG(XALPHAS,XNUS,XBS) )**(-XLBEXS)
+END IF
 !
 XLBEXG = 1.0/(XCXG-XBG)
 XLBG   = ( XAG*XCCG*MOMG(XALPHAG,XNUG,XBG))**(-XLBEXG)
@@ -340,21 +343,20 @@ XFSEDRI  = XC_I*GAMMA_X0D(XNUI+(XDI+XBI)/XALPHAI)/GAMMA_X0D(XNUI+XBI/XALPHAI)*
 XFSEDCI  = XC_I*GAMMA_X0D(XNUI+XDI/XALPHAI)/GAMMA_X0D(XNUI)*     &
             (ZRHO00)**XCEXVT
 !
-!OPER
-!XEXSEDS = (XBS+XDS-XCXS)/(XBS-XCXS)    ! Wurtz
-!XFSEDS  = XCS*XAS*XCCS*MOMG(XALPHAS,XNUS,XBS+XDS)*                         & ! Wurtz
-!            (XAS*XCCS*MOMG(XALPHAS,XNUS,XBS))**(-XEXSEDS)*(ZRHO00)**XCEXVT  ! Wurtz
-
+IF (LSNOW_T) THEN
 !HOUZE/HAIC
-!XEXSEDS = -XDS !(2*XBS+XDS) ! Modif Wurtz
-!XFSEDS  = XCS*MOMG(XALPHAS,XNUS,XBS+XDS)/(MOMG(XALPHAS,XNUS,XBS))    & ! Modif Wurtz
-!            *(ZRHO00)**XCEXVT
-
+   !XEXSEDS = -XDS !(2*XBS+XDS)
+   !XFSEDS  = XCS*MOMG(XALPHAS,XNUS,XBS+XDS)/(MOMG(XALPHAS,XNUS,XBS))    &
+   !            *(ZRHO00)**XCEXVT
 !LH_EXTENDED
-XEXSEDS = -XDS-XBS !-XNUS GAMMAGEN !(2*XBS+XDS) ! Modif Wurtz 
-XFSEDS  = XCS*MOMG(XALPHAS,XNUS,XBS+XDS)/(MOMG(XALPHAS,XNUS,XBS))    & ! Modif Wurtz
+   XEXSEDS = -XDS-XBS
+   XFSEDS  = XCS*MOMG(XALPHAS,XNUS,XBS+XDS)/(MOMG(XALPHAS,XNUS,XBS))    &
             *(ZRHO00)**XCEXVT
-!
+ELSE
+   XEXSEDS = (XBS+XDS-XCXS)/(XBS-XCXS)
+   XFSEDS  = XCS*XAS*XCCS*MOMG(XALPHAS,XNUS,XBS+XDS)*                         &
+            (XAS*XCCS*MOMG(XALPHAS,XNUS,XBS))**(-XEXSEDS)*(ZRHO00)**XCEXVT
+END IF
 !
 XEXSEDG = (XBG+XDG-XCXG)/(XBG-XCXG)
 XFSEDG  = XCG*XAG*XCCG*MOMG(XALPHAG,XNUG,XBG+XDG)*                         &
@@ -364,8 +366,6 @@ XEXSEDH = (XBH+XDH-XCXH)/(XBH-XCXH)
 XFSEDH  = XCH*XAH*XCCH*MOMG(XALPHAH,XNUH,XBH+XDH)*                         &
             (XAH*XCCH*MOMG(XALPHAH,XNUH,XBH))**(-XEXSEDH)*(ZRHO00)**XCEXVT
 !
-!
-!
 XLB(4)    = XLBI
 XLBEX(4)  = XLBEXI
 XD(4)     = XDI
@@ -601,7 +601,7 @@ XR1DEPIS = XC1DEPIS *(XAI*XDICNVS_LIM**XBI)
 !
 ! Harrington parameterization for snow to ice conversion
 !
-XLBDASCNVI_MAX = 6000.*XTRANS_MP_GAMMAS ! lbdas max after Field (1999) ! Wurtz
+XLBDASCNVI_MAX = 6000.*XTRANS_MP_GAMMAS ! lbdas max after Field (1999)
 !
 XDSCNVI_LIM    = 125.E-6  ! size in microns
 XLBDASCNVI_LIM = (50.0**(1.0/(XALPHAS)))/XDSCNVI_LIM  ! ZLBDAS Limitation
@@ -614,15 +614,10 @@ XR1DEPSI = XC1DEPSI *(XAS*XDSCNVI_LIM**XBS)
 !
 ! Vapor deposition on snow and graupel and hail
 !
-!X0DEPS = (4.0*XPI)*XCCS*XC1S*XF0S*MOMG(XALPHAS,XNUS,1.)
-!X1DEPS = (4.0*XPI)*XCCS*XC1S*XF1S*SQRT(XCS)*MOMG(XALPHAS,XNUS,0.5*XDS+1.5)
-X0DEPS = XLBS*(4.0*XPI)*XC1S*XF0S*MOMG(XALPHAS,XNUS,1.) ! Wurtz 
-X1DEPS = XLBS*(4.0*XPI)*XC1S*XF1S*SQRT(XCS)*MOMG(XALPHAS,XNUS,0.5*XDS+1.5) ! Wurtz 
-!XEX0DEPS = XCXS-1.0
-XEX0DEPS = XBS-1.0 ! Wurtz 
-!XEX1DEPS = XCXS-0.5*(XDS+3.0)
-XEX1DEPS = -0.5*(XDS+3.0) !GAMMAGEN !-0.5*(XDS+3.0)-XNUS ! Wurtz 
-
+X0DEPS = XLBS*(4.0*XPI)*XC1S*XF0S*MOMG(XALPHAS,XNUS,1.)
+X1DEPS = XLBS*(4.0*XPI)*XC1S*XF1S*SQRT(XCS)*MOMG(XALPHAS,XNUS,0.5*XDS+1.5)
+XEX0DEPS = XBS-1.0
+XEX1DEPS = -0.5*(XDS+3.0)
 !
 X0DEPG = (4.0*XPI)*XCCG*XC1G*XF0G*MOMG(XALPHAG,XNUG,1.)
 X1DEPG = (4.0*XPI)*XCCG*XC1G*XF1G*SQRT(XCG)*MOMG(XALPHAG,XNUG,0.5*XDG+1.5)
@@ -714,16 +709,13 @@ END IF
 !
 XDCSLIM  = 0.007 ! D_cs^lim = 7 mm as suggested by Farley et al. (1989)
 XCOLCS   = 1.0
-!XEXCRIMSS= XBS-XDS-2.0 ! Wurtz 
-XEXCRIMSS = -XDS-2.0 ! Wurtz GAMMAGEN
-!XCRIMSS  = (XPI/4.0)*XCOLCS*XCCS*XCS*(ZRHO00**XCEXVT)*MOMG(XALPHAS,XNUS,XDS+2.0)
-XCRIMSS  = XLBS * (XPI/4.0)*XCOLCS*XCS*(ZRHO00**XCEXVT)*MOMG(XALPHAS,XNUS,XDS+2.0) ! Wurtz
+XEXCRIMSS = -XDS-2.0
+XCRIMSS  = XLBS * (XPI/4.0)*XCOLCS*XCS*(ZRHO00**XCEXVT)*MOMG(XALPHAS,XNUS,XDS+2.0)
 
-XEXCRIMSG= XEXCRIMSS ! Done RHO
-XCRIMSG  = XCRIMSS   ! Done RHO
-!XSRIMCG  = XCCS*XAS*MOMG(XALPHAS,XNUS,XBS) 
-XSRIMCG  = XLBS*XAS*MOMG(XALPHAS,XNUS,XBS) ! Wurtz
-XEXSRIMCG=-XBS ! Wurtz GAMMAGEN
+XEXCRIMSG= XEXCRIMSS
+XCRIMSG  = XCRIMSS
+XSRIMCG  = XAS*MOMG(XALPHAS,XNUS,XBS)
+XEXSRIMCG= XCXS-XBS
 !
 GFLAG = .TRUE.
 IF (GFLAG) THEN
@@ -733,8 +725,8 @@ IF (GFLAG) THEN
 END IF
 !
 NGAMINC = 80
-XGAMINC_BOUND_MIN = (1000.*XTRANS_MP_GAMMAS*XDCSLIM)**XALPHAS !1.0E-1 ! Wurtz ! Minimal value of (Lbda * D_cs^lim)**alpha
-XGAMINC_BOUND_MAX = (50000.*XTRANS_MP_GAMMAS*XDCSLIM)**XALPHAS !1.0E7 ! Wurtz ! Maximal value of (Lbda * D_cs^lim)**alpha
+XGAMINC_BOUND_MIN = (1000.*XTRANS_MP_GAMMAS*XDCSLIM)**XALPHAS !1.0E-1 ! Minimal value of (Lbda * D_cs^lim)**alpha
+XGAMINC_BOUND_MAX = (50000.*XTRANS_MP_GAMMAS*XDCSLIM)**XALPHAS !1.0E7 ! Maximal value of (Lbda * D_cs^lim)**alpha
 ZRATE = EXP(LOG(XGAMINC_BOUND_MAX/XGAMINC_BOUND_MIN)/FLOAT(NGAMINC-1))
 !
 ALLOCATE( XGAMINC_RIM1(NGAMINC) )
@@ -781,15 +773,13 @@ XHMLINTP2 = 1.0 + XHMLINTP1*LOG( 25.E-6/(XGAMINC_HMC_BOUND_MIN)**(1.0/XALPHAC) )
 !
 !*       7.2    Constants for the accretion of raindrops onto aggregates
 !
-!XFRACCSS = ((XPI**2)/24.0)*XCCS*XCCR*XRHOLW*(ZRHO00**XCEXVT)
-XFRACCSS = XLBS*((XPI**2)/24.0)*XCCR*XRHOLW*(ZRHO00**XCEXVT) ! Wurtz 
+XFRACCSS = XLBS*((XPI**2)/24.0)*XCCR*XRHOLW*(ZRHO00**XCEXVT)
 !
 XLBRACCS1   =    MOMG(XALPHAS,XNUS,2.)*MOMG(XALPHAR,XNUR,3.)
 XLBRACCS2   = 2.*MOMG(XALPHAS,XNUS,1.)*MOMG(XALPHAR,XNUR,4.)
 XLBRACCS3   =                          MOMG(XALPHAR,XNUR,5.)
 !
-!XFSACCRG = (XPI/4.0)*XAS*XCCS*XCCR*(ZRHO00**XCEXVT)
-XFSACCRG = XLBS*(XPI/4.0)*XAS*XCCR*(ZRHO00**XCEXVT) ! Wurtz 
+XFSACCRG = XLBS*(XPI/4.0)*XAS*XCCR*(ZRHO00**XCEXVT)
 !
 XLBSACCR1   =    MOMG(XALPHAR,XNUR,2.)*MOMG(XALPHAS,XNUS,XBS)
 XLBSACCR2   = 2.*MOMG(XALPHAR,XNUR,1.)*MOMG(XALPHAS,XNUS,XBS+1.)
@@ -836,15 +826,15 @@ IF( (KACCLBDAS/=NACCLBDAS) .OR. (KACCLBDAR/=NACCLBDAR) .OR. (KND/=IND) .OR. &
     (PACCLBDAS_MIN/=XACCLBDAS_MIN) .OR. (PACCLBDAR_MIN/=XACCLBDAR_MIN) .OR. &
     (PFDINFTY/=ZFDINFTY)                                               ) THEN
   CALL RRCOLSS ( IND, XALPHAS, XNUS, XALPHAR, XNUR,                          &
-                 ZESR, XBR, XCS, XDS,XFVELOS, XCR, XDR,                      & ! Wurtz
+                 ZESR, XBR, XCS, XDS, XFVELOS, XCR, XDR,                     &
                  XACCLBDAS_MAX, XACCLBDAR_MAX, XACCLBDAS_MIN, XACCLBDAR_MIN, &
                  ZFDINFTY, XKER_RACCSS, XAG, XBS, XAS                        )
   CALL RZCOLX  ( IND, XALPHAS, XNUS, XALPHAR, XNUR,                          &
-                 ZESR, XBR, XCS, XDS, XFVELOS, XCR, XDR, 0.,                 & ! Wurtz
+                 ZESR, XBR, XCS, XDS, XFVELOS, XCR, XDR, 0.,                 &
                  XACCLBDAS_MAX, XACCLBDAR_MAX, XACCLBDAS_MIN, XACCLBDAR_MIN, &
                  ZFDINFTY, XKER_RACCS                                        )
   CALL RSCOLRG ( IND, XALPHAS, XNUS, XALPHAR, XNUR,                          &
-                 ZESR, XBS, XCS, XDS, XFVELOS, XCR, XDR,                     & ! Wurtz
+                 ZESR, XBS, XCS, XDS, XFVELOS, XCR, XDR,                     &
                  XACCLBDAS_MAX, XACCLBDAR_MAX, XACCLBDAS_MIN, XACCLBDAR_MIN, &
                  ZFDINFTY, XKER_SACCRG,XAG, XBS, XAS                         )
   WRITE(UNIT=ILUOUT0,FMT='("*****************************************")')
@@ -987,8 +977,7 @@ XCOLSG   = 0.01 ! Collection efficiency of S+G
 XCOLEXSG = 0.1  ! Temperature factor of the S+G collection efficiency
 WRITE (ILUOUT0, FMT=*) ' NEW Constants for the aggregate collection by the graupeln'
 WRITE (ILUOUT0, FMT=*) ' XCOLSG, XCOLEXSG  = ',XCOLSG,XCOLEXSG
-!XFSDRYG = (XPI/4.0)*XCOLSG*XCCG*XCCS*XAS*(ZRHO00**XCEXVT)
-XFSDRYG = XLBS*(XPI/4.0)*XCOLSG*XCCG*XAS*(ZRHO00**XCEXVT) ! Wurtz 
+XFSDRYG = XLBS*(XPI/4.0)*XCOLSG*XCCG*XAS*(ZRHO00**XCEXVT)
 !
 XLBSDRYG1   =    MOMG(XALPHAG,XNUG,2.)*MOMG(XALPHAS,XNUS,XBS)
 XLBSDRYG2   = 2.*MOMG(XALPHAG,XNUG,1.)*MOMG(XALPHAS,XNUS,XBS+1.)
@@ -1051,7 +1040,7 @@ IF( (KDRYLBDAG/=NDRYLBDAG) .OR. (KDRYLBDAS/=NDRYLBDAS) .OR. (KND/=IND) .OR. &
     (PDRYLBDAG_MIN/=XDRYLBDAG_MIN) .OR. (PDRYLBDAS_MIN/=XDRYLBDAS_MIN) .OR. &
     (PFDINFTY/=ZFDINFTY)                                               ) THEN
   CALL RZCOLX ( IND, XALPHAG, XNUG, XALPHAS, XNUS,                          &
-                ZEGS, XBS, XCG, XDG,0., XCS, XDS,XFVELOS,                   & !Wurtz
+                ZEGS, XBS, XCG, XDG, 0., XCS, XDS, XFVELOS,                 &
                 XDRYLBDAG_MAX, XDRYLBDAS_MAX, XDRYLBDAG_MIN, XDRYLBDAS_MIN, &
                 ZFDINFTY, XKER_SDRYG                                        )
   WRITE(UNIT=ILUOUT0,FMT='("*****************************************")')
diff --git a/src/MNH/ini_param_elec.f90 b/src/MNH/ini_param_elec.f90
index c1649f3f6..02a8b1578 100644
--- a/src/MNH/ini_param_elec.f90
+++ b/src/MNH/ini_param_elec.f90
@@ -85,6 +85,7 @@ END MODULE MODI_INI_PARAM_ELEC
 !!        J. Escobar 8/01/2016 bug , missing YDIR='XY' in READ 
 !!  Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
 !  P. Wautelet 26/04/2019: replace non-standard FLOAT function by REAL function
+!  J. Wurtz       03/2022: new snow characteristics
 !
 !-------------------------------------------------------------------------------
 !
diff --git a/src/MNH/ini_rain_ice.f90 b/src/MNH/ini_rain_ice.f90
index 754e25491..b89279296 100644
--- a/src/MNH/ini_rain_ice.f90
+++ b/src/MNH/ini_rain_ice.f90
@@ -102,6 +102,7 @@ END MODULE MODI_INI_RAIN_ICE
 !!      S. Riette 2016-11: new ICE3/ICE4 options
 !!      P. Wautelet 22/01/2019 bug correction: incorrect write
 !  P. Wautelet 26/04/2019: replace non-standard FLOAT function by REAL function
+!  J. Wurtz       03/2022: new snow characteristics
 !
 !-------------------------------------------------------------------------------
 !
diff --git a/src/MNH/ini_rain_ice_elec.f90 b/src/MNH/ini_rain_ice_elec.f90
index 6b132295f..d6ba287f1 100644
--- a/src/MNH/ini_rain_ice_elec.f90
+++ b/src/MNH/ini_rain_ice_elec.f90
@@ -87,6 +87,7 @@ END MODULE MODI_INI_RAIN_ICE_ELEC
 !!    Modifications:
 !!      C. Barthe   20/11/09   update to version 4.8.1
 !  P. Wautelet 26/04/2019: replace non-standard FLOAT function by REAL function
+!  J. Wurtz       03/2022: new snow characteristics
 !
 !-------------------------------------------------------------------------------
 !
diff --git a/src/MNH/init_aerosol_concentration.f90 b/src/MNH/init_aerosol_concentration.f90
index 90d16461f..32494739c 100644
--- a/src/MNH/init_aerosol_concentration.f90
+++ b/src/MNH/init_aerosol_concentration.f90
@@ -54,7 +54,7 @@ END MODULE MODI_INIT_AEROSOL_CONCENTRATION
 USE MODD_NSV
 USE MODD_PARAM_n,    ONLY : CCLOUD
 USE MODD_PARAM_LIMA, ONLY : LWARM, LACTI, NMOD_CCN, LSCAV, LAERO_MASS,      &
-                            XCCN_CONC, LCCN_HOM,                              &
+                            XCCN_CONC, LCCN_HOM,                            &
                             LCOLD, LNUCL, NMOD_IFN, LMEYERS,                &
                             XIFN_CONC, LIFN_HOM
 USE MODD_PARAMETERS, ONLY : JPVEXT
@@ -112,7 +112,7 @@ END IF ! LWARM AND LACTI
 ! Initialisation des concentrations en IFN
 !
 IF ( LCOLD .AND. LNUCL .AND. (.NOT. LMEYERS) ) THEN
-  DO JSV = NSV_LIMA_IFN_FREE, NSV_LIMA_IFN_NUCL+NMOD_IFN-1         
+  DO JSV = NSV_LIMA_IFN_FREE, NSV_LIMA_IFN_NUCL+NMOD_IFN-1
     PSVT(:,:,:,JSV) = 0.0
   ENDDO
   IKB = 1+JPVEXT
@@ -127,7 +127,7 @@ IF ( LCOLD .AND. LNUCL .AND. (.NOT. LMEYERS) ) THEN
   ELSE
 ! concentration décroissante selon z
     DO JSV = 1, NMOD_IFN
-      WHERE (PZZ(:,:,:) .LE. 1000.) 
+      WHERE (PZZ(:,:,:) .LE. 1000.)
          PSVT(:,:,:,NSV_LIMA_IFN_FREE+JSV-1) = XIFN_CONC(JSV)*1.0E3 / PRHODREF(:,:,:)
       ELSEWHERE (PZZ(:,:,:) .LE. 10000.)
          PSVT(:,:,:,NSV_LIMA_IFN_FREE+JSV-1) = XIFN_CONC(JSV)*1.0E3 &
diff --git a/src/MNH/lidar.f90 b/src/MNH/lidar.f90
index 4a3d987e8..1bee4acda 100644
--- a/src/MNH/lidar.f90
+++ b/src/MNH/lidar.f90
@@ -8,7 +8,7 @@
 !      #################
 !
 INTERFACE
-      SUBROUTINE LIDAR(HCLOUD,HVIEW,PALT,PWVL,PZZ,PRHO,PCLDFR,PRT, &
+      SUBROUTINE LIDAR(HCLOUD,HVIEW,PALT,PWVL,PZZ,PRHO,PT,PCLDFR,PRT, &
                        PLIDAROUT,PLIPAROUT,PCT,PDSTC,PDSTD,PDSTS)
 !
 CHARACTER(LEN=*),         INTENT(IN) :: HCLOUD  ! Name of the cloud scheme
@@ -17,6 +17,7 @@ REAL,                     INTENT(IN) :: PALT    ! Altitude of the lidar source
 REAL,                     INTENT(IN) :: PWVL    ! Wavelength of the lidar source
 REAL, DIMENSION(:,:,:),   INTENT(IN) :: PZZ     ! Altitude
 REAL, DIMENSION(:,:,:),   INTENT(IN) :: PRHO    ! Air density
+REAL, DIMENSION(:,:,:),   INTENT(IN) :: PT      ! Air temperature (C)
 REAL, DIMENSION(:,:,:),   INTENT(IN) :: PCLDFR  ! Cloud fraction
 REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PRT     ! Moist variables at t
 REAL, DIMENSION(:,:,:),  INTENT(OUT) :: PLIDAROUT ! Lidar output
@@ -36,7 +37,7 @@ END INTERFACE
 !
 END MODULE MODI_LIDAR
 !     #########################################################
-      SUBROUTINE LIDAR(HCLOUD,HVIEW,PALT,PWVL,PZZ,PRHO,PCLDFR,PRT, &
+      SUBROUTINE LIDAR(HCLOUD,HVIEW,PALT,PWVL,PZZ,PRHO,PT,PCLDFR,PRT, &
                        PLIDAROUT,PLIPAROUT,PCT,PDSTC,PDSTD,PDSTS)
 !     #########################################################
 !
@@ -109,9 +110,11 @@ USE MODD_ICE_C1R3_DESCR,  ONLY : XLBEXI,                      &
 USE MODD_PARAM_LIMA,      ONLY : URTMIN=>XRTMIN, UCTMIN=>XCTMIN, &
                                  UALPHAC=>XALPHAC,UNUC=>XNUC, &
                                  UALPHAR=>XALPHAR,UNUR=>XNUR, &
-                                 UALPHAI=>XALPHAI,UNUI=>XNUI 
+                                 UALPHAI=>XALPHAI,UNUI=>XNUI, &
+                                 USNOW_T=>LSNOW_T
 USE MODD_PARAM_LIMA_COLD, ONLY : UCCS=>XCCS, UCXS=>XCXS, ULBEXS=>XLBEXS, & 
-                                                         ULBS=>XLBS
+                                                         ULBS=>XLBS,     &
+                                 XLBDAS_MAX,XLBDAS_MIN,  UBS=>XBS
 USE MODD_PARAM_LIMA_MIXED,ONLY : UCCG=>XCCG, UCXG=>XCXG, ULBEXG=>XLBEXG, &
                                                          ULBG=>XLBG
 
@@ -130,6 +133,7 @@ REAL,                     INTENT(IN) :: PALT    ! Altitude of the lidar source
 REAL,                     INTENT(IN) :: PWVL    ! Wavelength of the lidar source
 REAL, DIMENSION(:,:,:),   INTENT(IN) :: PZZ     ! Altitude
 REAL, DIMENSION(:,:,:),   INTENT(IN) :: PRHO    ! Air density
+REAL, DIMENSION(:,:,:),   INTENT(IN) :: PT      ! Air temperature (C)
 REAL, DIMENSION(:,:,:),   INTENT(IN) :: PCLDFR  ! Cloud fraction
 REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PRT     ! Moist variables at t
 REAL, DIMENSION(:,:,:),  INTENT(OUT) :: PLIDAROUT ! Lidar output
@@ -523,9 +527,18 @@ SELECT CASE ( HCLOUD )
 !
             YDSD = 'MONOD'
             ZIWC    = PRHO(JI,JJ,JK)*PRT(JI,JJ,JK,5)
-            ZLBDAS  = ZLBS*(ZIWC)**ZLBEXS
+            IF (HCLOUD=='LIMA' .AND. USNOW_T) THEN
+               IF (PT(JI,JJ,JK)>-10.) THEN
+                  ZLBDAS = MAX(MIN(XLBDAS_MAX, 10**(14.554-0.0423*(PT(JI,JJ,JK)+273.15))),XLBDAS_MIN)
+               ELSE
+                  ZLBDAS = MAX(MIN(XLBDAS_MAX, 10**(6.226-0.0106*(PT(JI,JJ,JK)+273.15))),XLBDAS_MIN)
+               END IF
+               ZCONC=ZLBS*ZIWC*ZLBDAS**UBS
+            ELSE
+               ZLBDAS  = ZLBS*(ZIWC)**ZLBEXS
+               ZCONC   = ZCCS*(ZLBDAS)**ZCXS
+            END IF
             IF (ZLBDAS .GT. 0) THEN
-              ZCONC   = ZCCS*(ZLBDAS)**ZCXS
               ZRADIUS = 0.5*(3.0/ZLBDAS) ! Assume Marshall-Palmer law for Reff
               IANGLE  = 11
               CALL BHMIE_WATER( ZWAVE_LENGTH, ZZREFIND_ICE, YDSD, ZCONC,      &
diff --git a/src/MNH/lima_cold_slow_processes.f90 b/src/MNH/lima_cold_slow_processes.f90
index 5856b76d9..a95f91d46 100644
--- a/src/MNH/lima_cold_slow_processes.f90
+++ b/src/MNH/lima_cold_slow_processes.f90
@@ -80,6 +80,7 @@ END MODULE MODI_LIMA_COLD_SLOW_PROCESSES
 !  P. Wautelet 05/2016-04/2018: new data structures and calls for I/O
 !  P. Wautelet 28/05/2019: move COUNTJV function to tools.f90
 !  P. Wautelet    03/2020: use the new data structures and subroutines for budgets
+!  J. Wurtz       03/2022: new snow characteristics
 !
 !-------------------------------------------------------------------------------
 !
@@ -94,8 +95,8 @@ USE MODD_CST,             ONLY: XP00, XRD, XRV, XMV, XMD, XCPD, XCPV,        &
                                 XCL, XCI, XTT, XLSTT, XALPI, XBETAI, XGAMI
 USE MODD_NSV,             ONLY: NSV_LIMA_NI
 USE MODD_PARAMETERS,      ONLY: JPHEXT, JPVEXT
-USE MODD_PARAM_LIMA,      ONLY: LSNOW, XRTMIN, XCTMIN, XALPHAI, XALPHAS,     &
-                                XNUI, XNUS
+USE MODD_PARAM_LIMA,      ONLY: LSNOW, LSNOW_T, XRTMIN, XCTMIN,              &
+                                XALPHAI, XALPHAS, XNUI, XNUS
 USE MODD_PARAM_LIMA_COLD, ONLY: XLBI, XLBEXI, XLBS, XLBEXS, XBI, XCXS, XCCS, &
                                 XLBDAS_MAX, XDSCNVI_LIM, XLBDASCNVI_MAX,     &
                                 XC0DEPSI, XC1DEPSI, XR0DEPSI, XR1DEPSI,      &
@@ -317,16 +318,20 @@ IF( IMICRO >= 1 ) THEN
       ZLBDAI(:) = ( XLBI*ZCIT(:) / ZRIT(:) )**XLBEXI
    END WHERE
    ZLBDAS(:)  = 1.E10
-   WHERE (ZRST(:)>XRTMIN(5) )
-WHERE(ZZT(:)>263.15) 
-ZLBDAS(:) = MAX(MIN(XLBDAS_MAX, 10**(14.554-0.0423*ZZT(:))),XLBDAS_MIN)
-END WHERE
-WHERE(ZZT(:)<=263.15) 
-ZLBDAS(:) = MAX(MIN(XLBDAS_MAX, 10**(6.226-0.0106*ZZT(:))),XLBDAS_MIN)
-END WHERE
-   END WHERE
-ZLBDAS(:) = ZLBDAS(:) * XTRANS_MP_GAMMAS
-
+   IF (LSNOW_T) THEN
+      WHERE(ZZT(:)>263.15 .AND. ZRST(:)>XRTMIN(5)) 
+         ZLBDAS(:) = MAX(MIN(XLBDAS_MAX, 10**(14.554-0.0423*ZZT(:))),XLBDAS_MIN)
+      END WHERE
+      WHERE(ZZT(:)<=263.15 .AND. ZRST(:)>XRTMIN(5)) 
+         ZLBDAS(:) = MAX(MIN(XLBDAS_MAX, 10**(6.226-0.0106*ZZT(:))),XLBDAS_MIN)
+      END WHERE
+      ZLBDAS(:) = ZLBDAS(:) * XTRANS_MP_GAMMAS
+   ELSE
+      WHERE (ZRST(:)>XRTMIN(5) )
+         ZLBDAS(:) = MAX(MIN(XLBDAS_MAX,XLBS*( ZRHODREF(:)*ZRST(:) )**XLBEXS),XLBDAS_MIN)
+      END WHERE
+   END IF
+!
    ZKA(:) = 2.38E-2 + 0.0071E-2 * ( ZZT(:) - XTT )          ! k_a
    ZDV(:) = 0.211E-4 * (ZZT(:)/XTT)**1.94 * (XP00/ZPRES(:)) ! D_v
 !
@@ -349,21 +354,11 @@ ZLBDAS(:) = ZLBDAS(:) * XTRANS_MP_GAMMAS
           call Budget_store_init( tbudgets(NBUDGET_SV1 - 1 + nsv_lima_ni), 'CNVI', pcis(:, :, :) * prhodj(:, :, :) )
       end if
 
-!      WHERE ( ZRST(:)>XRTMIN(5) )
-!         ZLBDAS(:)  = MIN( XLBDAS_MAX,                                           &
-!                           XLBS*( ZRHODREF(:)*MAX( ZRST(:),XRTMIN(5) ) )**XLBEXS )
-!      END WHERE
       ZZW(:) = 0.0
       WHERE ( ZLBDAS(:)<XLBDASCNVI_MAX .AND. (ZRST(:)>XRTMIN(5)) &
                                        .AND. (ZSSI(:)<0.0)       )
          ZZW(:) = (ZLBDAS(:)*XDSCNVI_LIM)**(XALPHAS)
-       !  ZZX(:) = ( -ZSSI(:)/ZAI(:) ) * (XCCS*ZLBDAS(:)**XCXS) * (ZZW(:)**XNUI) &
-                                                              ! * EXP(-ZZW(:))
-         ZZX(:) = ( -ZSSI(:)/ZAI(:) ) * (XLBS*ZRST(:)*ZLBDAS(:)**XBS) * (ZZW(:)**XNUI)  &
-
-                                                              * EXP(-ZZW(:))
-!MODIF CNVI -> concentration de la neige en kg-1 et non en m-3
-
+         ZZX(:) = ( -ZSSI(:)/ZAI(:) ) * (XLBS*ZRST(:)*ZLBDAS(:)**XBS) * (ZZW(:)**XNUI) * EXP(-ZZW(:))
 !
          ZZW(:) = MIN( ( XR0DEPSI+XR1DEPSI*ZCJ(:) )*ZZX(:),ZRSS(:) )
          ZRIS(:) = ZRIS(:) + ZZW(:)
@@ -396,14 +391,10 @@ ZLBDAS(:) = ZLBDAS(:) * XTRANS_MP_GAMMAS
 
       ZZW(:) = 0.0
       WHERE ( (ZRST(:)>XRTMIN(5)) .AND. (ZRSS(:)>ZRTMIN(5)) )
-!Correction BVIE rhodref
-!         ZZW(:) = ( ZSSI(:)/(ZRHODREF(:)*ZAI(:)) ) *                               &
-     !    ZZW(:) = ( ZSSI(:)/(ZAI(:)) ) *                               &
-      !            ( X0DEPS*ZLBDAS(:)**XEX0DEPS + X1DEPS*ZCJ(:)*ZLBDAS(:)**XEX1DEPS )
-         ZZW(:) = ( ZRHODREF(:)*ZRST(:)*ZSSI(:)/(ZAI(:)) ) *                                & ! Wurtz
+         ZZW(:) = ( ZRST(:)*ZSSI(:)/(ZAI(:)) ) *                                            &
               ( X0DEPS*ZLBDAS(:)**XEX0DEPS +                                                &
               (X1DEPS*ZCJ(:)*(1+(XFVELOS/(2.*ZLBDAS))**XALPHAS)**(-XNUS+XEX1DEPS/XALPHAS) * &
-              (ZLBDAS(:))**(XEX1DEPS+XBS))) ! Wurtz
+              (ZLBDAS(:))**(XEX1DEPS+XBS)))
 
          ZZW(:) =    MIN( ZRVS(:),ZZW(:)      )*(0.5+SIGN(0.5,ZZW(:))) &
                    - MIN( ZRSS(:),ABS(ZZW(:)) )*(0.5-SIGN(0.5,ZZW(:)))
@@ -439,8 +430,6 @@ ZLBDAS(:) = ZLBDAS(:) * XTRANS_MP_GAMMAS
          ZZW(:) = (ZLBDAI(:)*XDICNVS_LIM)**(XALPHAI)
          ZZX(:) = ( ZSSI(:)/ZAI(:) )*ZCIT(:) * (ZZW(:)**XNUI) *EXP(-ZZW(:))
 !
-! Correction BVIE
-!         ZZW(:) = MAX( MIN( ( XR0DEPIS + XR1DEPIS*ZCJ(:) )*ZZX(:)/ZRHODREF(:) &
          ZZW(:) = MAX( MIN( ( XR0DEPIS + XR1DEPIS*ZCJ(:) )*ZZX(:) &
                             ,ZRIS(:) ) + ZRTMIN(5), ZRTMIN(5) ) - ZRTMIN(5)
          ZRIS(:) = ZRIS(:) - ZZW(:)
@@ -477,10 +466,6 @@ ZLBDAS(:) = ZLBDAS(:) * XTRANS_MP_GAMMAS
       WHERE ( (ZRIT(:)>XRTMIN(4)) .AND. (ZRST(:)>XRTMIN(5)) .AND. (ZRIS(:)>ZRTMIN(4)) &
                                                             .AND. (ZCIS(:)>ZCTMIN(4)) )
          ZZW1(:,3) = (ZLBDAI(:) / ZLBDAS(:))**3
-!         ZZW1(:,1) = (ZCIT(:)*(XCCS*ZLBDAS(:)**XCXS)*EXP( XCOLEXIS*(ZZT(:)-XTT) )) &
-!                                                    / (ZLBDAI(:)**3)
-!         ZZW1(:,1) = (ZCIT(:)*(XLBS*ZRHODREF(:)*ZRST(:)*ZLBDAS(:)**XBS)*EXP(XCOLEXIS*(ZZT(:)-XTT) )) &
-!                                                    / (ZLBDAI(:)**3)
          ZZW1(:,1) = (ZCIT(:)*(XLBS*ZRST(:)*ZLBDAS(:)**XBS)*EXP(XCOLEXIS*(ZZT(:)-XTT) )) &
                                                     / (ZLBDAI(:)**3)
          ZZW1(:,2) = MIN( ZZW1(:,1)*(XAGGS_CLARGE1+XAGGS_CLARGE2*ZZW1(:,3)),ZCIS(:) )
diff --git a/src/MNH/lima_conversion_melting_snow.f90 b/src/MNH/lima_conversion_melting_snow.f90
index b59c6f08c..307db0255 100644
--- a/src/MNH/lima_conversion_melting_snow.f90
+++ b/src/MNH/lima_conversion_melting_snow.f90
@@ -55,6 +55,7 @@ END MODULE MODI_LIMA_CONVERSION_MELTING_SNOW
 !!    -------------
 !!      Original             15/03/2018 
 !!
+!  J. Wurtz       03/2022: new snow characteristics
 !-------------------------------------------------------------------------------
 !
 !*       0.    DECLARATIONS
@@ -106,14 +107,10 @@ WHERE( (PRST(:)>XRTMIN(5)) .AND. (PT(:)>XTT) .AND. LDCOMPUTE(:) )
 !
 ! compute RSMLT
 !
-!   ZW(:)  = XFSCVMG*MAX( 0.0,( -ZW(:) *             &
-!                          ( X0DEPS*           PLBDS(:)**XEX0DEPS +     &
-!                            X1DEPS*PCJ(:)*PLBDS(:)**XEX1DEPS ) ))!-    &
-
-   ZW(:)  = XFSCVMG*MAX( 0.0,( -ZW(:) *             &   ! Wurtz
-                               PRHODREF(:)*PRST(:)*( X0DEPS*            PLBDS(:)**XEX0DEPS +     &
-                               X1DEPS*PCJ(:)*(1+(XFVELOS/(2.*PLBDS(:)))**XALPHAS)                &
-                               **(-XNUS+XEX1DEPS/XALPHAS)*(PLBDS(:))**(XEX1DEPS+XBS))))
+   ZW(:)  = XFSCVMG*MAX( 0.0,( -ZW(:) * PRHODREF(:) * PRST(:) *          &
+                               ( X0DEPS*PLBDS(:)**XEX0DEPS +             &
+                                 X1DEPS*PCJ(:)*PLBDS(:)**(XEX1DEPS+XBS)* &
+                                   (1+(XFVELOS/(2.*PLBDS(:)))**XALPHAS)**(-XNUS+XEX1DEPS/XALPHAS)) ))
 ! On ne tient pas compte de la collection de pluie et gouttelettes par la neige si T>0 !!!! 
 ! Note that no heat is exchanged because the graupeln produced are still icy!!!
    P_RS_CMEL(:) = - ZW(:)
diff --git a/src/MNH/lima_droplets_riming_snow.f90 b/src/MNH/lima_droplets_riming_snow.f90
index 132516a76..b1c4a8007 100644
--- a/src/MNH/lima_droplets_riming_snow.f90
+++ b/src/MNH/lima_droplets_riming_snow.f90
@@ -66,6 +66,7 @@ END MODULE MODI_LIMA_DROPLETS_RIMING_SNOW
 !!    -------------
 !!      Original             15/03/2018 
 !  P. Wautelet 26/04/2019: replace non-standard FLOAT function by REAL function
+!  J. Wurtz       03/2022: new snow characteristics
 !
 !-------------------------------------------------------------------------------
 !
@@ -171,10 +172,9 @@ WHERE( GRIM )
 !        4.     riming
 !
    ! Cloud droplets collected
-!   P_RC_RIM(:) = - XCRIMSS * PRCT(:) * PLBDS(:)**XEXCRIMSS * PRHODREF(:)**(-XCEXVT)
    P_RC_RIM(:) = - XCRIMSS  * PRCT(:) * PRST(:)*(1+(XFVELOS/PLBDS(:))**XALPHAS)**(-XNUS+XEXCRIMSS/XALPHAS) &
                                       * PRHODREF(:)**(-XCEXVT+1) &
-                                      * (PLBDS(:)) ** (XEXCRIMSS+XBS) ! Wurtz
+                                      * (PLBDS(:)) ** (XEXCRIMSS+XBS)
    P_CC_RIM(:) = P_RC_RIM(:) *(PCCT(:)/PRCT(:)) ! Lambda_c**3
    !
    ! Cloud droplets collected on small aggregates add to snow
@@ -184,8 +184,7 @@ WHERE( GRIM )
    P_RG_RIM(:) = - P_RC_RIM(:) - P_RS_RIM(:) 
    !
    ! Large aggregates collecting droplets add to graupel (instant process ???)
-!   ZZW3(:) = XSRIMCG * PLBDS(:)**XEXSRIMCG * (1.0 - ZZW2(:))/(PTSTEP*PRHODREF(:))
-   ZZW3(:) = XSRIMCG *PRHODREF(:)*PRST(:)* PLBDS(:)**(XEXSRIMCG+XBS) * (1.0 - ZZW2(:))!/(PTSTEP*PRHODREF(:)) ! Wurtz
+   ZZW3(:) = PRST(:)*(1.0 - ZZW2(:))/PTSTEP
    P_RS_RIM(:) = P_RS_RIM(:) - ZZW3(:) 
    P_RG_RIM(:) = P_RG_RIM(:) + ZZW3(:) 
    !
diff --git a/src/MNH/lima_graupel.f90 b/src/MNH/lima_graupel.f90
index 6806322cd..578f22bbf 100644
--- a/src/MNH/lima_graupel.f90
+++ b/src/MNH/lima_graupel.f90
@@ -129,6 +129,7 @@ END MODULE MODI_LIMA_GRAUPEL
 !!    -------------
 !!      Original             15/03/2018
 !  P. Wautelet 26/04/2019: replace non-standard FLOAT function by REAL function
+!  J. Wurtz       03/2022: new snow characteristics
 !
 !-------------------------------------------------------------------------------
 !
@@ -330,19 +331,12 @@ WHERE( GDRY )
        			                                    * (ZVEC1(:) - 1.0)
    ZZW(:) = ZVEC3(:)
 !
-!   ZZW3(:) = XFSDRYG * ZZW(:) * EXP( XCOLEXSG*(PT(:)-XTT) )  & ! RSDRYG - rs collected by graupel in dry mode
-!                    *( PLBDS(:)**(XCXS-XBS) )*( PLBDG(:)**XCXG )    &
-!                    *( PRHODREF(:)**(-XCEXVT-1.) )                      &
-!                    *( XLBSDRYG1/( PLBDG(:)**2                 ) + &
-!                       XLBSDRYG2/( PLBDG(:)   * PLBDS(:)   ) + &
-!                       XLBSDRYG3/(                  PLBDS(:)**2) )
-! Wurtz 
-   ZZW3(:) = XFSDRYG * ZZW(:) * EXP( XCOLEXSG*(PT(:)-XTT) )  & ! RSDRYG - rs collected by graupel in dry mode
-                    *( PRST(:))*( PLBDG(:)**XCXG )    &
-                    *( PRHODREF(:)**(-XCEXVT) )                      &
-                    *( XLBSDRYG1/( PLBDG(:)**2                 ) + &
-                       XLBSDRYG2/( PLBDG(:)   * PLBDS(:)   ) + &
-                       XLBSDRYG3/(                  PLBDS(:)**2) )
+   ZZW3(:) = XFSDRYG * ZZW(:) * EXP( XCOLEXSG*(PT(:)-XTT) )       & ! RSDRYG - rs collected by graupel in dry mode
+                    *( PRST(:))*( PLBDG(:)**XCXG )                &
+                    *( PRHODREF(:)**(-XCEXVT) )                   &
+                    *( XLBSDRYG1/( PLBDG(:)**2                ) + &
+                       XLBSDRYG2/( PLBDG(:)    * PLBDS(:)     ) + &
+                       XLBSDRYG3/(               PLBDS(:)**2) )
 END WHERE
 !
 !*           1.c  Collection of rr in the dry mode
diff --git a/src/MNH/lima_ice_aggregation_snow.f90 b/src/MNH/lima_ice_aggregation_snow.f90
index 6c7be6eee..24f74dca1 100644
--- a/src/MNH/lima_ice_aggregation_snow.f90
+++ b/src/MNH/lima_ice_aggregation_snow.f90
@@ -52,6 +52,7 @@ END MODULE MODI_LIMA_ICE_AGGREGATION_SNOW
 !!    MODIFICATIONS
 !!    -------------
 !!      Original             15/03/2018
+!  J. Wurtz       03/2022: new snow characteristics
 !!
 !-------------------------------------------------------------------------------
 !
@@ -101,10 +102,8 @@ P_CI_AGGS(:) = 0.
 !
 WHERE ( (PRIT(:)>XRTMIN(4)) .AND. (PRST(:)>XRTMIN(5)) .AND. LDCOMPUTE(:) )
    ZZW1(:) = (PLBDI(:) / PLBDS(:))**3
-!   ZZW2(:) = (PCIT(:)*(XCCS*PLBDS(:)**XCXS)/PRHODREF(:)*EXP( XCOLEXIS*(PT(:)-XTT) )) &
-!        / (PLBDI(:)**3)
-   ZZW2(:) = (PCIT(:)*(XLBS*PRST(:)*PLBDS(:)**XBS)*EXP(XCOLEXIS*(PT(:)-XTT) )) & ! Wurtz
-        / (PLBDI(:)**3)                                                          ! Wurtz
+   ZZW2(:) = (PCIT(:)*(XLBS*PRST(:)*PLBDS(:)**XBS)*EXP(XCOLEXIS*(PT(:)-XTT) )) &
+        / (PLBDI(:)**3)
    ZZW3(:) = ZZW2(:)*(XAGGS_CLARGE1+XAGGS_CLARGE2*ZZW1(:))
 !
    P_CI_AGGS(:) = - ZZW3(:)
diff --git a/src/MNH/lima_ice_snow_deposition.f90 b/src/MNH/lima_ice_snow_deposition.f90
index fd0261dc1..f31d3175d 100644
--- a/src/MNH/lima_ice_snow_deposition.f90
+++ b/src/MNH/lima_ice_snow_deposition.f90
@@ -78,6 +78,7 @@ SUBROUTINE LIMA_ICE_SNOW_DEPOSITION (PTSTEP, LDCOMPUTE,                        &
 !!    MODIFICATIONS
 !!    -------------
 !!      Original             15/03/2018
+!  J. Wurtz       03/2022: new snow characteristics
 !!
 !-------------------------------------------------------------------------------
 !
@@ -166,7 +167,6 @@ WHERE( GMICRO )
    WHERE ( PLBDS(:)<XLBDASCNVI_MAX .AND. (PRST(:)>XRTMIN(5)) &
                                    .AND. (PSSI(:)<0.0)       )
       ZZW(:) = (PLBDS(:)*XDSCNVI_LIM)**(XALPHAS)
-!      ZZX(:) = ( -PSSI(:)/PAI(:) ) * (XCCS*PLBDS(:)**XCXS)/PRHODREF(:) * (ZZW(:)**XNUS) * EXP(-ZZW(:))
       ZZX(:) = ( -PSSI(:)/PAI(:) ) * (XLBS*PRST(:)*PLBDS(:)**XBS) * (ZZW(:)**XNUS) * EXP(-ZZW(:))
 !
       ZZW(:) = ( XR0DEPSI+XR1DEPSI*PCJ(:) )*ZZX(:)
@@ -188,11 +188,10 @@ WHERE( GMICRO )
 !
    ZZW(:) = 0.0
    WHERE ( (PRST(:)>XRTMIN(5)) )
-!      ZZW(:) = ( PSSI(:)/(PAI(:))/PRHODREF(:) ) * &
-!           ( X0DEPS*PLBDS(:)**XEX0DEPS + X1DEPS*PCJ(:)*PLBDS(:)**XEX1DEPS )
-      ZZW(:) =10* ( PRHODREF(:) * PRST(:)*PSSI(:)/PAI(:)) *                               &       !Modif Wurtz
-           ( X0DEPS*PLBDS(:)**XEX0DEPS + (X1DEPS*PCJ(:)*(1+(XFVELOS/(2.*PLBDS(:)))**XALPHAS) &
-           **(-XNUS+XEX1DEPS/XALPHAS)*(PLBDS(:))**(XBS+XEX1DEPS))) ! Wurtz
+      ZZW(:) =( PRST(:)*PSSI(:)/PAI(:) ) *                               &
+           ( X0DEPS*PLBDS(:)**XEX0DEPS +                &
+           ( X1DEPS*PCJ(:)*(PLBDS(:))**(XBS+XEX1DEPS) * &
+                    (1+(XFVELOS/(2.*PLBDS(:)))**XALPHAS)**(-XNUS+XEX1DEPS/XALPHAS)))
       ZZW(:) =    ZZW(:)*(0.5+SIGN(0.5,ZZW(:))) - ABS(ZZW(:))*(0.5-SIGN(0.5,ZZW(:)))
    END WHERE
 !
diff --git a/src/MNH/lima_mixed.f90 b/src/MNH/lima_mixed.f90
index f8135b358..4d558eb1d 100644
--- a/src/MNH/lima_mixed.f90
+++ b/src/MNH/lima_mixed.f90
@@ -96,6 +96,7 @@ END MODULE MODI_LIMA_MIXED
 !  P. Wautelet    03/2020: use the new data structures and subroutines for budgets (no more call to budget in this subroutine)
 !  P. Wautelet 28/05/2020: bugfix: correct array start for PSVT and PSVS
 !  P. Wautelet 02/02/2021: budgets: add missing source terms for SV budgets in LIMA
+!  J. Wurtz       03/2022: new snow characteristics
 !-------------------------------------------------------------------------------
 !
 !*       0.    DECLARATIONS
@@ -108,7 +109,7 @@ USE MODD_CST,              ONLY: XP00, XRD, XRV, XMV, XMD, XCPD, XCPV,       &
 USE MODD_NSV
 USE MODD_PARAMETERS,       ONLY: JPHEXT, JPVEXT
 USE MODD_PARAM_LIMA,       ONLY: NMOD_IFN, XRTMIN, XCTMIN, LWARM, LCOLD,     &
-                                 NMOD_CCN, NMOD_IMM, LRAIN, LSNOW, LHAIL
+                                 NMOD_CCN, NMOD_IMM, LRAIN, LSNOW, LHAIL, LSNOW_T
 USE MODD_PARAM_LIMA_WARM,  ONLY: XLBC, XLBEXC, XLBR, XLBEXR
 USE MODD_PARAM_LIMA_COLD,  ONLY: XLBI, XLBEXI, XLBS, XLBEXS, XSCFAC, XLBDAS_MAX, XLBDAS_MIN, XTRANS_MP_GAMMAS
 USE MODD_PARAM_LIMA_MIXED, ONLY: XLBG, XLBEXG, XLBH, XLBEXH
@@ -467,20 +468,19 @@ IF( IMICRO >= 1 ) THEN
       ZLBDAI(:) = ( XLBI*ZCIT(:) / ZRIT(:) )**XLBEXI
    END WHERE
    ZLBDAS(:)  = 1.E10
- !  WHERE (ZRST(:)>XRTMIN(5) )
-  !    ZLBDAS(:) = XLBS*( ZRHODREF(:)*ZRST(:) )**XLBEXS
-   !END WHERE
-!Wurtz 
-      WHERE (ZRST(:)>XRTMIN(5) )
-                        WHERE(ZZT(:)>263.15) 
-                                ZLBDAS(:) = MAX(MIN(XLBDAS_MAX, 10**(14.554-0.0423*ZZT(:))),XLBDAS_MIN)
-                        END WHERE
-
-                       WHERE(ZZT(:)<=263.15) 
-                                ZLBDAS(:) = MAX(MIN(XLBDAS_MAX, 10**(6.226-0.0106*ZZT(:))),XLBDAS_MIN)
-                       END WHERE
+   IF (LSNOW_T) THEN
+      WHERE(ZZT(:)>263.15 .AND. ZRST(:)>XRTMIN(5)) 
+         ZLBDAS(:) = MAX(MIN(XLBDAS_MAX, 10**(14.554-0.0423*ZZT(:))),XLBDAS_MIN)
+      END WHERE
+      WHERE(ZZT(:)<=263.15 .AND. ZRST(:)>XRTMIN(5)) 
+         ZLBDAS(:) = MAX(MIN(XLBDAS_MAX, 10**(6.226-0.0106*ZZT(:))),XLBDAS_MIN)
       END WHERE
-   ZLBDAS(:) = ZLBDAS(:)*XTRANS_MP_GAMMAS
+      ZLBDAS(:) = ZLBDAS(:)*XTRANS_MP_GAMMAS
+   ELSE
+      WHERE (ZRST(:)>XRTMIN(5))
+         ZLBDAS(:) = MAX(MIN(XLBDAS_MAX, XLBS*( ZRHODREF(:)*ZRST(:) )**XLBEXS),XLBDAS_MIN)
+      END WHERE
+   END IF
    ZLBDAG(:)  = 1.E10
    WHERE (ZRGT(:)>XRTMIN(6) )
       ZLBDAG(:) = XLBG*( ZRHODREF(:)*ZRGT(:) )**XLBEXG
diff --git a/src/MNH/lima_mixed_fast_processes.f90 b/src/MNH/lima_mixed_fast_processes.f90
index a5c38e3e0..7cfeffec0 100644
--- a/src/MNH/lima_mixed_fast_processes.f90
+++ b/src/MNH/lima_mixed_fast_processes.f90
@@ -142,6 +142,7 @@ END MODULE MODI_LIMA_MIXED_FAST_PROCESSES
 !!      C. Barthe  * LACy *  jan. 2014    add budgets
 !  P. Wautelet 26/04/2019: replace non-standard FLOAT function by REAL function
 !  P. Wautelet    03/2020: use the new data structures and subroutines for budgets
+!  J. Wurtz       03/2022: new snow characteristics
 !-------------------------------------------------------------------------------
 !
 !*       0.    DECLARATIONS
@@ -296,17 +297,12 @@ IF( IGRIM>0 ) THEN
 !        1.1.4  riming of the small sized aggregates
 !
    WHERE ( GRIM(:) )
-    !  ZZW1(:,1) = MIN( ZRCS(:),                              &
-      !           XCRIMSS * ZZW(:) * ZRCT(:)                & ! RCRIMSS
-         !                           *   ZLBDAS(:)**XEXCRIMSS &
-    !		            * ZRHODREF(:)**(-XCEXVT) )
-!
-      ZZW1(:,1) = MIN( ZRCS(:),                              & ! Wurtz
-             XCRIMSS * ZZW(:) * ZRCT(:)                & ! RCRIMSS
+      ZZW1(:,1) = MIN( ZRCS(:),                              &
+             XCRIMSS * ZZW(:) * ZRCT(:)                &
                                       *  ZRST(:)*(1+(XFVELOS/ZLBDAS(:))**XALPHAS)**(-XNUS+XEXCRIMSS/XALPHAS) &
                                       * ZRHODREF(:)**(-XCEXVT+1) &
-                                      * (ZLBDAS(:)) ** (XEXCRIMSS+XBS) )! Wurtz
-
+                                      * (ZLBDAS(:)) ** (XEXCRIMSS+XBS) )
+!
       ZRCS(:) = ZRCS(:) - ZZW1(:,1)
       ZRSS(:) = ZRSS(:) + ZZW1(:,1)
       ZTHS(:) = ZTHS(:) + ZZW1(:,1)*(ZLSFACT(:)-ZLVFACT(:)) ! f(L_f*(RCRIMSS))
@@ -325,23 +321,14 @@ IF( IGRIM>0 ) THEN
 !
 !
    WHERE ( GRIM(:) .AND. (ZRSS(:)>XRTMIN(5)/PTSTEP) )
-     ! ZZW1(:,2) = MIN( ZRCS(:),                     &
-     !           XCRIMSG * ZRCT(:)                & ! RCRIMSG
-     !                  *  ZLBDAS(:)**XEXCRIMSG  &
-  	  !                 * ZRHODREF(:)**(-XCEXVT) &
-    !            - ZZW1(:,1)              )
-
-      ZZW1(:,2) = MIN( ZRCS(:),                     & ! Modif Wurtz
+      ZZW1(:,2) = MIN( ZRCS(:),                     & 
                    XCRIMSG * ZRCT(:)*  ZRST(:)                & ! RCRIMSG
                            *(1+(XFVELOS/ZLBDAS(:))**XALPHAS)**(-XNUS+XEXCRIMSG/XALPHAS)*ZLBDAS(:)**(XBS+XEXCRIMSG)  &
                            * ZRHODREF(:)**(-XCEXVT+1) &
                            - ZZW1(:,1)              )
 
       ZZW1(:,3) = MIN( ZRSS(:),                         &
-                     !  XSRIMCG * ZLBDAS(:)**XEXSRIMCG   & ! RSRIMCG
-                 !         * (1.0 - ZZW(:) )/(PTSTEP*ZRHODREF(:)))
-                       XSRIMCG * ZRHODREF(:) *ZRST(:)* ZLBDAS(:)**(XEXSRIMCG+XBS)   & ! RSRIMCG ! Wurtz
-                               * (1.0 - ZZW(:) )/(PTSTEP*ZRHODREF(:)))
+                       XSRIMCG * XLBS * ZRST(:) * (1.0 - ZZW(:))/PTSTEP )
 
       ZRCS(:) = ZRCS(:) - ZZW1(:,2)
       ZRSS(:) = ZRSS(:) - ZZW1(:,3)
@@ -493,12 +480,7 @@ IF( IGACC>0 .AND. LRAIN) THEN
 !        1.3.4  raindrop accretion on the small sized aggregates
 !
    WHERE ( GACC(:) )
- !     ZZW1(:,2) = ZCRT(:) *                                           & !! coef of RRACCS
-  !            XFRACCSS*( ZLBDAS(:)**XCXS )*( ZRHODREF(:)**(-XCEXVT-1.) ) &
-   !      *( XLBRACCS1/((ZLBDAS(:)**2)               ) +                  &
-    !        XLBRACCS2/( ZLBDAS(:)    * ZLBDAR(:)    ) +                  &
-     !       XLBRACCS3/(               (ZLBDAR(:)**2)) )/ZLBDAR(:)**3
-      ZZW1(:,2) = ZCRT(:) *                                           & !! coef of RRACCS ! Wurtz
+      ZZW1(:,2) = ZCRT(:) *                                           & !! coef of RRACCS
               XFRACCSS*( ZRST(:)*ZLBDAS(:)**XBS )*( ZRHODREF(:)**(-XCEXVT) ) &
          *( XLBRACCS1/((ZLBDAS(:)**2)               ) +                  &
             XLBRACCS2/( ZLBDAS(:)    * ZLBDAR(:)    ) +                  &
@@ -542,13 +524,7 @@ IF( IGACC>0 .AND. LRAIN) THEN
 !
    WHERE ( GACC(:) .AND. (ZRSS(:)>XRTMIN(5)/PTSTEP) )
       ZZW1(:,2) = MAX( MIN( ZRRS(:),ZZW1(:,2)-ZZW1(:,4) ) , 0. )      ! RRACCSG
-   !   ZZW1(:,3) = MIN( ZRSS(:),XFSACCRG*ZZW(:)*                     & ! RSACCRG
-    !        ( ZLBDAS(:)**(XCXS-XBS) )*( ZRHODREF(:)**(-XCEXVT-1.) ) &
-     !      *( XLBSACCR1/((ZLBDAR(:)**2)               ) +           &
-      !        XLBSACCR2/( ZLBDAR(:)    * ZLBDAS(:)    ) +           &
-       !       XLBSACCR3/(               (ZLBDAS(:)**2)) ) )
-
-      ZZW1(:,3) = MIN( ZRSS(:),XFSACCRG*ZZW(:)*                     & ! RSACCRG ! Wurtz
+      ZZW1(:,3) = MIN( ZRSS(:),XFSACCRG*ZZW(:)*                     & ! RSACCRG
             ( ZRST(:) )*( ZRHODREF(:)**(-XCEXVT) ) &
            *( XLBSACCR1/((ZLBDAR(:)**2)               ) +           &
               XLBSACCR2/( ZLBDAR(:)    * ZLBDAS(:)    ) +           &
@@ -600,14 +576,7 @@ WHERE( (ZRST(:)>XRTMIN(5)) .AND. (ZRSS(:)>XRTMIN(5)/PTSTEP) .AND. (ZZT(:)>XTT) )
 !
 ! compute RSMLT
 !
-  ! ZZW(:)  = MIN( ZRSS(:), XFSCVMG*MAX( 0.0,( -ZZW(:) *             &
-   !                       ( X0DEPS*       ZLBDAS(:)**XEX0DEPS +     &
-    !                        X1DEPS*ZCJ(:)*ZLBDAS(:)**XEX1DEPS ) -   &
-     !                               ( ZZW1(:,1)+ZZW1(:,4) ) *       &
-      !                       ( ZRHODREF(:)*XCL*(XTT-ZZT(:))) ) /    &
-       !                                     ( ZRHODREF(:)*XLMTT ) ) )
-
-   ZZW(:)  = MIN( ZRSS(:), XFSCVMG*MAX( 0.0,( -ZZW(:) *             & ! Wurtz
+   ZZW(:)  = MIN( ZRSS(:), XFSCVMG*MAX( 0.0,( -ZZW(:) *             &
                           ZRHODREF(:) * ZRST(:)*( X0DEPS*            ZLBDAS(:)**XEX0DEPS +     &
                           X1DEPS*ZCJ(:)*(1+(XFVELOS/(2.*ZLBDAS(:)))**XALPHAS) &
                           **(-XNUS+XEX1DEPS/XALPHAS)*(ZLBDAS(:))**(XEX1DEPS+XBS))-    &
@@ -775,22 +744,13 @@ IF( IGDRY>0 ) THEN
    ZZW(:) = UNPACK( VECTOR=ZVEC3(:),MASK=GDRY,FIELD=0.0 )
 !
    WHERE( GDRY(:) )  
-  !    ZZW1(:,3) = MIN( ZRSS(:),XFSDRYG*ZZW(:)                         & ! RSDRYG
-   !                                   * EXP( XCOLEXSG*(ZZT(:)-XTT) )  &
-    !                *( ZLBDAS(:)**(XCXS-XBS) )*( ZLBDAG(:)**XCXG )    &
-     !               *( ZRHODREF(:)**(-XCEXVT-1.) )                    &
-      !                   *( XLBSDRYG1/( ZLBDAG(:)**2              ) + &
-       !                     XLBSDRYG2/( ZLBDAG(:)   * ZLBDAS(:)   ) + &
-        !                    XLBSDRYG3/(               ZLBDAS(:)**2) ) )
-
-
       ZZW1(:,3) = MIN( ZRSS(:),XFSDRYG*ZZW(:)                         & ! RSDRYG
                                       * EXP( XCOLEXSG*(ZZT(:)-XTT) )  &
                     *( ZRST(:)) )*( ZLBDAG(:)**XCXG )    &
                     *( ZRHODREF(:)**(-XCEXVT) )                    &
                          *( XLBSDRYG1/( ZLBDAG(:)**2              ) + &
                             XLBSDRYG2/( ZLBDAG(:)   * ZLBDAS(:)   ) + &
-                            XLBSDRYG3/(               ZLBDAS(:)**2) ) 
+                            XLBSDRYG3/(               ZLBDAS(:)**2) )
    END WHERE
    DEALLOCATE(IVEC2)
    DEALLOCATE(IVEC1)
@@ -1229,14 +1189,6 @@ IF( IHAIL>0 ) THEN
       ZZW(:) = UNPACK( VECTOR=ZVEC3(:),MASK=GWET,FIELD=0.0 )
 !
       WHERE( GWET(:) )
-!        ZZW1(:,3) = MIN( ZRSS(:),XFSWETH*ZZW(:)                       & ! RSWETH
- !                     *( ZLBDAS(:)**(XCXS-XBS) )*( ZLBDAH(:)**XCXH )  &
-  !     	                 *( ZRHODREF(:)**(-XCEXVT-1.) )               &
-   !                      *( XLBSWETH1/( ZLBDAH(:)**2              ) + &
-     !                       XLBSWETH2/( ZLBDAH(:)   * ZLBDAS(:)   ) + &
-    !                        XLBSWETH3/(               ZLBDAS(:)**2) ) )
-
-! Wurtz 
         ZZW1(:,3) = MIN( ZRSS(:),XFSWETH*ZZW(:)                       & ! RSWETH
                       *( ZRST(:))*( ZLBDAH(:)**XCXH )  &
                         *( ZRHODREF(:)**(-XCEXVT) )               &
diff --git a/src/MNH/lima_rain_accr_snow.f90 b/src/MNH/lima_rain_accr_snow.f90
index 495117c2b..7c5c9bee5 100644
--- a/src/MNH/lima_rain_accr_snow.f90
+++ b/src/MNH/lima_rain_accr_snow.f90
@@ -60,6 +60,7 @@ END MODULE MODI_LIMA_RAIN_ACCR_SNOW
 !!    -------------
 !!      Original             15/03/2018
 !  P. Wautelet 26/04/2019: replace non-standard FLOAT function by REAL function
+!  J. Wurtz       03/2022: new snow characteristics
 !
 !-------------------------------------------------------------------------------
 !
@@ -212,52 +213,28 @@ WHERE( GACC )
 !      
 ! BVIE manque PCRT ???????????????????????????????????
 !      ZZW4(:) =                                            & !! coef of RRACCS and RRACCS
-!   ZZW4(:) = PCRT(:)                                     & !! coef of RRACCS and RRACCS
-!         *  XFRACCSS *( PLBDS(:)**XCXS )*( PRHODREF(:)**(-XCEXVT-1.) ) &
-!         *( XLBRACCS1/( PLBDS(:)**2 )                +                     &
-!            XLBRACCS2/( PLBDS(:) * PLBDR(:)    ) +                     &
-!            XLBRACCS3/(                PLBDR(:)**2 ) ) / PLBDR(:)**3
-   ZZW4(:) = PCRT(:)                                     & !! coef of RRACCS and RRACCS ! Wurtz
+   ZZW4(:) = PCRT(:)                                                       & !! coef of RRACCS and RRACCS
          *  XFRACCSS *( PRST(:)*PLBDS(:)**XBS )*( PRHODREF(:)**(-XCEXVT) ) &
-         *( XLBRACCS1/( PLBDS(:)**2 )                +                     &
-            XLBRACCS2/( PLBDS(:) * PLBDR(:)    ) +                     &
-            XLBRACCS3/(                PLBDR(:)**2 ) ) / PLBDR(:)**3
+         *( XLBRACCS1/( PLBDS(:)**2               ) +                      &
+            XLBRACCS2/( PLBDS(:)    * PLBDR(:)    ) +                      &
+            XLBRACCS3/(               PLBDR(:)**2 ) ) / PLBDR(:)**3
 
-!      ZRRS(:) = ZRRS(:) - ZZW1(:,4)
-!      ZRSS(:) = ZRSS(:) + ZZW1(:,4)
-!      ZTHS(:) = ZTHS(:) + ZZW1(:,4)*(ZLSFACT(:)-ZLVFACT(:)) ! f(L_f*(RRACCSS))
-!
-!      ZCRS(:) = MAX( ZCRS(:)-ZZW1(:,4)*(ZCRT(:)/ZRRT(:)),0.0 ) ! Lambda_r**3 
 !
 !        1.3.6  raindrop accretion-conversion of the large sized aggregates
 !               into graupeln
 !
-!   ZZW5(:) = XFSACCRG*ZZW3(:) *                             & ! RSACCRG
-!            ( PLBDS(:)**(XCXS-XBS) )*( PRHODREF(:)**(-XCEXVT-1.) ) &
-!           *( XLBSACCR1/((PLBDR(:)**2)               ) +           &
-!              XLBSACCR2/( PLBDR(:)    * PLBDS(:) ) +           &
-!              XLBSACCR3/(                  (PLBDS(:)**2)) )
-   ZZW5(:) = XFSACCRG*ZZW3(:) *                             & ! RSACCRG ! Wurtz
-            ( PRST(:) )*( PRHODREF(:)**(-XCEXVT) ) &
-           *( XLBSACCR1/((PLBDR(:)**2)               ) +           &
-              XLBSACCR2/( PLBDR(:)    * PLBDS(:) ) +           &
-              XLBSACCR3/(                  (PLBDS(:)**2)) )
-      !
-!      P_RR_ACC(:) = - ZZW4(:) * ZZW1(:)           ! RRACCSS
-!      P_CR_ACC(:) = P_RR_ACC(:) * PCRT(:)/PRRT(:) ! Lambda_r**3 
-!      P_RS_ACC(:) = - P_RR_ACC(:) 
-      !
-!      P_RR_ACC(:) = P_RR_ACC(:) - ( ZZW2(:)-P_RS_ACC(:) )
-!      P_CR_ACC(:) = P_CR_ACC(:) - ( ZZW2(:)-P_RS_ACC(:) ) * PCRT(:)/PRRT(:) ! Lambda_r**3
-!      P_RS_ACC(:) = P_RS_ACC(:) - ZZW5(:)
-!      P_RG_ACC(:) = ( ZZW2(:)-P_RS_ACC(:) ) + ZZW5(:)
-      !
+   ZZW5(:) = XFSACCRG*ZZW3(:) *                         & ! RSACCRG
+            ( PRST(:) )*( PRHODREF(:)**(-XCEXVT)    )   &
+           *( XLBSACCR1/( PLBDR(:)**2               ) + &
+              XLBSACCR2/( PLBDR(:)    * PLBDS(:)    ) + &
+              XLBSACCR3/(               PLBDS(:)**2 ) )
+!
    P_RR_ACC(:) = - ZZW4(:) *  ZZW2(:)
    P_CR_ACC(:) = P_RR_ACC(:) * PCRT(:)/PRRT(:)
    P_RS_ACC(:) = ZZW4(:) *  ZZW1(:) - ZZW5(:)
    P_RG_ACC(:) = ZZW4(:) * ( ZZW2(:) - ZZW1(:) ) + ZZW5(:)
    P_TH_ACC(:) = - P_RR_ACC(:) * (PLSFACT(:)-PLVFACT(:))
-   !
+!
 END WHERE
 !
 !
diff --git a/src/MNH/lima_sedimentation.f90 b/src/MNH/lima_sedimentation.f90
index 766c7c119..8d48b776d 100644
--- a/src/MNH/lima_sedimentation.f90
+++ b/src/MNH/lima_sedimentation.f90
@@ -66,6 +66,7 @@ END MODULE MODI_LIMA_SEDIMENTATION
 !  P. Wautelet 26/04/2019: replace non-standard FLOAT function by REAL function
 !  P. Wautelet 28/05/2019: move COUNTJV function to tools.f90
 !  B. Vie         03/2020: disable temperature change of droplets by air temperature
+!  J. Wurtz       03/2022: new snow characteristics
 !-------------------------------------------------------------------------------
 !
 !*       0.    DECLARATIONS
@@ -75,7 +76,7 @@ USE MODD_CST,              ONLY: XRHOLW, XCL, XCI
 USE MODD_PARAMETERS,       ONLY: JPHEXT, JPVEXT
 USE MODD_PARAM_LIMA,       ONLY: XCEXVT, XRTMIN, XCTMIN, NSPLITSED,           &
                                  XLB, XLBEX, XD, XFSEDR, XFSEDC,              &
-                                 XALPHAC, XNUC, XALPHAS, XNUS
+                                 XALPHAC, XNUC, XALPHAS, XNUS, LSNOW_T
 USE MODD_PARAM_LIMA_COLD,  ONLY: XLBEXI, XLBI, XDI, XLBDAS_MAX, XBS, XEXSEDS, &
                                  XLBDAS_MIN, XTRANS_MP_GAMMAS, XFVELOS
 
@@ -184,7 +185,7 @@ DO JN = 1 ,  NSPLITSED(KID)
          IF (KMOMENTS==2) ZCS(JL) = PCS(I1(JL),I2(JL),I3(JL))
       END DO
 !
-      IF (KID == 5) THEN
+      IF (KID == 5 .AND. LSNOW_T) THEN
          ZLBDA(:) = 1.E10
          WHERE(ZT(:)>263.15 .AND. ZRS(:)>XRTMIN(5))
             ZLBDA(:) = MAX(MIN(XLBDAS_MAX, 10**(14.554-0.0423*ZT(:))),XLBDAS_MIN)
diff --git a/src/MNH/lima_snow_deposition.f90 b/src/MNH/lima_snow_deposition.f90
index a883661f7..1f7e1dde4 100644
--- a/src/MNH/lima_snow_deposition.f90
+++ b/src/MNH/lima_snow_deposition.f90
@@ -63,6 +63,7 @@ SUBROUTINE LIMA_SNOW_DEPOSITION (LDCOMPUTE,                                &
 !!    -------------
 !!      Original             15/03/2018
 !!
+!  J. Wurtz       03/2022: new snow characteristics
 !-------------------------------------------------------------------------------
 !
 !*       0.    DECLARATIONS
@@ -132,7 +133,6 @@ WHERE( GMICRO )
    WHERE ( PLBDS(:)<XLBDASCNVI_MAX .AND. (PRST(:)>XRTMIN(5)) &
                                    .AND. (PSSI(:)<0.0)       )
       ZZW(:) = (PLBDS(:)*XDSCNVI_LIM)**(XALPHAS)
-!      ZZX(:) = ( -PSSI(:)/PAI(:) ) * (XCCS*PLBDS(:)**XCXS) * (ZZW(:)**XNUS) * EXP(-ZZW(:))
       ZZX(:) = ( -PSSI(:)/PAI(:) ) * (XLBS*PRST(:)*PLBDS(:)**XBS) * (ZZW(:)**XNUS) * EXP(-ZZW(:))
 !
       ZZW(:) = ( XR0DEPSI+XR1DEPSI*PCJ(:) )*ZZX(:)
@@ -150,11 +150,10 @@ WHERE( GMICRO )
 !
    ZZW(:) = 0.0
    WHERE ( (PRST(:)>XRTMIN(5)) )
-!      ZZW(:) = ( PSSI(:)/(PAI(:)*PRHODREF(:)) ) * &
-!           ( X0DEPS*PLBDS(:)**XEX0DEPS + X1DEPS*PCJ(:)*PLBDS(:)**XEX1DEPS )
-      ZZW(:) = ( PRHODREF(:) * PRST(:)*PSSI(:)/(PAI(:)) ) * & ! Wurtz
-           ( X0DEPS*PLBDS(:)**XEX0DEPS + (X1DEPS*PCJ(:)* &
-           (1+(XFVELOS/(2.*PLBDS))**XALPHAS)**(-XNUS+XEX1DEPS/XALPHAS)*(PLBDS(:))**(XEX1DEPS+XBS)) ) ! Wurtz
+      ZZW(:) = ( PRST(:)*PSSI(:)/(PAI(:)) ) *           &
+           ( X0DEPS*PLBDS(:)**XEX0DEPS +                &
+           ( X1DEPS*PCJ(:)*(PLBDS(:))**(XBS+XEX1DEPS) * &
+                (1+(XFVELOS/(2.*PLBDS))**XALPHAS)**(-XNUS+XEX1DEPS/XALPHAS)) )
       ZZW(:) =    ZZW(:)*(0.5+SIGN(0.5,ZZW(:))) - ABS(ZZW(:))*(0.5-SIGN(0.5,ZZW(:)))
    END WHERE
 !
diff --git a/src/MNH/lima_tendencies.f90 b/src/MNH/lima_tendencies.f90
index bc49940dc..96e68b530 100644
--- a/src/MNH/lima_tendencies.f90
+++ b/src/MNH/lima_tendencies.f90
@@ -231,6 +231,7 @@ SUBROUTINE LIMA_TENDENCIES (PTSTEP, LDCOMPUTE,
 !!    -------------
 !!      Original             15/03/2018
 !!
+!  J. Wurtz       03/2022: new snow characteristics
 !-------------------------------------------------------------------------------
 !
 !*       0.    DECLARATIONS
@@ -238,8 +239,8 @@ SUBROUTINE LIMA_TENDENCIES (PTSTEP, LDCOMPUTE,
 !
 USE MODD_CST,              ONLY : XP00, XRD, XRV, XMD, XMV, XCPD, XCPV, XCL, XCI, XLVTT, XLSTT, XTT, &
                                   XALPW, XBETAW, XGAMW, XALPI, XBETAI, XGAMI
-USE MODD_PARAM_LIMA,       ONLY : XRTMIN, XCTMIN,                                                    &
-                                  LCOLD, LNUCL, LSNOW, LHAIL, LWARM, LACTI, LRAIN,XNUS
+USE MODD_PARAM_LIMA,       ONLY : XRTMIN, XCTMIN, XNUS,                                    &
+                                  LCOLD, LNUCL, LSNOW, LHAIL, LWARM, LACTI, LRAIN, LSNOW_T
 USE MODD_PARAM_LIMA_WARM,  ONLY : XLBC, XLBEXC, XLBR, XLBEXR
 USE MODD_PARAM_LIMA_MIXED, ONLY : XLBG, XLBEXG, XLBH, XLBEXH, XLBDAG_MAX
 USE MODD_PARAM_LIMA_COLD,  ONLY : XSCFAC, XLBI, XLBEXI, XLBS, XLBEXS, XLBDAS_MAX, XTRANS_MP_GAMMAS,  &
@@ -503,18 +504,21 @@ WHERE (PRIT(:)>XRTMIN(4) .AND. PCIT(:)>XCTMIN(4) .AND. LDCOMPUTE(:))
    ZLBDI(:) = ( XLBI*PCIT(:) / PRIT(:) )**XLBEXI
 END WHERE
 ZLBDS(:)  = 1.E10
-!WHERE (PRST(:)>XRTMIN(5) .AND. LDCOMPUTE(:) )
-!   ZLBDS(:) = XLBS*( PRHODREF(:)*PRST(:) )**XLBEXS
-!END WHERE
-WHERE (PRST(:)>XRTMIN(5) .AND. LDCOMPUTE(:) )
-   WHERE(ZT(:)>263.15)
-      ZLBDS(:) = MAX(MIN(XLBDAS_MAX, 10**(14.554-0.0423*ZT(:))),XLBDAS_MIN)
+IF (LSNOW_T) THEN
+   WHERE (PRST(:)>XRTMIN(5) .AND. LDCOMPUTE(:) )
+      WHERE(ZT(:)>263.15)
+         ZLBDS(:) = MAX(MIN(XLBDAS_MAX, 10**(14.554-0.0423*ZT(:))),XLBDAS_MIN)
+      END WHERE
+      WHERE(ZT(:)<=263.15)
+         ZLBDS(:) = MAX(MIN(XLBDAS_MAX, 10**(6.226-0.0106*ZT(:))),XLBDAS_MIN)
+      END WHERE
    END WHERE
-   WHERE(ZT(:)<=263.15)
-      ZLBDS(:) = MAX(MIN(XLBDAS_MAX, 10**(6.226-0.0106*ZT(:))),XLBDAS_MIN)
+   ZLBDS(:) =  ZLBDS(:) * XTRANS_MP_GAMMAS
+ELSE
+   WHERE (PRST(:)>XRTMIN(5) .AND. LDCOMPUTE(:) )
+      ZLBDS(:) = XLBS*( PRHODREF(:)*PRST(:) )**XLBEXS
    END WHERE
-END WHERE                               ! Wurtz
-ZLBDS(:) =  ZLBDS(:) * XTRANS_MP_GAMMAS ! Wurtz
+END IF
 ZLBDG(:)  = 1.E10
 WHERE (PRGT(:)>XRTMIN(6) .AND. LDCOMPUTE(:) )
    ZLBDG(:) = XLBG*( PRHODREF(:)*PRGT(:) )**XLBEXG
diff --git a/src/MNH/modd_param_lima.f90 b/src/MNH/modd_param_lima.f90
index 66156a056..c32d7a606 100644
--- a/src/MNH/modd_param_lima.f90
+++ b/src/MNH/modd_param_lima.f90
@@ -80,6 +80,7 @@ REAL, DIMENSION(:),    SAVE, ALLOCATABLE :: XFRAC_REF       ! AP compostion in P
 !
 ! 1.3 Ice characteristics
 !
+LOGICAL, SAVE :: LSNOW_T                     ! TRUE to enable snow param. after Wurtz 2021
 CHARACTER(LEN=4), SAVE :: CPRISTINE_ICE_LIMA ! Pristine type PLAT, COLU or BURO
 CHARACTER(LEN=4), SAVE :: CHEVRIMED_ICE_LIMA ! Heavily rimed type GRAU or HAIL
 REAL,SAVE              :: XALPHAI,XNUI,    & ! Pristine ice   distribution parameters
diff --git a/src/MNH/modd_param_lima_cold.f90 b/src/MNH/modd_param_lima_cold.f90
index 37636b7ac..65cb7f9a6 100644
--- a/src/MNH/modd_param_lima_cold.f90
+++ b/src/MNH/modd_param_lima_cold.f90
@@ -19,6 +19,7 @@
 !!    MODIFICATIONS
 !!    -------------
 !!      Original             ??/??/13 
+!  J. Wurtz       03/2022: new snow characteristics
 !!
 !-------------------------------------------------------------------------------
 USE MODD_PARAMETERS, ONLY: JPSVNAMELGTMAX
@@ -52,8 +53,8 @@ REAL,SAVE ::                           XF0IS,XF1IS    ! (large Di vent. coef.)
 REAL,SAVE :: XAS,XBS,XCS,XDS,XCCS,XCXS,XF0S,XF1S,XC1S ! Snow/agg.      charact.
 !
 REAL,SAVE :: XLBDAS_MIN, XLBDAS_MAX   ! Max values allowed for the shape parameter of snow
-REAL,SAVE :: XFVELOS                  ! Wurtz
-REAL,SAVE :: XTRANS_MP_GAMMAS         ! Wurtz
+REAL,SAVE :: XFVELOS                  ! Wurtz - snow fall speed parameterizaed after Thompson 2008
+REAL,SAVE :: XTRANS_MP_GAMMAS         ! Wurtz - change between lambda value for MP and gen. gamma
 !
 CHARACTER(LEN=JPSVNAMELGTMAX),DIMENSION(5),PARAMETER &
                               :: CLIMA_COLD_NAMES=(/'CICE    ','CIFNFREE','CIFNNUCL', &
diff --git a/src/MNH/modd_param_lima_mixed.f90 b/src/MNH/modd_param_lima_mixed.f90
index 2a2b75491..57ea4d1a5 100644
--- a/src/MNH/modd_param_lima_mixed.f90
+++ b/src/MNH/modd_param_lima_mixed.f90
@@ -14,6 +14,7 @@
 !!    MODIFICATIONS
 !!    -------------
 !!      Original             ??/??/13 
+!  J. Wurtz       03/2022: new snow characteristics
 !!
 !-------------------------------------------------------------------------------
 !
diff --git a/src/MNH/modn_param_lima.f90 b/src/MNH/modn_param_lima.f90
index f86b1add0..e8e56a409 100644
--- a/src/MNH/modn_param_lima.f90
+++ b/src/MNH/modn_param_lima.f90
@@ -20,7 +20,7 @@ IMPLICIT NONE
 NAMELIST/NAM_PARAM_LIMA/LCOLD, LNUCL, LSEDI, LSNOW, LHAIL, LHHONI, LMEYERS,&
                         NMOD_IFN, XIFN_CONC, LIFN_HOM,                     &
                         CIFN_SPECIES, CINT_MIXING, NMOD_IMM, NIND_SPECIE,  &
-                        CPRISTINE_ICE_LIMA, CHEVRIMED_ICE_LIMA,            &
+                        LSNOW_T, CPRISTINE_ICE_LIMA, CHEVRIMED_ICE_LIMA,   &
                         XALPHAI, XNUI, XALPHAS, XNUS, XALPHAG, XNUG,       &
                         XFACTNUC_DEP, XFACTNUC_CON, NPHILLIPS,             &
 !
diff --git a/src/MNH/profilern.f90 b/src/MNH/profilern.f90
index 433a5c6ef..e2025c578 100644
--- a/src/MNH/profilern.f90
+++ b/src/MNH/profilern.f90
@@ -112,10 +112,11 @@ USE MODE_FGAU,             ONLY : GAULAG
 USE MODE_FSCATTER,         ONLY: QEPSW,QEPSI,BHMIE,MOMG,MG
 USE MODD_PARAM_LIMA,       ONLY: XALPHAR_L=>XALPHAR,XNUR_L=>XNUR,XALPHAS_L=>XALPHAS,XNUS_L=>XNUS,&
                                  XALPHAG_L=>XALPHAG,XNUG_L=>XNUG, XALPHAI_L=>XALPHAI,XNUI_L=>XNUI,&
-                                 XRTMIN_L=>XRTMIN,XALPHAC_L=>XALPHAC,XNUC_L=>XNUC
+                                 XRTMIN_L=>XRTMIN,XALPHAC_L=>XALPHAC,XNUC_L=>XNUC, LSNOW_T_L=>LSNOW_T
 USE MODD_PARAM_LIMA_COLD,  ONLY: XDI_L=>XDI,XLBEXI_L=>XLBEXI,XLBI_L=>XLBI,XAI_L=>XAI,XBI_L=>XBI,XC_I_L=>XC_I,&
                                  XLBEXS_L=>XLBEXS,XLBS_L=>XLBS,XCCS_L=>XCCS,&
-                                 XAS_L=>XAS,XBS_L=>XBS,XCXS_L=>XCXS
+                                 XAS_L=>XAS,XBS_L=>XBS,XCXS_L=>XCXS,        &
+                                 XLBDAS_MIN,XLBDAS_MAX
 USE MODD_PARAM_LIMA_MIXED, ONLY: XDG_L=>XDG,XLBEXG_L=>XLBEXG,XLBG_L=>XLBG,XCCG_L=>XCCG,&
                                  XAG_L=>XAG,XBG_L=>XBG,XCXG_L=>XCXG,XCG_L=>XCG
 USE MODD_PARAM_LIMA_WARM,  ONLY: XLBEXR_L=>XLBEXR,XLBR_L=>XLBR,XBR_L=>XBR,XAR_L=>XAR,&
@@ -242,6 +243,7 @@ REAL                           :: ZAETOT,ZAETMP,ZREFLOC,ZQSCA,ZQBACK,ZQEXT ! tem
 REAL,DIMENSION(:),ALLOCATABLE  :: ZAELOC,ZZMZ ! temporary arrays
 INTEGER                        :: JPTS_GAULAG=9 ! number of points for Gauss-Laguerre quadrature
 REAL                           :: ZLBDA   ! slope distribution parameter
+REAL                           :: ZN   ! number cocentration
 REAL                           :: ZFRAC_ICE  ! ice water fraction
 REAL                           :: ZDELTA_EQUIV ! mass-equivalent Gauss-Laguerre point
 REAL                           :: ZFW ! liquid fraction
@@ -724,7 +726,17 @@ IF (GSTORE) THEN
                     ZLB=( ZA*ZCC*MOMG(ZALPHA,ZNU,ZB) )**(-ZLBEX)
                   ENDIF
             END SELECT
-            ZLBDA=ZLB*(ZRHODREFZ(JK)*ZRZ(JK,JLOOP))**ZLBEX
+            IF (JLOOP.EQ.5 .AND. CCLOUD=='LIMA' .AND. LSNOW_T_L) THEN
+               IF (ZTEMPZ(JK)>-10.) THEN
+                  ZLBDA = MAX(MIN(XLBDAS_MAX, 10**(14.554-0.0423*(ZTEMPZ(JK)+273.15))),XLBDAS_MIN)
+               ELSE
+                  ZLBDA = MAX(MIN(XLBDAS_MAX, 10**(6.226-0.0106*(ZTEMPZ(JK)+273.15))),XLBDAS_MIN)
+               END IF
+               ZN=XLBS_L*ZRHODREFZ(JK)*ZRZ(JK,JLOOP)*ZLBDA**ZB
+            ELSE
+               ZLBDA=ZLB*(ZRHODREFZ(JK)*ZRZ(JK,JLOOP))**ZLBEX
+               ZN=ZCC*ZLBDA**ZCX
+            END IF
             ZREFLOC=0.
             ZAETMP=0.
             DO JJ=1,JPTS_GAULAG ! Gauss-Laguerre quadrature
@@ -754,8 +766,8 @@ IF (GSTORE) THEN
                ZREFLOC=ZREFLOC+ZQBACK*ZX(JJ)**(ZNU-1)*ZDELTA_EQUIV**2*ZW(JJ)
                ZAETMP =ZAETMP +ZQEXT *ZX(JJ)**(ZNU-1)*ZDELTA_EQUIV**2*ZW(JJ)
              END DO
-             ZREFLOC=ZREFLOC*(XLAM_CRAD/XPI)**4*ZCC*ZLBDA**ZCX/(4.*GAMMA(ZNU)*.93)
-             ZAETMP=ZAETMP  *           XPI    *ZCC*ZLBDA**ZCX/(4.*GAMMA(ZNU))
+             ZREFLOC=ZREFLOC*(XLAM_CRAD/XPI)**4*ZN/(4.*GAMMA(ZNU)*.93)
+             ZAETMP=ZAETMP  *           XPI    *ZN/(4.*GAMMA(ZNU))
              TPROFILER%CRARE(IN,JK,I)=TPROFILER%CRARE(IN,JK,I)+ZREFLOC
              ZAELOC(JK)=ZAELOC(JK)+ZAETMP
            END IF
diff --git a/src/MNH/radar_rain_ice.f90 b/src/MNH/radar_rain_ice.f90
index eddac2294..b42848f56 100644
--- a/src/MNH/radar_rain_ice.f90
+++ b/src/MNH/radar_rain_ice.f90
@@ -114,7 +114,8 @@ USE MODD_PARAM_LIMA_WARM, ONLY: XLBEXR_L=>XLBEXR,XLBR_L=>XLBR,XBR_L=>XBR,XAR_L=>
                                 XBC_L=>XBC,XAC_L=>XAC,XCR_L=>XCR,XDR_L=>XDR
 USE MODD_PARAM_LIMA_COLD, ONLY: XDI_L=>XDI,XLBEXI_L=>XLBEXI,XLBI_L=>XLBI,XAI_L=>XAI,XBI_L=>XBI,XC_I_L=>XC_I,&
                                 XLBEXS_L=>XLBEXS,XLBS_L=>XLBS,XCCS_L=>XCCS,&
-                                XAS_L=>XAS,XBS_L=>XBS,XCXS_L=>XCXS,XCS_L=>XCS,XDS_L=>XDS
+                                XAS_L=>XAS,XBS_L=>XBS,XCXS_L=>XCXS,XCS_L=>XCS,XDS_L=>XDS,&
+                                XLBDAS_MIN,XLBDAS_MAX
 
 USE MODD_PARAM_LIMA_MIXED, ONLY:XDG_L=>XDG,XLBEXG_L=>XLBEXG,XLBG_L=>XLBG,XCCG_L=>XCCG,&
                                 XAG_L=>XAG,XBG_L=>XBG,XCXG_L=>XCXG,XCG_L=>XCG,&
@@ -123,7 +124,7 @@ USE MODD_PARAM_LIMA_MIXED, ONLY:XDG_L=>XDG,XLBEXG_L=>XLBEXG,XLBG_L=>XLBG,XCCG_L=
 
 USE MODD_PARAM_LIMA, ONLY: XALPHAR_L=>XALPHAR,XNUR_L=>XNUR,XALPHAS_L=>XALPHAS,XNUS_L=>XNUS,&
                            XALPHAG_L=>XALPHAG,XNUG_L=>XNUG, XALPHAI_L=>XALPHAI,XNUI_L=>XNUI,&
-                           XRTMIN_L=>XRTMIN,XALPHAC_L=>XALPHAC,XNUC_L=>XNUC                      
+                           XRTMIN_L=>XRTMIN,XALPHAC_L=>XALPHAC,XNUC_L=>XNUC, LSNOW_T_L=>LSNOW_T                      
 USE MODD_PARAMETERS
 USE MODD_PARAM_n, ONLY : CCLOUD
 USE MODD_LUNIT
@@ -169,6 +170,8 @@ REAL :: ZRHO00                        ! Surface reference air density
 LOGICAL, DIMENSION(SIZE(PTEMP,1),SIZE(PTEMP,2),SIZE(PTEMP,3)) :: GRAIN
 REAL,    DIMENSION(SIZE(PTEMP,1),SIZE(PTEMP,2),SIZE(PTEMP,3)) :: ZLBDA 
                                       ! slope distribution parameter
+REAL,    DIMENSION(SIZE(PTEMP,1),SIZE(PTEMP,2),SIZE(PTEMP,3)) :: ZN 
+                                      ! number concentration
 REAL,    DIMENSION(SIZE(PTEMP,1),SIZE(PTEMP,2),SIZE(PTEMP,3)) :: ZW
 REAL,    DIMENSION(SIZE(PTEMP,1),SIZE(PTEMP,2),SIZE(PTEMP,3)) :: ZREFL_MELT_CONV
 INTEGER                                                       :: JLBDA
@@ -325,7 +328,24 @@ END IF
 !               ---------------
 !
 IF (SIZE(PRT,4) >= 5) THEN
-  IF (CCLOUD=='LIMA') THEN
+  IF (CCLOUD=='LIMA' .AND. LSNOW_T_L) THEN
+    ZDMELT_FACT = ( (6.0*XAS_L)/(XPI*XRHOLW) )**(2.0)
+    ZEXP = 2.0*XBS_L
+    WHERE(PTEMP(:,:,:)>-10. .AND. PRT(:,:,:,5).GT.XRTMIN_L(5))
+       ZLBDA(:,:,:) = MAX(MIN(XLBDAS_MAX, 10**(14.554-0.0423*(PTEMP(:,:,:)+273.15))),XLBDAS_MIN)
+    END WHERE
+    WHERE(PTEMP(:,:,:)<=-10 .AND. PRT(:,:,:,5).GT.XRTMIN_L(5))
+       ZLBDA(:,:,:) = MAX(MIN(XLBDAS_MAX, 10**(6.226-0.0106*(PTEMP(:,:,:)+273.15))),XLBDAS_MIN)
+    END WHERE
+    ZN(:,:,:)=XLBS_L*PRHODREF(:,:,:)*PRT(:,:,:,5)*ZLBDA(:,:,:)**XBS_L
+    WHERE( PRT(:,:,:,5).GT.XRTMIN_L(5) )
+      ZW(:,:,:) = ZEQICE*ZDMELT_FACT                                             &
+                  *1.E18*ZN(:,:,:)*(ZLBDA(:,:,:)**(-ZEXP))*MOMG(XALPHAS_L,XNUS_L,ZEXP)
+      PVDOP(:,:,:) = PVDOP(:,:,:)+ZEQICE*ZDMELT_FACT*MOMG(XALPHAS_L,XNUS_L,ZEXP+XDS_L) &
+                     *1.E18*ZN(:,:,:)*XCS_L*(ZLBDA(:,:,:)**(-ZEXP-XDS_L))
+      PRARE(:,:,:) = PRARE(:,:,:) + ZW(:,:,:)
+    END WHERE
+  ELSEIF (CCLOUD=='LIMA') THEN
     ZDMELT_FACT = ( (6.0*XAS_L)/(XPI*XRHOLW) )**(2.0)
     ZEXP = 2.0*XBS_L
     WHERE( PRT(:,:,:,5).GT.XRTMIN_L(5) )
diff --git a/src/MNH/radar_scattering.f90 b/src/MNH/radar_scattering.f90
index dc3ddc1e5..413bbbb92 100644
--- a/src/MNH/radar_scattering.f90
+++ b/src/MNH/radar_scattering.f90
@@ -119,12 +119,13 @@ USE MODD_RAIN_ICE_DESCR, ONLY: XALPHAR_I=>XALPHAR,XNUR_I=>XNUR,XDR_I=>XDR,XLBEXR
 !!LIMA         
 USE MODD_PARAM_LIMA_WARM, ONLY: XDR_L=>XDR,XLBEXR_L=>XLBEXR,XLBR_L=>XLBR,XBR_L=>XBR,XCR_L=>XCR
 USE MODD_PARAM_LIMA_COLD, ONLY: XDI_L=>XDI,XLBEXI_L=>XLBEXI,XLBI_L=>XLBI,XAI_L=>XAI,XBI_L=>XBI,XC_I_L=>XC_I,&
-                                XDS_L=>XDS,XLBEXS_L=>XLBEXS,XLBS_L=>XLBS,XCCS_L=>XCCS,XAS_L=>XAS,XBS_L=>XBS,XCXS_L=>XCXS,XCS_L=>XCS
+                                XDS_L=>XDS,XLBEXS_L=>XLBEXS,XLBS_L=>XLBS,XCCS_L=>XCCS,XAS_L=>XAS,XBS_L=>XBS,&
+                                XCXS_L=>XCXS,XCS_L=>XCS,XLBDAS_MIN,XLBDAS_MAX
 
 USE MODD_PARAM_LIMA_MIXED, ONLY:XDG_L=>XDG,XLBEXG_L=>XLBEXG,XLBG_L=>XLBG,XCCG_L=>XCCG,XAG_L=>XAG,XBG_L=>XBG,XCXG_L=>XCXG,XCG_L=>XCG
 USE MODD_PARAM_LIMA, ONLY: XALPHAR_L=>XALPHAR,XNUR_L=>XNUR,XALPHAS_L=>XALPHAS,XNUS_L=>XNUS,&
                            XALPHAG_L=>XALPHAG,XNUG_L=>XNUG, XALPHAI_L=>XALPHAI,XNUI_L=>XNUI,&
-                           XRTMIN_L=>XRTMIN
+                           XRTMIN_L=>XRTMIN, LSNOW_T_L=>LSNOW_T
 !!LIMA
 USE MODD_RADAR, ONLY:XLAM_RAD,XSTEP_RAD,NBELEV,NDIFF,LATT,NPTS_GAULAG,LQUAD,XVALGROUND,NDGS, &
      LFALL,LWBSCS,LWREFL,XREFLVDOPMIN,XREFLMIN,LSNRT,XSNRMIN
@@ -192,6 +193,7 @@ REAL :: ZDMELT_FACT ! factor used to compute the equivalent melted diameter
 REAL :: ZEQICE=0.224! factor used to convert the ice crystals reflectivity into an equivalent  liquid water reflectivity (from Smith, JCAM 84)
 REAL :: ZEXP        ! anciliary parameter
 REAL :: ZLBDA   ! slope distribution parameter
+REAL :: ZN      ! Number concentration
 REAL :: ZFRAC_ICE,ZD,ZDE ! auxiliary variables
 REAL :: ZQSCA
 REAL,DIMENSION(2) :: ZQEXT
@@ -1282,28 +1284,37 @@ DO JI=1,INBRAD
                     ZDMELT_FACT=6.*ZAS/(XPI*.92*XRHOLW)
                     ZEXP=2.*ZBS !XBS = 1.9 in ini_radar.f90 (bj tab 2.1 p24)
                     !dans ini_rain_ice.f90 :
-                    ZLBDA= ZLBS*(ZM)**ZLBEXS
-
+                    IF (GLIMA .AND. LSNOW_T_L) THEN
+                       IF (PT_RAY(JI,JEL,JAZ,JL,JH,JV)>-10.) THEN
+                          ZLBDA = MAX(MIN(XLBDAS_MAX, 10**(14.554-0.0423*(PT_RAY(JI,JEL,JAZ,JL,JH,JV)+273.15))),XLBDAS_MIN)
+                       ELSE
+                          ZLBDA = MAX(MIN(XLBDAS_MAX, 10**(6.226-0.0106*(PT_RAY(JI,JEL,JAZ,JL,JH,JV)+273.15))),XLBDAS_MIN)
+                       END IF
+                       ZN=ZLBS*ZM*ZLBDA**ZBS
+                    ELSE
+                       ZLBDA= ZLBS*(ZM)**ZLBEXS
+                       ZN=ZCCS*ZLBDA**ZCXS
+                    END IF
                     ! Rayleigh or Rayleigh-Gans or Rayleigh with 6th order for attenuation
                     IF(NDIFF==0.OR.NDIFF==3.OR.NDIFF==4) THEN
-                      ZREFLOC(1:2)=ZEQICE*.92**2*ZDMELT_FACT**2*1.E18*ZCCS*ZLBDA**(ZCXS-ZEXP)*MOMG(ZALPHAS,ZNUS,ZEXP)
+                      ZREFLOC(1:2)=ZEQICE*.92**2*ZDMELT_FACT**2*1.E18*ZN*ZLBDA**(ZEXP)*MOMG(ZALPHAS,ZNUS,ZEXP)
                       ZREFLOC(3)=0.
                       IF(LWREFL) THEN ! weighting by reflectivities
                         ZREFL(JI,JEL,JAZ,JL,JH,JV,IVDOP)=ZREFL(JI,JEL,JAZ,JL,JH,JV,IVDOP)&
                                -ZCS*SIN(PELEV(JI,JEL,JL,JV))*ZEQICE*.92**2*ZDMELT_FACT**2&
-                               *1.E18*ZCCS*ZLBDA**(ZCXS-ZEXP-ZDS)*MOMG(ZALPHAS,ZNUS,ZEXP+ZDS)
+                               *1.E18*ZN*ZLBDA**(ZEXP-ZDS)*MOMG(ZALPHAS,ZNUS,ZEXP+ZDS)
                       ELSE
-                        ZREFL(JI,JEL,JAZ,JL,JH,JV,IMAX)=ZREFL(JI,JEL,JAZ,JL,JH,JV,IMAX)+ZCCS*ZLBDA**ZCXS
+                        ZREFL(JI,JEL,JAZ,JL,JH,JV,IMAX)=ZREFL(JI,JEL,JAZ,JL,JH,JV,IMAX)+ZN
                         ZREFL(JI,JEL,JAZ,JL,JH,JV,IVDOP)=ZREFL(JI,JEL,JAZ,JL,JH,JV,IVDOP)&
                                                             -ZCS*SIN(PELEV(JI,JEL,JL,JV))&
-                                           *ZCCS*ZLBDA**(ZCXS-ZDS)*MOMG(ZALPHAS,ZNUS,ZDS)
+                                           *ZN*ZLBDA**(ZDS)*MOMG(ZALPHAS,ZNUS,ZDS)
                       END IF
                       IF(LATT) THEN
                         IF(NDIFF==0.OR.NDIFF==3) THEN
-                          ZAETMP(:)=ZCCS*ZLBDA**ZCXS*(ZDMELT_FACT*XPI**2/XLAM_RAD(JI)*AIMAG(ZQK)&
+                          ZAETMP(:)=ZN*(ZDMELT_FACT*XPI**2/XLAM_RAD(JI)*AIMAG(ZQK)&
                                     *MOMG(ZALPHAS,ZNUS,ZBS)/ZLBDA**ZBS)
                         ELSE
-                          ZAETMP(:)=ZCCS*ZLBDA**ZCXS*(ZDMELT_FACT*XPI**2/XLAM_RAD(JI)*AIMAG(ZQK)      &
+                          ZAETMP(:)=ZN*(ZDMELT_FACT*XPI**2/XLAM_RAD(JI)*AIMAG(ZQK)      &
                                     *MOMG(ZALPHAS,ZNUS,ZBS)/ZLBDA**ZBS                               &
                                     +ZDMELT_FACT**(5./3.)*XPI**4/15./XLAM_RAD(JI)**3                 &
                                     *AIMAG(ZQK**2*(ZQMI**4+27.*ZQMI**2+38.)                             &
@@ -1355,7 +1366,7 @@ DO JI=1,INBRAD
                       ZREFLOC(3)=180.E3/XPI*XLAM_RAD(JI)*ZRE_S22FMS11F
                       ZREFL(JI,JEL,JAZ,JL,JH,JV,IVDOP)=PVDOP_RAY(JI,JEL,JAZ,JL,JH,JV)*ZREFLOC(1) &
                                                         -ZCS*SIN(PELEV(JI,JEL,JL,JV))*ZREFLOC(1) &
-                                         *1.E18*(XLAM_RAD(JI)/XPI)**4/.93*ZCCS/4./ZLBDA**(3+ZDS)
+                                         *1.E18*(XLAM_RAD(JI)/XPI)**4/.93*(ZN*ZLBDA**(-ZCXS))/4./ZLBDA**(3+ZDS)
                       IF(LATT) THEN
                         ZAETMP(1)=ZIM_S22FT*XLAM_RAD(JI)*2
                         ZAETMP(2)=ZIM_S11FT*XLAM_RAD(JI)*2
@@ -1374,15 +1385,15 @@ DO JI=1,INBRAD
                         ZREFLOC(4)=ZREFLOC(4)+ZQBACK(1)*ZX(JJ)**(ZNUS-1.+2.*ZBS/3./ZALPHAS+ZDS/ZALPHAS)*ZW(JJ)
                         IF(LATT) ZAETMP(:)=ZAETMP(:)+ZQEXT(:)*ZX(JJ)**(ZNUS-1.+2.*ZBS/3./ZALPHAS)*ZW(JJ)
                       END DO ! ****** end loop Gauss-Laguerre quadrature
-                      ZREFLOC(1:2)=1.E18*(XLAM_RAD(JI)/XPI)**4*ZCCS*ZLBDA**(ZCXS-2.*ZBS/3.)/&
+                      ZREFLOC(1:2)=1.E18*(XLAM_RAD(JI)/XPI)**4*ZN*ZLBDA**(-2.*ZBS/3.)/&
                                   (4.*GAMMA(ZNUS)*.93)*ZDMELT_FACT**(2./3.)*ZREFLOC(1:2)               
                       ZREFL(JI,JEL,JAZ,JL,JH,JV,IVDOP)=ZREFL(JI,JEL,JAZ,JL,JH,JV,IVDOP)&
                                             +PVDOP_RAY(JI,JEL,JAZ,JL,JH,JV)*ZREFLOC(1) &
                                             -ZCS*SIN(PELEV(JI,JEL,JL,JV))*ZREFLOC(4)   &
-                                            *1.E18*(XLAM_RAD(JI)/XPI)**4*ZCCS          &
-                                            *ZLBDA**(ZCXS-2.*ZBS/3.-ZDS)/              &
+                                            *1.E18*(XLAM_RAD(JI)/XPI)**4*ZN          &
+                                            *ZLBDA**(2.*ZBS/3.-ZDS)/              &
                                             (4.*GAMMA(ZNUS)*.93)*ZDMELT_FACT**(2./3.)                              
-                      IF(LATT) ZAETMP(:)=ZAETMP(:)*XPI*ZCCS*ZLBDA**(ZCXS-2.*ZBS/3.)/(4.*GAMMA(ZNUS))&
+                      IF(LATT) ZAETMP(:)=ZAETMP(:)*XPI*ZN*ZLBDA**(-2.*ZBS/3.)/(4.*GAMMA(ZNUS))&
                                          *ZDMELT_FACT**(2./3.)            
                       ZRE_S22S11_S=0
                       ZIM_S22S11_S=0
diff --git a/src/MNH/rrcolss.f90 b/src/MNH/rrcolss.f90
index ccb2db545..00e7b6cc9 100644
--- a/src/MNH/rrcolss.f90
+++ b/src/MNH/rrcolss.f90
@@ -28,7 +28,7 @@ REAL, INTENT(IN) :: PESR      ! Efficiency of aggregates collecting rain
 REAL, INTENT(IN) :: PEXMASSR  ! Mass exponent of rain 
 REAL, INTENT(IN) :: PFALLS    ! Fall speed constant of aggregates
 REAL, INTENT(IN) :: PEXFALLS  ! Fall speed exponent of aggregates
-REAL, INTENT(IN) :: PFALLEXPS ! Fall speed exponential of aggregates ! Wurtz
+REAL, INTENT(IN) :: PFALLEXPS ! Fall speed exponential of aggregates (Thompson 2008)
 REAL, INTENT(IN) :: PFALLR    ! Fall speed constant of rain 
 REAL, INTENT(IN) :: PEXFALLR  ! Fall speed exponent of rain 
 REAL, INTENT(IN) :: PLBDASMAX ! Maximun slope of size distribution of aggregates
@@ -118,6 +118,7 @@ END INTERFACE
 !!      Original    8/11/95
 !!
 !  P. Wautelet 26/04/2019: replace non-standard FLOAT function by REAL function
+!  J. Wurtz       03/2022: new snow characteristics
 !
 !-------------------------------------------------------------------------------
 !
@@ -152,7 +153,7 @@ REAL, INTENT(IN) :: PESR      ! Efficiency of aggregates collecting rain
 REAL, INTENT(IN) :: PEXMASSR  ! Mass exponent of rain 
 REAL, INTENT(IN) :: PFALLS    ! Fall speed constant of aggregates
 REAL, INTENT(IN) :: PEXFALLS  ! Fall speed exponent of aggregates
-REAL, INTENT(IN) :: PFALLEXPS ! Fall speed exponential of aggregates ! Wurtz
+REAL, INTENT(IN) :: PFALLEXPS ! Fall speed exponential of aggregates  (Thompson 2008)
 REAL, INTENT(IN) :: PFALLR    ! Fall speed constant of rain 
 REAL, INTENT(IN) :: PEXFALLR  ! Fall speed exponent of rain 
 REAL, INTENT(IN) :: PLBDASMAX ! Maximun slope of size distribution of aggregates
@@ -279,11 +280,11 @@ DO JLBDAS = 1,SIZE(PRRCOLSS(:,:),1)
           DO JDR = 1,INR-1
             ZDR = ZDDCOLLR * REAL(JDR)
             ZCOLLR = ZCOLLR + (ZDS+ZDR)**2 * ZDR**PEXMASSR                     &
-                       * PESR * ABS(PFALLS*ZDS**PEXFALLS * EXP(-(PFALLEXPS*XDS)**PALPHAS)-PFALLR*ZDR**PEXFALLR) & ! Wurtz
+                       * PESR * ABS(PFALLS*ZDS**PEXFALLS * EXP(-(PFALLEXPS*XDS)**PALPHAS)-PFALLR*ZDR**PEXFALLR) &
                                       * GENERAL_GAMMA(PALPHAR,PNUR,ZLBDAR,ZDR)
           END DO
           ZCOLLDRMAX = (ZDS+ZDRMAX)**2 * ZDRMAX**PEXMASSR                      &
-                    * PESR * ABS(PFALLS*ZDS**PEXFALLS* EXP(-(PFALLEXPS*XDS)**PALPHAS)-PFALLR*ZDRMAX**PEXFALLR) & ! Wurtz
+                    * PESR * ABS(PFALLS*ZDS**PEXFALLS* EXP(-(PFALLEXPS*XDS)**PALPHAS)-PFALLR*ZDRMAX**PEXFALLR) &
                                    * GENERAL_GAMMA(PALPHAR,PNUR,ZLBDAR,ZDRMAX)
           ZCOLLR = (ZCOLLR + 0.5*ZCOLLDRMAX)*(ZDDCOLLR/ZDDSCALR)
 !
diff --git a/src/MNH/rscolrg.f90 b/src/MNH/rscolrg.f90
index a1508649b..210df0380 100644
--- a/src/MNH/rscolrg.f90
+++ b/src/MNH/rscolrg.f90
@@ -10,7 +10,7 @@
 INTERFACE
 !
       SUBROUTINE RSCOLRG( KND, PALPHAS, PZNUS, PALPHAR, PNUR,                &
-                         PESR, PEXMASSS, PFALLS, PEXFALLS, PFALLEXPS, PFALLR, PEXFALLR, & ! Wurtz
+                         PESR, PEXMASSS, PFALLS, PEXFALLS, PFALLEXPS, PFALLR, PEXFALLR, &
                          PLBDASMAX, PLBDARMAX, PLBDASMIN, PLBDARMIN,         &
                          PDINFTY, PRSCOLRG,PAG, PBS, PAS                     )
 !
@@ -28,7 +28,7 @@ REAL, INTENT(IN) :: PESR      ! Efficiency of the aggregates collecting rain
 REAL, INTENT(IN) :: PEXMASSS  ! Mass exponent of the aggregates
 REAL, INTENT(IN) :: PFALLS    ! Fall speed constant of the aggregates
 REAL, INTENT(IN) :: PEXFALLS  ! Fall speed exponent of the aggregates
-REAL, INTENT(IN) :: PFALLEXPS  ! Fall speed exponential constant of the aggregates ! Wurtz
+REAL, INTENT(IN) :: PFALLEXPS  ! Fall speed exponential constant of the aggregates
 REAL, INTENT(IN) :: PFALLR    ! Fall speed constant of rain 
 REAL, INTENT(IN) :: PEXFALLR  ! Fall speed exponent of rain 
 REAL, INTENT(IN) :: PLBDASMAX ! Maximun slope of size distribution of the aggregates
@@ -118,6 +118,7 @@ END INTERFACE
 !!      Original    8/11/95
 !!
 !  P. Wautelet 26/04/2019: replace non-standard FLOAT function by REAL function
+!  J. Wurtz       03/2022: new snow characteristics
 !
 !-------------------------------------------------------------------------------
 !
@@ -150,7 +151,7 @@ REAL, INTENT(IN) :: PESR      ! Efficiency of the aggregates collecting rain
 REAL, INTENT(IN) :: PEXMASSS  ! Mass exponent of the aggregates
 REAL, INTENT(IN) :: PFALLS    ! Fall speed constant of the aggregates
 REAL, INTENT(IN) :: PEXFALLS  ! Fall speed exponent of the aggregates
-REAL, INTENT(IN) :: PFALLEXPS ! Fall speed exponential constant of the aggregates ! Wurtz
+REAL, INTENT(IN) :: PFALLEXPS ! Fall speed exponential constant of the aggregates
 REAL, INTENT(IN) :: PFALLR    ! Fall speed constant of rain 
 REAL, INTENT(IN) :: PEXFALLR  ! Fall speed exponent of rain 
 REAL, INTENT(IN) :: PLBDASMAX ! Maximun slope of size distribution of the aggregates
@@ -273,12 +274,12 @@ DO JLBDAR = 1,SIZE(PRSCOLRG(:,:),1)
             ZDR = ZDDCOLLR * REAL(JDR) + ZDRMIN
             ZCOLLR = ZCOLLR + (ZDS+ZDR)**2                                     &
                        * GENERAL_GAMMA(PALPHAR,PNUR,ZLBDAR,ZDR)                &
-                         * PESR * ABS(PFALLS*ZDS**PEXFALLS*EXP(-(ZDS*PFALLEXPS)**PALPHAS)-PFALLR*ZDR**PEXFALLR) ! Wurtz
+                         * PESR * ABS(PFALLS*ZDS**PEXFALLS*EXP(-(ZDS*PFALLEXPS)**PALPHAS)-PFALLR*ZDR**PEXFALLR)
           END DO
           IF( ZDRMIN>0.0 ) THEN
             ZCOLLDRMIN = (ZDS+ZDRMIN)**2                                       &
                       * GENERAL_GAMMA(PALPHAR,PNUR,ZLBDAR,ZDRMIN)              &
-                      * PESR * ABS(PFALLS*ZDS**PEXFALLS*EXP(-(ZDS*PFALLEXPS)**PALPHAS)-PFALLR*ZDRMIN**PEXFALLR) ! Wurtz
+                      * PESR * ABS(PFALLS*ZDS**PEXFALLS*EXP(-(ZDS*PFALLEXPS)**PALPHAS)-PFALLR*ZDRMIN**PEXFALLR)
             ELSE
             ZCOLLDRMIN = 0.0
           END IF 
diff --git a/src/MNH/rzcolx.f90 b/src/MNH/rzcolx.f90
index 3887c24b7..d3977af63 100644
--- a/src/MNH/rzcolx.f90
+++ b/src/MNH/rzcolx.f90
@@ -11,7 +11,7 @@ INTERFACE
 !
       SUBROUTINE RZCOLX( KND, PALPHAX, PNUX, PALPHAZ, PNUZ,                  &
                          PEXZ, PEXMASSZ, PFALLX, PEXFALLX, PFALLEXPX,        &
-                         PFALLZ, PEXFALLZ, PFALLEXPZ,                        & ! Wurtz
+                         PFALLZ, PEXFALLZ, PFALLEXPZ,                        &
 		         PLBDAXMAX, PLBDAZMAX, PLBDAXMIN, PLBDAZMIN,         &
 		         PDINFTY, PRZCOLX                                    )
 !
@@ -30,10 +30,10 @@ REAL, INTENT(IN) :: PEXZ      ! Efficiency of specy X collecting specy Z
 REAL, INTENT(IN) :: PEXMASSZ  ! Mass exponent of specy Z
 REAL, INTENT(IN) :: PFALLX    ! Fall speed constant of specy X
 REAL, INTENT(IN) :: PEXFALLX  ! Fall speed exponent of specy X
-REAL, INTENT(IN) :: PFALLEXPX ! Fall speed exponential constant of specy X ! Wurtz
+REAL, INTENT(IN) :: PFALLEXPX ! Fall speed exponential constant of specy X
 REAL, INTENT(IN) :: PFALLZ    ! Fall speed constant of specy Z
 REAL, INTENT(IN) :: PEXFALLZ  ! Fall speed exponent of specy Z
-REAL, INTENT(IN) :: PFALLEXPZ ! Fall speed exponent of specy Z  ! Wurtz
+REAL, INTENT(IN) :: PFALLEXPZ ! Fall speed exponential constant of specy Z
 REAL, INTENT(IN) :: PLBDAXMAX ! Maximun slope of size distribution of specy X
 REAL, INTENT(IN) :: PLBDAZMAX ! Maximun slope of size distribution of specy Z
 REAL, INTENT(IN) :: PLBDAXMIN ! Minimun slope of size distribution of specy X
@@ -53,7 +53,7 @@ END INTERFACE
 !     ########################################################################
       SUBROUTINE RZCOLX( KND, PALPHAX, PNUX, PALPHAZ, PNUZ,                  &
                          PEXZ, PEXMASSZ, PFALLX, PEXFALLX, PFALLEXPX,        &
-                         PFALLZ, PEXFALLZ, PFALLEXPZ,                        & ! Wurtz
+                         PFALLZ, PEXFALLZ, PFALLEXPZ,                        &
 		         PLBDAXMAX, PLBDAZMAX, PLBDAXMIN, PLBDAZMIN,         &
 		         PDINFTY, PRZCOLX                                    )
 !     ########################################################################
@@ -125,6 +125,7 @@ END INTERFACE
 !!      Original    8/11/95
 !!
 !  P. Wautelet 26/04/2019: replace non-standard FLOAT function by REAL function
+!  J. Wurtz       03/2022: new snow characteristics
 !
 !-------------------------------------------------------------------------------
 !
@@ -156,10 +157,10 @@ REAL, INTENT(IN) :: PEXZ      ! Efficiency of specy X collecting specy Z
 REAL, INTENT(IN) :: PEXMASSZ  ! Mass exponent of specy Z
 REAL, INTENT(IN) :: PFALLX    ! Fall speed constant of specy X
 REAL, INTENT(IN) :: PEXFALLX  ! Fall speed exponent of specy X
-REAL, INTENT(IN) :: PFALLEXPX ! Fall speed exponential constant of specy X ! Wurtz
+REAL, INTENT(IN) :: PFALLEXPX ! Fall speed exponential constant of specy X
 REAL, INTENT(IN) :: PFALLZ    ! Fall speed constant of specy Z
 REAL, INTENT(IN) :: PEXFALLZ  ! Fall speed exponent of specy Z
-REAL, INTENT(IN) :: PFALLEXPZ ! Fall speed exponent of specy Z  ! Wurtz
+REAL, INTENT(IN) :: PFALLEXPZ ! Fall speed exponential constant of specy Z
 REAL, INTENT(IN) :: PLBDAXMAX ! Maximun slope of size distribution of specy X
 REAL, INTENT(IN) :: PLBDAZMAX ! Maximun slope of size distribution of specy Z
 REAL, INTENT(IN) :: PLBDAXMIN ! Minimun slope of size distribution of specy X
diff --git a/src/MNH/write_lfifm1_for_diag.f90 b/src/MNH/write_lfifm1_for_diag.f90
index 967b74d5d..7ca40d7a4 100644
--- a/src/MNH/write_lfifm1_for_diag.f90
+++ b/src/MNH/write_lfifm1_for_diag.f90
@@ -3926,20 +3926,20 @@ IF (LLIDAR) THEN
     ZTMP3(:,:,:,1)=ZSIG_DST(:,:,:,IACCMODE)
     SELECT CASE ( CCLOUD )
     CASE('KESS''ICE3','ICE4')
-      CALL LIDAR(CCLOUD, YVIEW, XALT_LIDAR, XWVL_LIDAR, XZZ, XRHODREF, XCLDFR, &
+      CALL LIDAR(CCLOUD, YVIEW, XALT_LIDAR, XWVL_LIDAR, XZZ, XRHODREF, ZTEMP, XCLDFR, &
                  XRT, ZWORK31, ZWORK32,                                        &
                  PDSTC=ZTMP1,                                                  &
                  PDSTD=ZTMP2,                                                  &
                  PDSTS=ZTMP3)
     CASE('C2R2')
-      CALL LIDAR(CCLOUD, YVIEW, XALT_LIDAR, XWVL_LIDAR, XZZ, XRHODREF, XCLDFR, &
+      CALL LIDAR(CCLOUD, YVIEW, XALT_LIDAR, XWVL_LIDAR, XZZ, XRHODREF, ZTEMP, XCLDFR, &
                  XRT, ZWORK31, ZWORK32,                                        &
                  PCT=XSVT(:,:,:,NSV_C2R2BEG+1:NSV_C2R2END),                    &
                  PDSTC=ZTMP1,                                                  &
                  PDSTD=ZTMP2,                                                  &
                  PDSTS=ZTMP3)
     CASE('C3R5')
-      CALL LIDAR(CCLOUD, YVIEW, XALT_LIDAR, XWVL_LIDAR, XZZ, XRHODREF, XCLDFR, &
+      CALL LIDAR(CCLOUD, YVIEW, XALT_LIDAR, XWVL_LIDAR, XZZ, XRHODREF, ZTEMP, XCLDFR, &
                  XRT, ZWORK31, ZWORK32,                                        &
                  PCT=XSVT(:,:,:,NSV_C2R2BEG+1:NSV_C1R3END-1),                  &
                  PDSTC=ZTMP1,                                                  &
@@ -3953,7 +3953,7 @@ IF (LLIDAR) THEN
        ZTMP4(:,:,:,3)=XSVT(:,:,:,NSV_LIMA_NR)
        ZTMP4(:,:,:,4)=XSVT(:,:,:,NSV_LIMA_NI)
 !
-       CALL LIDAR(CCLOUD, YVIEW, XALT_LIDAR, XWVL_LIDAR, XZZ, XRHODREF, MAX(XCLDFR,XICEFR),&
+       CALL LIDAR(CCLOUD, YVIEW, XALT_LIDAR, XWVL_LIDAR, XZZ, XRHODREF, ZTEMP, MAX(XCLDFR,XICEFR),&
             XRT, ZWORK31, ZWORK32,                                  &
             PCT=ZTMP4,                            &
             PDSTC=ZTMP1,                          &
@@ -3964,14 +3964,14 @@ IF (LLIDAR) THEN
   ELSE
     SELECT CASE ( CCLOUD )
     CASE('KESS','ICE3','ICE4')
-      CALL LIDAR(CCLOUD, YVIEW, XALT_LIDAR, XWVL_LIDAR, XZZ, XRHODREF, XCLDFR, &
+      CALL LIDAR(CCLOUD, YVIEW, XALT_LIDAR, XWVL_LIDAR, XZZ, XRHODREF, ZTEMP, XCLDFR, &
            XRT, ZWORK31, ZWORK32)
     CASE('C2R2')
-      CALL LIDAR(CCLOUD, YVIEW, XALT_LIDAR, XWVL_LIDAR, XZZ, XRHODREF, XCLDFR, &
+      CALL LIDAR(CCLOUD, YVIEW, XALT_LIDAR, XWVL_LIDAR, XZZ, XRHODREF, ZTEMP, XCLDFR, &
            XRT, ZWORK31, ZWORK32,                                  &
            PCT=XSVT(:,:,:,NSV_C2R2BEG+1:NSV_C2R2END))
     CASE('C3R5')
-      CALL LIDAR(CCLOUD, YVIEW, XALT_LIDAR, XWVL_LIDAR, XZZ, XRHODREF, XCLDFR, &
+      CALL LIDAR(CCLOUD, YVIEW, XALT_LIDAR, XWVL_LIDAR, XZZ, XRHODREF, ZTEMP, XCLDFR, &
            XRT, ZWORK31, ZWORK32,                                  &
            PCT=XSVT(:,:,:,NSV_C2R2BEG+1:NSV_C1R3END-1))
     CASE('LIMA')
@@ -3982,7 +3982,7 @@ IF (LLIDAR) THEN
        ZTMP4(:,:,:,3)=XSVT(:,:,:,NSV_LIMA_NR)
        ZTMP4(:,:,:,4)=XSVT(:,:,:,NSV_LIMA_NI)
 !
-       CALL LIDAR(CCLOUD, YVIEW, XALT_LIDAR, XWVL_LIDAR, XZZ, XRHODREF, MAX(XCLDFR,XICEFR),&
+       CALL LIDAR(CCLOUD, YVIEW, XALT_LIDAR, XWVL_LIDAR, XZZ, XRHODREF, ZTEMP, MAX(XCLDFR,XICEFR),&
             XRT, ZWORK31, ZWORK32,                                  &
             PCT=ZTMP4)
     END SELECT
-- 
GitLab