Index: LMDZ6/trunk/libf/phylmd/lmdz_blosno_ini.f90
===================================================================
--- LMDZ6/trunk/libf/phylmd/lmdz_blosno_ini.f90	(revision 6028)
+++ LMDZ6/trunk/libf/phylmd/lmdz_blosno_ini.f90	(revision 6028)
@@ -0,0 +1,84 @@
+module lmdz_blosno_ini
+
+   implicit none
+
+   real, save, protected :: RCPD, RV, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG, RPI
+   real, save, protected :: coef_sub_bs, fallv_bs, zeta_bs, c_esalt_bs
+   real, save, protected :: prt_bs, pbst_bs, qbst_bs, r_bs
+   integer, save, protected :: iflag_saltation_bs, iflag_sedim_bs, iflag_sublim_bs
+
+   !$OMP THREADPRIVATE(RCPD, RV, RLSTT, RLVTT, RLMLT, RVTMP2, RTT,RD,RG, RPI)
+   !$OMP THREADPRIVATE(coef_sub_bs, fallv_bs, r_bs, zeta_bs, c_esalt_bs)
+   !$OMP THREADPRIVATE(pbst_bs, prt_bs, qbst_bs)
+   !$OMP THREADPRIVATE(iflag_saltation_bs,iflag_sedim_bs, iflag_sublim_bs)
+
+   real, save, protected :: tbsmelt = 278.15    ! parameter to calculate melting fraction of BS sedimentation
+   real, save, protected :: taumeltbs0 = 600.0  ! Melting time scale of blowing snow at 273.15K
+   real, save, protected :: qbmin = 1.E-10      ! Minimum blowing snow specific content
+   !$OMP THREADPRIVATE(tbsmelt, taumeltbs0, qbmin)
+
+   real, save, protected :: tau_dens0_bs = 864000.      ! 10 days by default, in s
+   real, save, protected :: tau_densmin_bs = 21600.    ! 1/4 days according to in situ obs by C. Amory during blowing snow +
+   ! Marshall et al. 1999 (snow densification during rain)
+   real, save, protected :: tau_eqsalt_bs = 10.        ! Nemoto and Nishimura 2004 show that steady-state saltation is achieved within a time tau_eqsalt of about 10s
+   real, save, protected :: rhofresh_bs = 300.0       ! fresh snow density kg/m3
+   real, save, protected :: rhohard_bs = 450.0       ! hard snow density kg/m3
+   real, save, protected :: rhoice_bs = 920.0         ! ice density kg/m3
+   real, save, protected :: rhobs = 900.0               ! blowing snow density (kg/m3) following Bintanja et al. 2001 part I
+   !$OMP THREADPRIVATE(rhoice_bs, rhofresh_bs, rhohard_bs, tau_dens0_bs, tau_densmin_bs, tau_eqsalt_bs, rhobs)
+
+contains
+
+   subroutine blosno_ini(RCPD_in, RLSTT_in, RLVTT_in, RLMLT_in, &
+                               RVTMP2_in, RTT_in, RD_in, RG_in, RV_in, RPI_in)
+
+      USE ioipsl_getin_p_mod, ONLY: getin_p
+
+      real, intent(in) :: RCPD_in, RLSTT_in, RLVTT_in, RLMLT_in, RPI_in
+      real, intent(in) ::  RVTMP2_in, RTT_in, RD_in, RG_in, RV_in
+
+      RG = RG_in
+      RD = RD_in
+      RV = RV_in
+      RCPD = RCPD_in
+      RLVTT = RLVTT_in
+      RLSTT = RLSTT_in
+      RLMLT = RLMLT_in
+      RTT = RTT_in
+      RG = RG_in
+      RVTMP2 = RVTMP2_in
+      RPI = RPI_in
+
+      c_esalt_bs = 3.25
+      CALL getin_p('c_esalt_bs', c_esalt_bs)
+
+      qbst_bs = 0.001
+      CALL getin_p('qbst_bs', qbst_bs)
+
+      pbst_bs = 0.00003
+      CALL getin_p('pbst_bs', pbst_bs)
+
+      prt_bs = 0.00003
+      CALL getin_p('prt_bs', prt_bs)
+
+      zeta_bs = 1.
+      CALL getin_p('zeta_bs', zeta_bs)
+
+      fallv_bs = 0.5
+      CALL getin_p('fallv_bs', fallv_bs)
+
+      coef_sub_bs = 0.01
+      CALL getin_p('coef_sub_bs', coef_sub_bs)
+
+      iflag_sublim_bs = 1
+      CALL getin_p('iflag_sublim_bs', iflag_sublim_bs)
+
+      iflag_sedim_bs = 1
+      CALL getin_p('iflag_sedim_bs', iflag_sedim_bs)
+
+      r_bs = 50.0e-6
+      CALL getin_p('r_bs', r_bs)
+
+   end subroutine blosno_ini
+
+end module lmdz_blosno_ini
Index: LMDZ6/trunk/libf/phylmd/lmdz_blosno_sublimsedim.f90
===================================================================
--- LMDZ6/trunk/libf/phylmd/lmdz_blosno_sublimsedim.f90	(revision 6028)
+++ LMDZ6/trunk/libf/phylmd/lmdz_blosno_sublimsedim.f90	(revision 6028)
@@ -0,0 +1,292 @@
+module lmdz_blosno_sublimsedim
+
+contains
+
+!==============================================================================
+   subroutine blosno_sublimsedim(ngrid,nlay,dtime,ustar,temp,qv,qb,pplay,paprs,dtemp_bs,dqv_bs,dqb_bs,bsfl,precip_bs)
+
+!==============================================================================
+! Routine that calculates the sublimation, melting and sedimentation
+! of blowing snow
+! Reference: Vignon et al 2025, GMD, https://doi.org/10.5194/egusphere-2025-2871
+!
+! Etienne Vignon, October 2022
+!==============================================================================
+
+      use lmdz_blosno_ini, only: iflag_sublim_bs, iflag_sedim_bs, coef_sub_bs, RTT, RD, RG, fallv_bs
+      use lmdz_blosno_ini, only: qbmin, RCPD, RLSTT, RLMLT, RLVTT, RVTMP2, RV, RPI
+      use lmdz_blosno_ini, only: tbsmelt, taumeltbs0, rhobs, r_bs
+      USE lmdz_lscp_tools, only: calc_qsat_ecmwf
+
+      implicit none
+
+!++++++++++++++++++++++++++++++++++++++++++++++++++++
+! Declarations
+!++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+!INPUT
+!=====
+      integer, intent(in)                     :: ngrid ! number of horizontal grid points
+      integer, intent(in)                     :: nlay  ! number of vertical layers
+      real, intent(in)                        :: dtime ! physics time step [s]
+      real, intent(in), dimension(ngrid)      :: ustar ! surface friction velocity [m/s]
+      real, intent(in), dimension(ngrid, nlay) :: temp  ! temperature of the air [K]
+      real, intent(in), dimension(ngrid, nlay) :: qv    ! specific content of water [kg/kg]
+      real, intent(in), dimension(ngrid, nlay) :: qb    ! specific content of blowing snow [kg/kg]
+      real, intent(in), dimension(ngrid, nlay) :: pplay ! pressure at the middle of layers [Pa]
+      real, intent(in), dimension(ngrid, nlay + 1) :: paprs ! pressure at the layer bottom interface [Pa]
+
+! OUTPUT
+!========
+      real, intent(out), dimension(ngrid, nlay) :: dtemp_bs ! temperature tendency [K/s]
+      real, intent(out), dimension(ngrid, nlay) :: dqv_bs   ! water vapor tendency [kg/kg/s]
+      real, intent(out), dimension(ngrid, nlay) :: dqb_bs   ! blowing snow mass tendancy [kg/kg/s]
+      real, intent(out), dimension(ngrid, nlay + 1) :: bsfl   ! vertical profile of blowing snow vertical flux [kg/m2/s]
+      real, intent(out), dimension(ngrid)      :: precip_bs ! surface sedimentation flux of blowing snow [kg/s/m2]
+
+! LOCAL
+!======
+
+      integer                                  :: k, i, n
+      real                                     :: cpd, cpw, dqbsub, maxdqbsub, dqbmelt, zmair
+      real                                     :: dqsedim, precbs, dqvmelt, zmelt, taumeltbs
+      real                                     :: maxdqvmelt, rhoair, dz, dqbsedim
+      real                                     :: delta_p, b_p, a_p, c_p, c_sub, qvsub
+      real                                     :: p0, T0, Dv, Aprime, Bprime, Ka, radius0, radius
+      real, dimension(ngrid)                   :: ztemp, zqv, zqb, zpres, qsi, dqsi, qsl, dqsl, qzero, sedim
+      real, dimension(ngrid)                   :: zpaprsup, zpaprsdn, ztemp_up, zqb_up, zvelo
+      real, dimension(ngrid, nlay)              :: temp_seri, qb_seri, qv_seri, zz
+
+!++++++++++++++++++++++++++++++++++++++++++++++++++
+! Initialisation
+!++++++++++++++++++++++++++++++++++++++++++++++++++
+
+      qzero(:) = 0.
+      dtemp_bs(:, :) = 0.
+      dqv_bs(:, :) = 0.
+      dqb_bs(:, :) = 0.
+      zvelo(:) = 0.
+      sedim(:) = 0.
+      precip_bs(:) = 0.
+      bsfl(:, :) = 0.
+      zz(:, :) = 0.
+
+      p0 = 101325.0    ! ref pressure
+
+      DO k = 1, nlay
+         DO i = 1, ngrid
+            temp_seri(i, k) = temp(i, k)
+            qv_seri(i, k) = qv(i, k)
+            qb_seri(i, k) = qb(i, k)
+            zz(i, k) = zz(i, k) + (paprs(i, k) - paprs(i, k + 1))/(pplay(i, k)/RD/temp(i, k)*RG)
+         END DO
+      END DO
+
+! Sedimentation scheme
+!----------------------
+
+      IF (iflag_sedim_bs .GT. 0) THEN
+! begin of top-down loop
+         DO k = nlay, 1, -1
+
+            DO i = 1, ngrid
+               ztemp(i) = temp_seri(i, k)
+               zqv(i) = qv_seri(i, k)
+               zqb(i) = qb_seri(i, k)
+               zpres(i) = pplay(i, k)
+               zpaprsup(i) = paprs(i, k + 1)
+               zpaprsdn(i) = paprs(i, k)
+            END DO
+
+            ! thermalization of blowing snow precip coming from above
+            IF (k .LE. nlay - 1) THEN
+
+               DO i = 1, ngrid
+                  zmair = (zpaprsdn(i) - zpaprsup(i))/RG
+                  ! RVTMP2=rcpv/rcpd-1
+                  cpd = RCPD*(1.0 + RVTMP2*zqv(i))
+                  cpw = RCPD*RVTMP2
+                  ! zqb_up: blowing snow mass that has to be thermalized with
+                  ! layer's air so that precipitation at the ground has the
+                  ! same temperature as the lowermost layer
+                  zqb_up(i) = sedim(i)*dtime/zmair
+                  ztemp_up(i) = temp_seri(i, k + 1)
+
+                  ! t(i,k+1)+d_t(i,k+1): new temperature of the overlying layer
+                  ztemp(i) = (ztemp_up(i)*zqb_up(i)*cpw + cpd*ztemp(i)) &
+                             /(cpd + zqb_up(i)*cpw)
+               END DO
+
+            END IF
+
+            DO i = 1, ngrid
+
+               rhoair = zpres(i)/ztemp(i)/RD
+               dz = (zpaprsdn(i) - zpaprsup(i))/(rhoair*RG)
+               ! BS fall velocity assumed to be constant for now
+               zvelo(i) = fallv_bs
+               ! dqb/dt_sedim=1/rho(sedim/dz - rho w qb/dz)
+               ! implicit resolution
+               dqbsedim = (sedim(i)/rhoair/dz*dtime + zqb(i))/(1.+zvelo(i)*dtime/dz) - zqb(i)
+               ! flux and dqb update
+               zqb(i) = zqb(i) + dqbsedim
+               !sedim(i) = sedim(i)-dqbsedim/dtime*rhoair*dz
+               sedim(i) = rhoair*zvelo(i)*zqb(i)
+
+               ! variables update:
+               bsfl(i, k) = sedim(i)
+
+               qb_seri(i, k) = zqb(i)
+               qv_seri(i, k) = zqv(i)
+               temp_seri(i, k) = ztemp(i)
+
+            END DO ! Loop on ngrid
+
+         END DO ! vertical loop
+
+!surface bs flux
+         DO i = 1, ngrid
+            precip_bs(i) = sedim(i)
+         END DO
+
+      END IF
+
+!+++++++++++++++++++++++++++++++++++++++++++++++
+! Sublimation and melting
+!++++++++++++++++++++++++++++++++++++++++++++++
+      IF (iflag_sublim_bs .GT. 0) THEN
+
+         DO k = 1, nlay
+
+            DO i = 1, ngrid
+               ztemp(i) = temp_seri(i, k)
+               zqv(i) = qv_seri(i, k)
+               zqb(i) = qb_seri(i, k)
+               zpres(i) = pplay(i, k)
+               zpaprsup(i) = paprs(i, k + 1)
+               zpaprsdn(i) = paprs(i, k)
+
+               ! computation of blowing snow mean radius. Either constant value = r_bs if >0
+               ! or use of the height-dependent formulation from Sauter et al. 2013 and Saigger et al. 2024
+               IF (r_bs .GT. 0.) THEN
+                  radius = r_bs
+               ELSE
+                  radius0 = 0.5*(ustar(i)*(7.8e-6)/0.036 + 31.e-6)
+                  radius = radius0*zz(i, k)**(-0.258)
+               END IF
+
+            END DO
+
+            ! calulation saturation specific humidity
+            CALL CALC_QSAT_ECMWF(ngrid, ztemp(:), qzero(:), zpres(:), RTT, 2, .false., qsi(:), dqsi(:))
+
+            DO i = 1, ngrid
+
+               rhoair = zpres(i)/ztemp(i)/RD
+               dz = (zpaprsdn(i) - zpaprsup(i))/(rhoair*RG)
+               ! BS fall velocity assumed to be constant for now
+               zvelo(i) = fallv_bs
+
+               IF (ztemp(i) .GT. RTT) THEN
+
+                  ! if temperature is positive, we assume that part of the blowing snow
+                  ! already present  melts and evaporates with a typical time
+                  ! constant taumeltbs
+
+                  taumeltbs = taumeltbs0*exp(-max(0., (ztemp(i) - RTT)/(tbsmelt - RTT)))
+                  dqvmelt = min(zqb(i), -1.*zqb(i)*(exp(-dtime/taumeltbs) - 1.))
+                  maxdqvmelt = max(RCPD*(1.0 + RVTMP2*(zqv(i) + zqb(i)))*(ztemp(i) - RTT)/(RLMLT + RLVTT), 0.)
+                  dqvmelt = min(dqvmelt, maxdqvmelt)
+                  ! qv update, melting + evaporation
+                  zqv(i) = zqv(i) + dqvmelt
+                  ! temp update melting + evaporation
+                  ztemp(i) = ztemp(i) - dqvmelt*(RLMLT + RLVTT)/RCPD/(1.0 + RVTMP2*(zqv(i) + zqb(i)))
+                  ! qb update melting + evaporation
+                  zqb(i) = zqb(i) - dqvmelt
+
+               ELSE
+                  ! negative celcius temperature
+                  ! Sublimation scheme
+
+                  ! Sublimation formulation for ice crystals from Pruppacher & Klett, Rutledge & Hobbs 1983
+                  ! assuming monodispered crystal distrib
+                  ! dqb/dt_sub=-coef_sub_bs*(1-qv/qsi)*nc*8*radius/(Aprime+Bprime)
+                  ! assuming Mi=rhobs*4/3*pi*radius**3
+                  ! qb=nc*Mi -> nc=qb/Mi
+                  ! dqb/dt_sub=-coef_sub_bs*(1-qv/qsi)*6*qb/(rhobs*pi*radius**2)/(Aprime+Bprime)
+                  ! dqb/dt_sub=-c_sub(1-qv/qsi)*qb
+                  ! c_sub=coef_sub_bs*6/(rhobs*pi*radius**2)/(Aprime+Bprime)
+                  !
+                  ! Note the strong coupling between specific contents of water vapor and blowing snow during sublimation
+                  ! equations dqv/dt_sub and dqb/dt_sub must be solved jointly to prevent irrealistic increase of water vapor
+                  ! at typical physics time steps
+                  ! we thus solve the differential equation using an implicit scheme for both qb and qv
+
+                  ! we do not consider deposition, only sublimation
+                  IF (zqv(i) .LT. qsi(i)) THEN
+                     rhoair = zpres(i)/ztemp(i)/RD
+                     Dv = 0.0001*0.211*(p0/zpres(i))*((ztemp(i)/RTT)**1.94)           ! water vapor diffusivity in air, SI
+                     Ka = (5.69 + 0.017*(ztemp(i) - RTT))*1.e-5*100.*4.184                ! thermal conductivity of the air, SI
+                     Aprime = RLSTT/Ka/ztemp(i)*(RLSTT/RV/ztemp(i) - 1.)
+                     Bprime = 1./(rhoair*Dv*qsi(i))
+                     c_sub = coef_sub_bs*3./(rhobs*radius**2)/(Aprime + Bprime)
+                     c_p = -zqb(i)
+                     b_p = 1.+c_sub*dtime - c_sub*dtime/qsi(i)*zqb(i) - c_sub*dtime/qsi(i)*zqv(i)
+                     a_p = c_sub*dtime/qsi(i)
+                     delta_p = (b_p**2) - 4.*a_p*c_p
+                     dqbsub = (-b_p + sqrt(delta_p))/(2.*a_p) - zqb(i)
+                     dqbsub = MIN(0.0, MAX(dqbsub, -zqb(i)))
+                     ! Sublimation limit: we ensure that the whole mesh does not exceed saturation wrt ice
+                     maxdqbsub = MAX(0.0, qsi(i) - zqv(i))
+                     dqbsub = MAX(dqbsub, -maxdqbsub)
+                  ELSE
+                     dqbsub = 0.
+                  END IF
+
+                  ! vapor, temperature, precip fluxes update following sublimation
+                  zqv(i) = zqv(i) - dqbsub
+                  zqb(i) = zqb(i) + dqbsub
+                  ztemp(i) = ztemp(i) + dqbsub*RLSTT/RCPD/(1.0 + RVTMP2*(zqv(i) + zqb(i)))
+
+               END IF
+
+               ! if qb<qbmin, sublimate or melt and evaporate qb
+               ! see Gerber et al. 2023, JGR Atmos for the choice of qbmin
+
+               IF (zqb(i) .LT. qbmin) THEN
+                  zqv(i) = zqv(i) + zqb(i)
+                  IF (ztemp(i) .LT. RTT) THEN
+                     ztemp(i) = ztemp(i) - zqb(i)*RLSTT/RCPD/(1.0 + RVTMP2*(zqv(i)))
+                  ELSE
+                     ztemp(i) = ztemp(i) - zqb(i)*(RLVTT + RLMLT)/RCPD/(1.0 + RVTMP2*(zqv(i)))
+                  END IF
+                  zqb(i) = 0.
+               END IF
+
+               ! variables update
+               temp_seri(i, k) = ztemp(i)
+               qv_seri(i, k) = zqv(i)
+               qb_seri(i, k) = zqb(i)
+            END DO
+         END DO
+
+      END IF
+
+! OUTPUTS
+!++++++++++
+
+! 3D variables
+      DO k = 1, nlay
+         DO i = 1, ngrid
+            dqb_bs(i, k) = qb_seri(i, k) - qb(i, k)
+            dqv_bs(i, k) = qv_seri(i, k) - qv(i, k)
+            dtemp_bs(i, k) = temp_seri(i, k) - temp(i, k)
+         END DO
+      END DO
+
+      return
+
+   end subroutine blosno_sublimsedim
+!==============================================================================
+
+end module lmdz_blosno_sublimsedim
Index: LMDZ6/trunk/libf/phylmd/lmdz_blowing_snow_ini.f90
===================================================================
--- LMDZ6/trunk/libf/phylmd/lmdz_blowing_snow_ini.f90	(revision 6027)
+++ 	(revision )
@@ -1,84 +1,0 @@
-module lmdz_blowing_snow_ini
-
-   implicit none
-
-   real, save, protected :: RCPD, RV, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG, RPI
-   real, save, protected :: coef_sub_bs, fallv_bs, zeta_bs, c_esalt_bs
-   real, save, protected :: prt_bs, pbst_bs, qbst_bs, r_bs
-   integer, save, protected :: iflag_saltation_bs, iflag_sedim_bs, iflag_sublim_bs
-
-   !$OMP THREADPRIVATE(RCPD, RV, RLSTT, RLVTT, RLMLT, RVTMP2, RTT,RD,RG, RPI)
-   !$OMP THREADPRIVATE(coef_sub_bs, fallv_bs, r_bs, zeta_bs, c_esalt_bs)
-   !$OMP THREADPRIVATE(pbst_bs, prt_bs, qbst_bs)
-   !$OMP THREADPRIVATE(iflag_saltation_bs,iflag_sedim_bs, iflag_sublim_bs)
-
-   real, save, protected :: tbsmelt = 278.15    ! parameter to calculate melting fraction of BS sedimentation
-   real, save, protected :: taumeltbs0 = 600.0  ! Melting time scale of blowing snow at 273.15K
-   real, save, protected :: qbmin = 1.E-10      ! Minimum blowing snow specific content
-   !$OMP THREADPRIVATE(tbsmelt, taumeltbs0, qbmin)
-
-   real, save, protected :: tau_dens0_bs = 864000.      ! 10 days by default, in s
-   real, save, protected :: tau_densmin_bs = 21600.    ! 1/4 days according to in situ obs by C. Amory during blowing snow +
-   ! Marshall et al. 1999 (snow densification during rain)
-   real, save, protected :: tau_eqsalt_bs = 10.        ! Nemoto and Nishimura 2004 show that steady-state saltation is achieved within a time tau_eqsalt of about 10s
-   real, save, protected :: rhofresh_bs = 300.0       ! fresh snow density kg/m3
-   real, save, protected :: rhohard_bs = 450.0       ! hard snow density kg/m3
-   real, save, protected :: rhoice_bs = 920.0         ! ice density kg/m3
-   real, save, protected :: rhobs = 900.0               ! blowing snow density (kg/m3) following Bintanja et al. 2001 part I
-   !$OMP THREADPRIVATE(rhoice_bs, rhofresh_bs, rhohard_bs, tau_dens0_bs, tau_densmin_bs, tau_eqsalt_bs, rhobs)
-
-contains
-
-   subroutine blowing_snow_ini(RCPD_in, RLSTT_in, RLVTT_in, RLMLT_in, &
-                               RVTMP2_in, RTT_in, RD_in, RG_in, RV_in, RPI_in)
-
-      USE ioipsl_getin_p_mod, ONLY: getin_p
-
-      real, intent(in) :: RCPD_in, RLSTT_in, RLVTT_in, RLMLT_in, RPI_in
-      real, intent(in) ::  RVTMP2_in, RTT_in, RD_in, RG_in, RV_in
-
-      RG = RG_in
-      RD = RD_in
-      RV = RV_in
-      RCPD = RCPD_in
-      RLVTT = RLVTT_in
-      RLSTT = RLSTT_in
-      RLMLT = RLMLT_in
-      RTT = RTT_in
-      RG = RG_in
-      RVTMP2 = RVTMP2_in
-      RPI = RPI_in
-
-      c_esalt_bs = 3.25
-      CALL getin_p('c_esalt_bs', c_esalt_bs)
-
-      qbst_bs = 0.001
-      CALL getin_p('qbst_bs', qbst_bs)
-
-      pbst_bs = 0.00003
-      CALL getin_p('pbst_bs', pbst_bs)
-
-      prt_bs = 0.00003
-      CALL getin_p('prt_bs', prt_bs)
-
-      zeta_bs = 1.
-      CALL getin_p('zeta_bs', zeta_bs)
-
-      fallv_bs = 0.5
-      CALL getin_p('fallv_bs', fallv_bs)
-
-      coef_sub_bs = 0.01
-      CALL getin_p('coef_sub_bs', coef_sub_bs)
-
-      iflag_sublim_bs = 1
-      CALL getin_p('iflag_sublim_bs', iflag_sublim_bs)
-
-      iflag_sedim_bs = 1
-      CALL getin_p('iflag_sedim_bs', iflag_sedim_bs)
-
-      r_bs = 50.0e-6
-      CALL getin_p('r_bs', r_bs)
-
-   end subroutine blowing_snow_ini
-
-end module lmdz_blowing_snow_ini
Index: LMDZ6/trunk/libf/phylmd/lmdz_blowing_snow_sublim_sedim.f90
===================================================================
--- LMDZ6/trunk/libf/phylmd/lmdz_blowing_snow_sublim_sedim.f90	(revision 6027)
+++ 	(revision )
@@ -1,292 +1,0 @@
-module lmdz_blowing_snow_sublim_sedim
-
-contains
-
-!==============================================================================
-   subroutine blowing_snow_sublim_sedim(ngrid,nlay,dtime,ustar,temp,qv,qb,pplay,paprs,dtemp_bs,dqv_bs,dqb_bs,bsfl,precip_bs)
-
-!==============================================================================
-! Routine that calculates the sublimation, melting and sedimentation
-! of blowing snow
-! Reference: Vignon et al 2025, GMD, https://doi.org/10.5194/egusphere-2025-2871
-!
-! Etienne Vignon, October 2022
-!==============================================================================
-
-      use lmdz_blowing_snow_ini, only: iflag_sublim_bs, iflag_sedim_bs, coef_sub_bs, RTT, RD, RG, fallv_bs
-      use lmdz_blowing_snow_ini, only: qbmin, RCPD, RLSTT, RLMLT, RLVTT, RVTMP2, RV, RPI
-      use lmdz_blowing_snow_ini, only: tbsmelt, taumeltbs0, rhobs, r_bs
-      USE lmdz_lscp_tools, only: calc_qsat_ecmwf
-
-      implicit none
-
-!++++++++++++++++++++++++++++++++++++++++++++++++++++
-! Declarations
-!++++++++++++++++++++++++++++++++++++++++++++++++++++
-
-!INPUT
-!=====
-      integer, intent(in)                     :: ngrid ! number of horizontal grid points
-      integer, intent(in)                     :: nlay  ! number of vertical layers
-      real, intent(in)                        :: dtime ! physics time step [s]
-      real, intent(in), dimension(ngrid)      :: ustar ! surface friction velocity [m/s]
-      real, intent(in), dimension(ngrid, nlay) :: temp  ! temperature of the air [K]
-      real, intent(in), dimension(ngrid, nlay) :: qv    ! specific content of water [kg/kg]
-      real, intent(in), dimension(ngrid, nlay) :: qb    ! specific content of blowing snow [kg/kg]
-      real, intent(in), dimension(ngrid, nlay) :: pplay ! pressure at the middle of layers [Pa]
-      real, intent(in), dimension(ngrid, nlay + 1) :: paprs ! pressure at the layer bottom interface [Pa]
-
-! OUTPUT
-!========
-      real, intent(out), dimension(ngrid, nlay) :: dtemp_bs ! temperature tendency [K/s]
-      real, intent(out), dimension(ngrid, nlay) :: dqv_bs   ! water vapor tendency [kg/kg/s]
-      real, intent(out), dimension(ngrid, nlay) :: dqb_bs   ! blowing snow mass tendancy [kg/kg/s]
-      real, intent(out), dimension(ngrid, nlay + 1) :: bsfl   ! vertical profile of blowing snow vertical flux [kg/m2/s]
-      real, intent(out), dimension(ngrid)      :: precip_bs ! surface sedimentation flux of blowing snow [kg/s/m2]
-
-! LOCAL
-!======
-
-      integer                                  :: k, i, n
-      real                                     :: cpd, cpw, dqbsub, maxdqbsub, dqbmelt, zmair
-      real                                     :: dqsedim, precbs, dqvmelt, zmelt, taumeltbs
-      real                                     :: maxdqvmelt, rhoair, dz, dqbsedim
-      real                                     :: delta_p, b_p, a_p, c_p, c_sub, qvsub
-      real                                     :: p0, T0, Dv, Aprime, Bprime, Ka, radius0, radius
-      real, dimension(ngrid)                   :: ztemp, zqv, zqb, zpres, qsi, dqsi, qsl, dqsl, qzero, sedim
-      real, dimension(ngrid)                   :: zpaprsup, zpaprsdn, ztemp_up, zqb_up, zvelo
-      real, dimension(ngrid, nlay)              :: temp_seri, qb_seri, qv_seri, zz
-
-!++++++++++++++++++++++++++++++++++++++++++++++++++
-! Initialisation
-!++++++++++++++++++++++++++++++++++++++++++++++++++
-
-      qzero(:) = 0.
-      dtemp_bs(:, :) = 0.
-      dqv_bs(:, :) = 0.
-      dqb_bs(:, :) = 0.
-      zvelo(:) = 0.
-      sedim(:) = 0.
-      precip_bs(:) = 0.
-      bsfl(:, :) = 0.
-      zz(:, :) = 0.
-
-      p0 = 101325.0    ! ref pressure
-
-      DO k = 1, nlay
-         DO i = 1, ngrid
-            temp_seri(i, k) = temp(i, k)
-            qv_seri(i, k) = qv(i, k)
-            qb_seri(i, k) = qb(i, k)
-            zz(i, k) = zz(i, k) + (paprs(i, k) - paprs(i, k + 1))/(pplay(i, k)/RD/temp(i, k)*RG)
-         END DO
-      END DO
-
-! Sedimentation scheme
-!----------------------
-
-      IF (iflag_sedim_bs .GT. 0) THEN
-! begin of top-down loop
-         DO k = nlay, 1, -1
-
-            DO i = 1, ngrid
-               ztemp(i) = temp_seri(i, k)
-               zqv(i) = qv_seri(i, k)
-               zqb(i) = qb_seri(i, k)
-               zpres(i) = pplay(i, k)
-               zpaprsup(i) = paprs(i, k + 1)
-               zpaprsdn(i) = paprs(i, k)
-            END DO
-
-            ! thermalization of blowing snow precip coming from above
-            IF (k .LE. nlay - 1) THEN
-
-               DO i = 1, ngrid
-                  zmair = (zpaprsdn(i) - zpaprsup(i))/RG
-                  ! RVTMP2=rcpv/rcpd-1
-                  cpd = RCPD*(1.0 + RVTMP2*zqv(i))
-                  cpw = RCPD*RVTMP2
-                  ! zqb_up: blowing snow mass that has to be thermalized with
-                  ! layer's air so that precipitation at the ground has the
-                  ! same temperature as the lowermost layer
-                  zqb_up(i) = sedim(i)*dtime/zmair
-                  ztemp_up(i) = temp_seri(i, k + 1)
-
-                  ! t(i,k+1)+d_t(i,k+1): new temperature of the overlying layer
-                  ztemp(i) = (ztemp_up(i)*zqb_up(i)*cpw + cpd*ztemp(i)) &
-                             /(cpd + zqb_up(i)*cpw)
-               END DO
-
-            END IF
-
-            DO i = 1, ngrid
-
-               rhoair = zpres(i)/ztemp(i)/RD
-               dz = (zpaprsdn(i) - zpaprsup(i))/(rhoair*RG)
-               ! BS fall velocity assumed to be constant for now
-               zvelo(i) = fallv_bs
-               ! dqb/dt_sedim=1/rho(sedim/dz - rho w qb/dz)
-               ! implicit resolution
-               dqbsedim = (sedim(i)/rhoair/dz*dtime + zqb(i))/(1.+zvelo(i)*dtime/dz) - zqb(i)
-               ! flux and dqb update
-               zqb(i) = zqb(i) + dqbsedim
-               !sedim(i) = sedim(i)-dqbsedim/dtime*rhoair*dz
-               sedim(i) = rhoair*zvelo(i)*zqb(i)
-
-               ! variables update:
-               bsfl(i, k) = sedim(i)
-
-               qb_seri(i, k) = zqb(i)
-               qv_seri(i, k) = zqv(i)
-               temp_seri(i, k) = ztemp(i)
-
-            END DO ! Loop on ngrid
-
-         END DO ! vertical loop
-
-!surface bs flux
-         DO i = 1, ngrid
-            precip_bs(i) = sedim(i)
-         END DO
-
-      END IF
-
-!+++++++++++++++++++++++++++++++++++++++++++++++
-! Sublimation and melting
-!++++++++++++++++++++++++++++++++++++++++++++++
-      IF (iflag_sublim_bs .GT. 0) THEN
-
-         DO k = 1, nlay
-
-            DO i = 1, ngrid
-               ztemp(i) = temp_seri(i, k)
-               zqv(i) = qv_seri(i, k)
-               zqb(i) = qb_seri(i, k)
-               zpres(i) = pplay(i, k)
-               zpaprsup(i) = paprs(i, k + 1)
-               zpaprsdn(i) = paprs(i, k)
-
-               ! computation of blowing snow mean radius. Either constant value = r_bs if >0
-               ! or use of the height-dependent formulation from Sauter et al. 2013 and Saigger et al. 2024
-               IF (r_bs .GT. 0.) THEN
-                  radius = r_bs
-               ELSE
-                  radius0 = 0.5*(ustar(i)*(7.8e-6)/0.036 + 31.e-6)
-                  radius = radius0*zz(i, k)**(-0.258)
-               END IF
-
-            END DO
-
-            ! calulation saturation specific humidity
-            CALL CALC_QSAT_ECMWF(ngrid, ztemp(:), qzero(:), zpres(:), RTT, 2, .false., qsi(:), dqsi(:))
-
-            DO i = 1, ngrid
-
-               rhoair = zpres(i)/ztemp(i)/RD
-               dz = (zpaprsdn(i) - zpaprsup(i))/(rhoair*RG)
-               ! BS fall velocity assumed to be constant for now
-               zvelo(i) = fallv_bs
-
-               IF (ztemp(i) .GT. RTT) THEN
-
-                  ! if temperature is positive, we assume that part of the blowing snow
-                  ! already present  melts and evaporates with a typical time
-                  ! constant taumeltbs
-
-                  taumeltbs = taumeltbs0*exp(-max(0., (ztemp(i) - RTT)/(tbsmelt - RTT)))
-                  dqvmelt = min(zqb(i), -1.*zqb(i)*(exp(-dtime/taumeltbs) - 1.))
-                  maxdqvmelt = max(RCPD*(1.0 + RVTMP2*(zqv(i) + zqb(i)))*(ztemp(i) - RTT)/(RLMLT + RLVTT), 0.)
-                  dqvmelt = min(dqvmelt, maxdqvmelt)
-                  ! qv update, melting + evaporation
-                  zqv(i) = zqv(i) + dqvmelt
-                  ! temp update melting + evaporation
-                  ztemp(i) = ztemp(i) - dqvmelt*(RLMLT + RLVTT)/RCPD/(1.0 + RVTMP2*(zqv(i) + zqb(i)))
-                  ! qb update melting + evaporation
-                  zqb(i) = zqb(i) - dqvmelt
-
-               ELSE
-                  ! negative celcius temperature
-                  ! Sublimation scheme
-
-                  ! Sublimation formulation for ice crystals from Pruppacher & Klett, Rutledge & Hobbs 1983
-                  ! assuming monodispered crystal distrib
-                  ! dqb/dt_sub=-coef_sub_bs*(1-qv/qsi)*nc*8*radius/(Aprime+Bprime)
-                  ! assuming Mi=rhobs*4/3*pi*radius**3
-                  ! qb=nc*Mi -> nc=qb/Mi
-                  ! dqb/dt_sub=-coef_sub_bs*(1-qv/qsi)*6*qb/(rhobs*pi*radius**2)/(Aprime+Bprime)
-                  ! dqb/dt_sub=-c_sub(1-qv/qsi)*qb
-                  ! c_sub=coef_sub_bs*6/(rhobs*pi*radius**2)/(Aprime+Bprime)
-                  !
-                  ! Note the strong coupling between specific contents of water vapor and blowing snow during sublimation
-                  ! equations dqv/dt_sub and dqb/dt_sub must be solved jointly to prevent irrealistic increase of water vapor
-                  ! at typical physics time steps
-                  ! we thus solve the differential equation using an implicit scheme for both qb and qv
-
-                  ! we do not consider deposition, only sublimation
-                  IF (zqv(i) .LT. qsi(i)) THEN
-                     rhoair = zpres(i)/ztemp(i)/RD
-                     Dv = 0.0001*0.211*(p0/zpres(i))*((ztemp(i)/RTT)**1.94)           ! water vapor diffusivity in air, SI
-                     Ka = (5.69 + 0.017*(ztemp(i) - RTT))*1.e-5*100.*4.184                ! thermal conductivity of the air, SI
-                     Aprime = RLSTT/Ka/ztemp(i)*(RLSTT/RV/ztemp(i) - 1.)
-                     Bprime = 1./(rhoair*Dv*qsi(i))
-                     c_sub = coef_sub_bs*3./(rhobs*radius**2)/(Aprime + Bprime)
-                     c_p = -zqb(i)
-                     b_p = 1.+c_sub*dtime - c_sub*dtime/qsi(i)*zqb(i) - c_sub*dtime/qsi(i)*zqv(i)
-                     a_p = c_sub*dtime/qsi(i)
-                     delta_p = (b_p**2) - 4.*a_p*c_p
-                     dqbsub = (-b_p + sqrt(delta_p))/(2.*a_p) - zqb(i)
-                     dqbsub = MIN(0.0, MAX(dqbsub, -zqb(i)))
-                     ! Sublimation limit: we ensure that the whole mesh does not exceed saturation wrt ice
-                     maxdqbsub = MAX(0.0, qsi(i) - zqv(i))
-                     dqbsub = MAX(dqbsub, -maxdqbsub)
-                  ELSE
-                     dqbsub = 0.
-                  END IF
-
-                  ! vapor, temperature, precip fluxes update following sublimation
-                  zqv(i) = zqv(i) - dqbsub
-                  zqb(i) = zqb(i) + dqbsub
-                  ztemp(i) = ztemp(i) + dqbsub*RLSTT/RCPD/(1.0 + RVTMP2*(zqv(i) + zqb(i)))
-
-               END IF
-
-               ! if qb<qbmin, sublimate or melt and evaporate qb
-               ! see Gerber et al. 2023, JGR Atmos for the choice of qbmin
-
-               IF (zqb(i) .LT. qbmin) THEN
-                  zqv(i) = zqv(i) + zqb(i)
-                  IF (ztemp(i) .LT. RTT) THEN
-                     ztemp(i) = ztemp(i) - zqb(i)*RLSTT/RCPD/(1.0 + RVTMP2*(zqv(i)))
-                  ELSE
-                     ztemp(i) = ztemp(i) - zqb(i)*(RLVTT + RLMLT)/RCPD/(1.0 + RVTMP2*(zqv(i)))
-                  END IF
-                  zqb(i) = 0.
-               END IF
-
-               ! variables update
-               temp_seri(i, k) = ztemp(i)
-               qv_seri(i, k) = zqv(i)
-               qb_seri(i, k) = zqb(i)
-            END DO
-         END DO
-
-      END IF
-
-! OUTPUTS
-!++++++++++
-
-! 3D variables
-      DO k = 1, nlay
-         DO i = 1, ngrid
-            dqb_bs(i, k) = qb_seri(i, k) - qb(i, k)
-            dqv_bs(i, k) = qv_seri(i, k) - qv(i, k)
-            dtemp_bs(i, k) = temp_seri(i, k) - temp(i, k)
-         END DO
-      END DO
-
-      return
-
-   end subroutine blowing_snow_sublim_sedim
-!==============================================================================
-
-end module lmdz_blowing_snow_sublim_sedim
Index: LMDZ6/trunk/libf/phylmd/lmdz_call_blosno.f90
===================================================================
--- LMDZ6/trunk/libf/phylmd/lmdz_call_blosno.f90	(revision 6028)
+++ LMDZ6/trunk/libf/phylmd/lmdz_call_blosno.f90	(revision 6028)
@@ -0,0 +1,69 @@
+!$gpum horizontal klon ngrid
+module lmdz_call_blosno
+
+contains
+
+   subroutine call_blosno(ngrid, nlay, abortphy, flag_inhib_tend, itap, dtime, ustar, temp, q, qb, pplay, paprs, &
+                                dtemp_bs, dq_bs, dqb_bs, bsfl, precip_bs)
+
+!==================================================================================================
+! call_blosno in the main interface between the LMDZ physics monitor physiq_mod
+! and the blowing snow parameterizations
+!
+! contact: Etienne Vignon etienne.vignon@lmd.ipsl.fr
+!================================================================================================
+
+      use lmdz_blosno_sublimsedim, only: blosno_sublimsedim
+      USE add_phys_tend_mod, ONLY: add_phys_tend, prt_enerbil
+
+      implicit none
+!======================================================
+! Declarations
+!======================================================
+
+!INPUT
+!=====
+
+      INTEGER, INTENT(IN)                     :: abortphy, flag_inhib_tend, itap ! flag and physics time counter for add_phys_tend
+      integer, intent(in)                     :: ngrid ! number of horizontal grid points
+      integer, intent(in)                     :: nlay  ! number of vertical layers
+      real, intent(in)                        :: dtime ! physics time step [s]
+      real, intent(in), dimension(ngrid)      :: ustar ! surface friction velocity [m/s]
+      real, intent(in), dimension(ngrid, nlay) :: temp  ! temperature of the air [K]
+      real, intent(in), dimension(ngrid, nlay) :: q     ! specific content of water [kg/kg]
+      real, intent(in), dimension(ngrid, nlay) :: qb    ! specific content of blowing snow [kg/kg]
+      real, intent(in), dimension(ngrid, nlay) :: pplay ! pressure at the middle of layers [Pa]
+      real, intent(in), dimension(ngrid, nlay + 1) :: paprs ! pressure at the layer bottom interface [Pa]
+
+! OUTPUT
+!========
+      real, intent(out), dimension(ngrid, nlay) :: dtemp_bs ! temperature tendency [K/s]
+      real, intent(out), dimension(ngrid, nlay) :: dq_bs    ! water vapor tendency [kg/kg/s]
+      real, intent(out), dimension(ngrid, nlay) :: dqb_bs   ! blowing snow mass tendancy [kg/kg/s]
+      real, intent(out), dimension(ngrid, nlay + 1) :: bsfl   ! vertical profile of blowing snow vertical flux [kg/m2/s]
+      real, intent(out), dimension(ngrid)      :: precip_bs ! surface sedimentation flux of blowing snow [kg/s/m2]
+
+! LOCAL
+!========
+      real, dimension(ngrid, nlay)              :: du0, dv0, dql0, dqi0
+
+!=================================================================
+! Call to main routine of blowing snow sublimation and sedim.
+!=================================================================
+      call blosno_sublimsedim(ngrid, nlay, dtime, ustar, temp, q, qb, pplay, paprs, &
+                                     dtemp_bs, dq_bs, dqb_bs, bsfl, precip_bs)
+
+!=================================================================
+! Add tendencies
+!=================================================================
+      du0(:, :) = 0.
+      dv0(:, :) = 0.
+      dql0(:, :) = 0.
+      dqi0(:, :) = 0.
+
+      CALL add_phys_tend(du0, dv0, dtemp_bs, dq_bs, dql0, dqi0, dqb_bs, paprs, &
+                         'bsss', abortphy, flag_inhib_tend, itap, 0)
+
+   end subroutine call_blosno
+
+end module lmdz_call_blosno
Index: LMDZ6/trunk/libf/phylmd/lmdz_call_blowing_snow.f90
===================================================================
--- LMDZ6/trunk/libf/phylmd/lmdz_call_blowing_snow.f90	(revision 6027)
+++ 	(revision )
@@ -1,69 +1,0 @@
-!$gpum horizontal klon ngrid
-module lmdz_call_blowing_snow
-
-contains
-
-   subroutine call_blowing_snow(ngrid, nlay, abortphy, flag_inhib_tend, itap, dtime, ustar, temp, q, qb, pplay, paprs, &
-                                dtemp_bs, dq_bs, dqb_bs, bsfl, precip_bs)
-
-!==================================================================================================
-! call_blowing_snow in the main interface between the LMDZ physics monitor physiq_mod
-! and the blowing snow parameterizations
-!
-! contact: Etienne Vignon etienne.vignon@lmd.ipsl.fr
-!================================================================================================
-
-      use lmdz_blowing_snow_sublim_sedim, only: blowing_snow_sublim_sedim
-      USE add_phys_tend_mod, ONLY: add_phys_tend, prt_enerbil
-
-      implicit none
-!======================================================
-! Declarations
-!======================================================
-
-!INPUT
-!=====
-
-      INTEGER, INTENT(IN)                     :: abortphy, flag_inhib_tend, itap ! flag and physics time counter for add_phys_tend
-      integer, intent(in)                     :: ngrid ! number of horizontal grid points
-      integer, intent(in)                     :: nlay  ! number of vertical layers
-      real, intent(in)                        :: dtime ! physics time step [s]
-      real, intent(in), dimension(ngrid)      :: ustar ! surface friction velocity [m/s]
-      real, intent(in), dimension(ngrid, nlay) :: temp  ! temperature of the air [K]
-      real, intent(in), dimension(ngrid, nlay) :: q     ! specific content of water [kg/kg]
-      real, intent(in), dimension(ngrid, nlay) :: qb    ! specific content of blowing snow [kg/kg]
-      real, intent(in), dimension(ngrid, nlay) :: pplay ! pressure at the middle of layers [Pa]
-      real, intent(in), dimension(ngrid, nlay + 1) :: paprs ! pressure at the layer bottom interface [Pa]
-
-! OUTPUT
-!========
-      real, intent(out), dimension(ngrid, nlay) :: dtemp_bs ! temperature tendency [K/s]
-      real, intent(out), dimension(ngrid, nlay) :: dq_bs    ! water vapor tendency [kg/kg/s]
-      real, intent(out), dimension(ngrid, nlay) :: dqb_bs   ! blowing snow mass tendancy [kg/kg/s]
-      real, intent(out), dimension(ngrid, nlay + 1) :: bsfl   ! vertical profile of blowing snow vertical flux [kg/m2/s]
-      real, intent(out), dimension(ngrid)      :: precip_bs ! surface sedimentation flux of blowing snow [kg/s/m2]
-
-! LOCAL
-!========
-      real, dimension(ngrid, nlay)              :: du0, dv0, dql0, dqi0
-
-!=================================================================
-! Call to main routine of blowing snow sublimation and sedim.
-!=================================================================
-      call blowing_snow_sublim_sedim(ngrid, nlay, dtime, ustar, temp, q, qb, pplay, paprs, &
-                                     dtemp_bs, dq_bs, dqb_bs, bsfl, precip_bs)
-
-!=================================================================
-! Add tendencies
-!=================================================================
-      du0(:, :) = 0.
-      dv0(:, :) = 0.
-      dql0(:, :) = 0.
-      dqi0(:, :) = 0.
-
-      CALL add_phys_tend(du0, dv0, dtemp_bs, dq_bs, dql0, dqi0, dqb_bs, paprs, &
-                         'bsss', abortphy, flag_inhib_tend, itap, 0)
-
-   end subroutine call_blowing_snow
-
-end module lmdz_call_blowing_snow
Index: LMDZ6/trunk/libf/phylmd/lmdz_call_lscp.f90
===================================================================
--- LMDZ6/trunk/libf/phylmd/lmdz_call_lscp.f90	(revision 6027)
+++ LMDZ6/trunk/libf/phylmd/lmdz_call_lscp.f90	(revision 6028)
@@ -24,6 +24,6 @@
                         pfrac_impa, pfrac_nucl, pfrac_1nucl, &
                         prfl, psfl, rhcl, qta, fraca, &
-                        tv, pspsk, tla, thl, wth, iflag_cld_th, &
-                        iflag_ice_thermo, distcltop, temp_cltop, &
+                        tv, pspsk, tla, thl, wth, &
+                        distcltop, temp_cltop, &
                         tke, tke_dissip, &
                         fm_therm, entr_therm, detr_therm, &
@@ -66,5 +66,4 @@
       INTEGER, INTENT(IN)   :: nbsrf, is_ter, is_lic ! number of subgrid tiles and indices for land and landice
       INTEGER, INTENT(IN)   :: abortphy, flag_inhib_tend, itap ! flag and physics time counter for add_phys_tend
-      INTEGER, INTENT(IN)   :: iflag_ice_thermo! flag to activate the ice thermodynamics
       INTEGER, INTENT(IN)   :: nqo             ! number of water species
       LOGICAL, INTENT(IN)   :: ok_new_lscp     ! flag that controls the version of lscp code used
@@ -91,5 +90,4 @@
       !Inputs associated with thermal plumes
 
-      INTEGER, INTENT(IN)   :: iflag_cld_th    ! flag that determines the distribution of convective clouds
       REAL, DIMENSION(klon, klev), INTENT(IN)   :: tv              ! virtual potential temperature [K]
       REAL, DIMENSION(klon, klev), INTENT(IN)   :: qta             ! specific humidity within thermals [kg/kg]
@@ -250,6 +248,5 @@
       CALL ratqs_main_first(klon, cell_area)
 
-      CALL ratqs_main(klon, klev, nbsrf, is_ter, is_lic, &
-                      iflag_cld_th, dtime, &
+      CALL ratqs_main(klon, klev, nbsrf, is_ter, is_lic, dtime, &
                       pctsrf, s_pblh, zstd, wake_s, wake_deltaq, &
                       paprs, pplay, temp, qt, &
@@ -282,6 +279,6 @@
                    frac_impa, frac_nucl, beta, &
                    prfl, psfl, rhcl, &
-                   qta, fraca, tv, pspsk, tla, thl, wth, iflag_cld_th, &
-                   iflag_ice_thermo, distcltop, temp_cltop, &
+                   qta, fraca, tv, pspsk, tla, thl, wth, &
+                   distcltop, temp_cltop, &
                    tke, tke_dissip, &
                    entr_therm, detr_therm, &
@@ -310,6 +307,5 @@
                        frac_impa, frac_nucl, beta, &
                        prfl, psfl, rhcl, &
-                       qta, fraca, tv, pspsk, tla, thl, iflag_cld_th, &
-                       iflag_ice_thermo, &
+                       qta, fraca, tv, pspsk, tla, thl, & 
                        cloudth_sth, cloudth_senv, cloudth_sigmath, cloudth_sigmaenv)
 
Index: LMDZ6/trunk/libf/phylmd/lmdz_lscp_ini.f90
===================================================================
--- LMDZ6/trunk/libf/phylmd/lmdz_lscp_ini.f90	(revision 6027)
+++ LMDZ6/trunk/libf/phylmd/lmdz_lscp_ini.f90	(revision 6028)
@@ -34,4 +34,10 @@
    INTEGER, SAVE, PROTECTED :: niter_lscp = 5      ! number of iterations to calculate autoconversion to precipitation
    !$OMP THREADPRIVATE(niter_lscp)
+
+   INTEGER, SAVE, PROTECTED :: iflag_ice_thermo    ! activation of ice thermodynamics
+   !$OMP THREADPRIVATE(iflag_ice_thermo)
+
+   INTEGER, SAVE, PROTECTED :: iflag_cld_th        ! controls the properties of clouds in presence of thermals
+   !$OMP THREADPRIVATE(iflag_cld_th)
 
    INTEGER, SAVE, PROTECTED :: iflag_evap_prec = 1 ! precipitation evaporation flag. 0: nothing, 1: "old way",
@@ -398,6 +404,7 @@
 CONTAINS
 
-      SUBROUTINE lscp_ini(dtime, klon, klev, iflag_thermals, lunout_in, &
-                       prt_level_in, ok_ice_supersat_in, fl_cor_ebil_in, &
+      SUBROUTINE lscp_ini(dtime, klon, klev, iflag_thermals, iflag_cld_th_in,  &
+                       iflag_ice_thermo_in, lunout_in, prt_level_in,           &
+                       ok_ice_supersat_in, fl_cor_ebil_in, &
                        RCPD_in, RLSTT_in, RLVTT_in, RLMLT_in, RVTMP2_in, &
                        RTT_in, RD_in, RV_in, RG_in, RPI_in, EPS_W_in)
@@ -407,4 +414,5 @@
       REAL, INTENT(IN)      :: dtime
       INTEGER, INTENT(IN)   :: klon, klev, iflag_thermals
+      INTEGER, INTENT(IN)   :: iflag_cld_th_in, iflag_ice_thermo_in
       INTEGER, INTENT(IN)   :: lunout_in, prt_level_in, fl_cor_ebil_in
       LOGICAL, INTENT(IN)   :: ok_ice_supersat_in
@@ -421,4 +429,6 @@
       fl_cor_ebil = fl_cor_ebil_in
       ok_ice_supersat = ok_ice_supersat_in
+      iflag_cld_th=iflag_cld_th_in
+      iflag_ice_thermo=iflag_ice_thermo_in
 
       RG = RG_in
@@ -608,4 +618,7 @@
       WRITE (lunout, *) 'lscp_ini, temp_nowater', temp_nowater
       WRITE (lunout, *) 'lscp_ini, ok_bug_phase_lscp', ok_bug_phase_lscp
+      WRITE (lunout, *) 'lscp_ini, iflag_cld_th', iflag_cld_th
+      WRITE (lunout, *) 'lscp_ini, iflag_ice_thermo', iflag_ice_thermo
+
       ! for poprecip
       WRITE (lunout, *) 'lscp_ini, ok_poprecip', ok_poprecip
Index: LMDZ6/trunk/libf/phylmd/lmdz_lscp_main.f90
===================================================================
--- LMDZ6/trunk/libf/phylmd/lmdz_lscp_main.f90	(revision 6027)
+++ LMDZ6/trunk/libf/phylmd/lmdz_lscp_main.f90	(revision 6028)
@@ -18,6 +18,6 @@
                    frac_impa, frac_nucl, beta, &
                    prfl, psfl, rhcl, qta, fraca, &
-                   tv, pspsk, tla, thl, wth, iflag_cld_th, &
-                   iflag_ice_thermo, distcltop, temp_cltop, &
+                   tv, pspsk, tla, thl, wth, &
+                   distcltop, temp_cltop, &
                    tke, tke_dissip, &
                    entr_therm, detr_therm, &
@@ -83,4 +83,5 @@
       USE lmdz_lscp_ini, ONLY: ok_radocond_snow, a_tr_sca
       USE lmdz_lscp_ini, ONLY: iflag_cloudth_vert, iflag_t_glace, iflag_fisrtilp_qsat
+      USE lmdz_lscp_ini, ONLY: iflag_cld_th, iflag_ice_thermo
       USE lmdz_lscp_ini, ONLY: min_frac_th_cld, temp_nowater
       USE lmdz_lscp_ini, ONLY: RCPD, RLSTT, RLVTT, RVTMP2, RTT, RD, RG
@@ -99,5 +100,4 @@
 
       INTEGER, INTENT(IN)   :: klon, klev       ! number of horizontal grid points and vertical levels
-      INTEGER, INTENT(IN)   :: iflag_ice_thermo! flag to activate the ice thermodynamics
 
       REAL, INTENT(IN)   :: dtime           ! time step [s]
@@ -120,5 +120,4 @@
       !Inputs associated with thermal plumes
 
-      INTEGER, INTENT(IN)   :: iflag_cld_th    ! flag that determines the distribution of convective clouds
       REAL, DIMENSION(klon, klev), INTENT(IN)   :: tv                  ! virtual potential temperature [K]
       REAL, DIMENSION(klon, klev), INTENT(IN)   :: qta                 ! specific humidity within thermals [kg/kg]
Index: LMDZ6/trunk/libf/phylmd/lmdz_lscp_old.f90
===================================================================
--- LMDZ6/trunk/libf/phylmd/lmdz_lscp_old.f90	(revision 6027)
+++ LMDZ6/trunk/libf/phylmd/lmdz_lscp_old.f90	(revision 6028)
@@ -63,6 +63,5 @@
                        frac_impa, frac_nucl, beta, &
                        prfl, psfl, rhcl, zqta, fraca, &
-                       ztv, zpspsk, ztla, zthl, iflag_cld_th, &
-                       iflag_ice_thermo, &
+                       ztv, zpspsk, ztla, zthl, &
                        cloudth_sth, cloudth_senv, cloudth_sigmath, cloudth_sigmaenv)
 
@@ -75,4 +74,5 @@
       USE lmdz_lscp_ini, ONLY: fl_cor_ebil
       USE lmdz_lscp_ini, ONLY: iflag_t_glace, t_glace_min, t_glace_max, exposant_glace
+      USE lmdz_lscp_ini, ONLY: iflag_cld_th, iflag_ice_thermo
       USE lmdz_lscp_ini, ONLY: seuil_neb, rain_int_min, iflag_evap_prec, iflag_oldbug_fisrtilp, a_tr_sca
       USE lmdz_lscp_ini, ONLY: iflag_cloudth_vert, iflag_rain_incloud_vol
@@ -119,6 +119,4 @@
       REAL, DIMENSION(klon, klev), INTENT(IN)   :: q      ! humidite specifique (kg/kg)
       LOGICAL, DIMENSION(klon, klev), INTENT(IN)   :: ptconv ! points ou le schema de conv. prof. est actif
-      INTEGER, INTENT(IN)   :: iflag_cld_th
-      INTEGER, INTENT(IN)   :: iflag_ice_thermo
       !
       ! Inputs lies aux thermiques
Index: LMDZ6/trunk/libf/phylmd/lmdz_lscp_subgridvarq.f90
===================================================================
--- LMDZ6/trunk/libf/phylmd/lmdz_lscp_subgridvarq.f90	(revision 6027)
+++ LMDZ6/trunk/libf/phylmd/lmdz_lscp_subgridvarq.f90	(revision 6028)
@@ -31,6 +31,5 @@
 
 !=======================================================================
-   SUBROUTINE ratqs_main(klon, klev, nbsrf, is_ter, is_lic, &
-                         iflag_cld_th, pdtphys, &
+   SUBROUTINE ratqs_main(klon, klev, nbsrf, is_ter, is_lic, pdtphys, &
                          pctsrf, s_pblh, zstd, wake_s, wake_deltaq, &
                          paprs, pplay, t_seri, q_seri, &
@@ -40,5 +39,5 @@
                          ratqsc, ratqs, ratqs_inter_, sigma_qtherm)
 
-      USE lmdz_lscp_ini, ONLY: prt_level, lunout, iflag_ratqs
+      USE lmdz_lscp_ini, ONLY: prt_level, lunout, iflag_ratqs, iflag_cld_th
       USE lmdz_lscp_ini, ONLY: ratqsbas, ratqshaut
       USE lmdz_lscp_ini, ONLY: tau_ratqs, ratqsp0, ratqsdp
@@ -71,5 +70,4 @@
       integer, intent(in) :: klon, klev           ! horizontal and vertical dimensions
       integer, intent(in) :: nbsrf, is_ter, is_lic ! number of subgrid tiles and indices for land and landice
-      integer, intent(in) :: iflag_cld_th        ! flag that controls cloud properties in presence of thermals
       real, intent(in)     :: pdtphys             ! physics time step [s]
       real, dimension(klon, klev + 1), intent(in) :: paprs ! pressure at layer interfaces [Pa]
Index: LMDZ6/trunk/libf/phylmd/pbl_surface_subsrf_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/pbl_surface_subsrf_mod.F90	(revision 6027)
+++ LMDZ6/trunk/libf/phylmd/pbl_surface_subsrf_mod.F90	(revision 6028)
@@ -199,5 +199,5 @@
          frac_tersrf, z0m_tersrf, ratio_z0m_z0h_tersrf !AM 
     use phys_output_var_mod, only: tkt, tks, taur, sss
-    use lmdz_blowing_snow_ini, only : zeta_bs
+    use lmdz_blosno_ini, only : zeta_bs
     USE dimsoil_mod_h, ONLY: nsoilmx
     USE surf_param_mod, ONLY: eff_surf_param  !AM
Index: LMDZ6/trunk/libf/phylmd/physiq_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/physiq_mod.F90	(revision 6027)
+++ LMDZ6/trunk/libf/phylmd/physiq_mod.F90	(revision 6028)
@@ -95,5 +95,5 @@
     use wxios_mod, ONLY: g_ctx, wxios_set_context
     USE lmdz_call_cloud_optics_prop, ONLY : call_cloud_optics_prop, call_cloud_optics_prop_post
-    USE lmdz_call_blowing_snow, ONLY : call_blowing_snow
+    USE lmdz_call_blosno, ONLY : call_blosno
     USE lmdz_call_lscp, ONLY: call_lscp
     USE lmdz_call_reevap, ONLY : call_reevap
@@ -109,5 +109,5 @@
     USE calltherm_mod, ONLY : calltherm
     USE lmdz_thermcell_dtke, ONLY : thermcell_dtke
-    USE lmdz_blowing_snow_ini, ONLY : blowing_snow_ini , qbst_bs
+    USE lmdz_blosno_ini, ONLY : blosno_ini , qbst_bs
     USE lmdz_lscp_ini, ONLY : lscp_ini, ratqsbas
     USE lmdz_cloud_optics_prop_ini, ONLY : cloud_optics_prop_ini
@@ -1873,8 +1873,8 @@
        CALL thermcell_ini(iflag_thermals,prt_level,tau_thermals,lunout, &
             &    RG,RD,RCPD,RKAPPA,RLVTT,RETV)
-       CALL lscp_ini(pdtphys,klon,klev,iflag_thermals,lunout,prt_level,       &
-                     ok_ice_supersat,fl_cor_ebil,                             &
-                     RCPD,RLSTT,RLVTT,RLMLT,RVTMP2,RTT,RD,RV,RG,RPI,EPS_W)
-       CALL blowing_snow_ini(RCPD, RLSTT, RLVTT, RLMLT, &
+        CALL lscp_ini(pdtphys,klon,klev,iflag_thermals,iflag_cld_th,iflag_ice_thermo, & 
+                  lunout,prt_level,ok_ice_supersat,fl_cor_ebil,                       &
+                  RCPD,RLSTT,RLVTT,RLMLT,RVTMP2,RTT,RD,RV,RG,RPI,EPS_W)
+       CALL blosno_ini(RCPD, RLSTT, RLVTT, RLMLT, &
             RVTMP2, RTT,RD,RG, RV, RPI)
        ! Test de coherence sur oc_cdnc utilisé uniquement par cloud_optics_prop
@@ -3083,5 +3083,5 @@
     IF (ok_bs) THEN
 
-       CALL call_blowing_snow(klon,klev,abortphy,flag_inhib_tend,itap,phys_tstep,zustar, &
+       CALL call_blosno(klon,klev,abortphy,flag_inhib_tend,itap,phys_tstep,zustar, &
                               t_seri,q_seri,qbs_seri,pplay,paprs, &
                                d_t_bsss,d_q_bsss,d_qbs_bsss,bsfl,bs_fall)
@@ -3974,6 +3974,6 @@
      pfrac_impa, pfrac_nucl, pfrac_1nucl,               &
      prfl, psfl, rhcl, zqasc, fraca,                    &
-     ztv, zpspsk, ztla, zthl, zw2, iflag_cld_th,        &
-     iflag_ice_thermo, distcltop, temp_cltop,           &
+     ztv, zpspsk, ztla, zthl, zw2,                      &
+     distcltop, temp_cltop,                             &
      pbl_tke(:,:,is_ave), pbl_eps(:,:,is_ave),          &
      fm_therm, entr_therm, detr_therm,                  &
Index: LMDZ6/trunk/libf/phylmd/surf_landice_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/surf_landice_mod.F90	(revision 6027)
+++ LMDZ6/trunk/libf/phylmd/surf_landice_mod.F90	(revision 6028)
@@ -89,6 +89,6 @@
     USE clesphys_mod_h
     USE yomcst_mod_h
-    USE lmdz_blowing_snow_ini, ONLY : c_esalt_bs, zeta_bs, pbst_bs, prt_bs, rhoice_bs, rhohard_bs
-    USE lmdz_blowing_snow_ini, ONLY : rhofresh_bs, tau_eqsalt_bs, tau_dens0_bs, tau_densmin_bs
+    USE lmdz_blosno_ini, ONLY : c_esalt_bs, zeta_bs, pbst_bs, prt_bs, rhoice_bs, rhohard_bs
+    USE lmdz_blosno_ini, ONLY : rhofresh_bs, tau_eqsalt_bs, tau_dens0_bs, tau_densmin_bs
     USE surf_inlandsis_mod,  ONLY : surf_inlandsis
 
Index: LMDZ6/trunk/libf/phylmdiso/lmdz_blosno_ini.f90
===================================================================
--- LMDZ6/trunk/libf/phylmdiso/lmdz_blosno_ini.f90	(revision 6028)
+++ LMDZ6/trunk/libf/phylmdiso/lmdz_blosno_ini.f90	(revision 6028)
@@ -0,0 +1,1 @@
+link ../phylmd/lmdz_blosno_ini.f90
Index: LMDZ6/trunk/libf/phylmdiso/lmdz_blowing_snow_ini.f90
===================================================================
--- LMDZ6/trunk/libf/phylmdiso/lmdz_blowing_snow_ini.f90	(revision 6027)
+++ 	(revision )
@@ -1,1 +1,0 @@
-link ../phylmd/lmdz_blowing_snow_ini.f90
Index: LMDZ6/trunk/libf/phylmdiso/lmdz_blowing_snow_sublim_sedim.f90
===================================================================
--- LMDZ6/trunk/libf/phylmdiso/lmdz_blowing_snow_sublim_sedim.f90	(revision 6027)
+++ 	(revision )
@@ -1,1 +1,0 @@
-link ../phylmd/lmdz_blowing_snow_sublim_sedim.f90
Index: LMDZ6/trunk/libf/phylmdiso/pbl_surface_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmdiso/pbl_surface_mod.F90	(revision 6027)
+++ LMDZ6/trunk/libf/phylmdiso/pbl_surface_mod.F90	(revision 6028)
@@ -487,5 +487,5 @@
          frac_tersrf, z0m_tersrf, ratio_z0m_z0h_tersrf, albedo_tersrf !AM 
     USE phys_output_var_mod, ONLY : tkt, tks, taur, sss
-    USE lmdz_blowing_snow_ini, ONLY : zeta_bs
+    USE lmdz_blosno_ini, ONLY : zeta_bs
     USE wxios_mod, ONLY : missing_val_xios => missing_val, using_xios
     USE netcdf, ONLY : missing_val_netcdf => nf90_fill_real
Index: LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90	(revision 6027)
+++ LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90	(revision 6028)
@@ -97,5 +97,4 @@
     USE lmdz_call_cloud_optics_prop, ONLY : call_cloud_optics_prop, call_cloud_optics_prop_post
     USE lmdz_lscp_old, ONLY : fisrtilp, fisrtilp_first
-    !USE lmdz_call_blowing_snow, ONLY : call_blowing_snow
     USE calwake_mod, ONLY : calwake, calwake_first
     USE lmdz_wake_ini, ONLY : wake_ini
@@ -107,5 +106,5 @@
     USE calltherm_mod, ONLY : calltherm
     USE lmdz_thermcell_dtke, ONLY : thermcell_dtke
-    USE lmdz_blowing_snow_ini, ONLY : blowing_snow_ini , qbst_bs
+    USE lmdz_blosno_ini, ONLY : blosno_ini, qbst_bs
     USE lmdz_lscp_ini, ONLY : lscp_ini, ratqsbas
     USE lmdz_lscp_subgridvarq, ONLY : ratqs_main, ratqs_main_first
@@ -2024,9 +2023,10 @@
        CALL thermcell_ini(iflag_thermals,prt_level,tau_thermals,lunout, &
    &    RG,RD,RCPD,RKAPPA,RLVTT,RETV)
-       CALL lscp_ini(pdtphys,klon,klev,iflag_thermals,lunout,prt_level,       &
-                     ok_ice_supersat,fl_cor_ebil,                             &
-                     RCPD,RLSTT,RLVTT,RLMLT,RVTMP2,RTT,RD,RV,RG,RPI,EPS_W)
-       CALL blowing_snow_ini(RCPD, RLSTT, RLVTT, RLMLT, &
-                             RVTMP2, RTT,RD,RG, RV, RPI)
+       CALL lscp_ini(pdtphys,klon,klev,iflag_thermals,iflag_cld_th,iflag_ice_thermo, &
+                  lunout,prt_level,ok_ice_supersat,fl_cor_ebil,                       &
+                  RCPD,RLSTT,RLVTT,RLMLT,RVTMP2,RTT,RD,RV,RG,RPI,EPS_W)
+       CALL blosno_ini(RCPD, RLSTT, RLVTT, RLMLT, &
+            RVTMP2, RTT,RD,RG, RV, RPI)
+
        ! Test de coherence sur oc_cdnc utilisé uniquement par cloud_optics_prop
        IF (ok_newmicro) then
@@ -3617,5 +3617,5 @@
     !IF (ok_bs) THEN
 
-    ! CALL call_blowing_snow(klon,klev,abortphy,flag_inhib_tend,itap,phys_tstep,zustar, &
+    ! CALL call_blosno(klon,klev,abortphy,flag_inhib_tend,itap,phys_tstep,zustar, &
     !                          t_seri,q_seri,qbs_seri,pplay,paprs, &
     !                           d_t_bsss,d_q_bsss,d_qbs_bsss,bsfl,bs_fall)
Index: LMDZ6/trunk/libf/phylmdiso/surf_landice_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmdiso/surf_landice_mod.F90	(revision 6027)
+++ LMDZ6/trunk/libf/phylmdiso/surf_landice_mod.F90	(revision 6028)
@@ -51,6 +51,6 @@
     USE yomcst_mod_h
     USE ioipsl_getin_p_mod, ONLY : getin_p
-    USE lmdz_blowing_snow_ini, ONLY : c_esalt_bs, zeta_bs, pbst_bs, prt_bs, rhoice_bs, rhohard_bs
-    USE lmdz_blowing_snow_ini, ONLY : rhofresh_bs, tau_eqsalt_bs, tau_dens0_bs, tau_densmin_bs
+    USE lmdz_blosno_ini, ONLY : c_esalt_bs, zeta_bs, pbst_bs, prt_bs, rhoice_bs, rhohard_bs
+    USE lmdz_blosno_ini, ONLY : rhofresh_bs, tau_eqsalt_bs, tau_dens0_bs, tau_densmin_bs
     USE surf_inlandsis_mod,  ONLY : surf_inlandsis