Index: /LMDZ6/trunk/DefLists/field_def_lmdz.xml
===================================================================
--- /LMDZ6/trunk/DefLists/field_def_lmdz.xml	(revision 4834)
+++ /LMDZ6/trunk/DefLists/field_def_lmdz.xml	(revision 4835)
@@ -117,4 +117,5 @@
 	<field id="rhosnow_lic"    long_name="snow density lic"                unit="kg/m3" />
         <field id="ustart_lic"    long_name="ustar threshold for snow erosion"                unit="m/s" />
+	<field id="qsalt_lic"    long_name="blowing snow concentration in saltation layer"    unit="kg/kg" />
 	<field id="evap_ter"    long_name="evaporation at surface ter"    unit="kg/(s*m2)" />
         <field id="evap_lic"    long_name="evaporation at surface lic"    unit="kg/(s*m2)" />
Index: /LMDZ6/trunk/libf/phylmd/lmdz_blowing_snow_ini.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/lmdz_blowing_snow_ini.F90	(revision 4834)
+++ /LMDZ6/trunk/libf/phylmd/lmdz_blowing_snow_ini.F90	(revision 4835)
@@ -3,20 +3,28 @@
 implicit none
 
-   real, save, protected :: RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT,RD,RG
-   real, save, protected :: coef_eva_bs, expo_eva_bs, fallv_bs, zeta_bs
-   real, save, protected :: prt_bs, pbst_bs, qbst_bs
-   
-   integer, save, protected :: iflag_saltation_bs
+   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, RLSTT, RLVTT, RLMLT, RVTMP2, RTT,RD,RG)
-   !$OMP THREADPRIVATE(coef_eva_bs, expo_eva_bs, fallv_bs, zeta_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)
+   !$OMP THREADPRIVATE(iflag_sedim_bs, iflag_sublim_bs)
 
-   real, save, protected :: tbsmelt=278.15     ! parameter to calculate melting fraction of BS sedimentation
-   real, save, protected :: taumeltbs0=1800.0  ! Melting time scale of blowing snow at 273.15K
-   real, save, protected :: qbsmin=1.E-10      ! Minimum blowing snow specific content
+   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)
 
-   !$OMP THREADPRIVATE(tbsmelt, taumeltbs0, qbsmin)
+   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)
 
 
@@ -24,14 +32,15 @@
 
    subroutine blowing_snow_ini(RCPD_in, RLSTT_in, RLVTT_in, RLMLT_in,&
-                                  RVTMP2_in, RTT_in,RD_in,RG_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
-         real, intent(in) ::  RVTMP2_in, RTT_in, RD_in, RG_in
+         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
@@ -41,5 +50,8 @@
          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
@@ -58,13 +70,15 @@
          CALL getin_p('fallv_bs',fallv_bs)
 
-         coef_eva_bs =  2.e-5 
-         CALL getin_p('coef_eva_bs',coef_eva_bs)
+         coef_sub_bs =  0.1 
+         CALL getin_p('coef_sub_bs',coef_sub_bs)
 
-         expo_eva_bs = 0.5
-         CALL getin_p('expo_eva_bs',expo_eva_bs)
+         iflag_sublim_bs=1
+         CALL getin_p('iflag_sublim_bs',iflag_sublim_bs)
 
-         iflag_saltation_bs=0
-         CALL getin_p('iflag_saltation_bs',iflag_saltation_bs)
+         iflag_sedim_bs=1
+         CALL getin_p('iflag_sedim_bs',iflag_sedim_bs)
 
+         r_bs=150.0e-6
+         CALL getin_p('r_bs',r_bs)
 
       end subroutine 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 4834)
+++ /LMDZ6/trunk/libf/phylmd/lmdz_blowing_snow_sublim_sedim.F90	(revision 4835)
@@ -2,5 +2,5 @@
 
 contains
-subroutine blowing_snow_sublim_sedim(ngrid,nlay,dtime,temp,qv,qbs,pplay,paprs,dtemp_bs,dq_bs,dqbs_bs,bsfl,precip_bs)
+subroutine blowing_snow_sublim_sedim(ngrid,nlay,dtime,temp,qv,qb,pplay,paprs,dtemp_bs,dqv_bs,dqb_bs,bsfl,precip_bs)
 
 !==============================================================================
@@ -11,6 +11,6 @@
 
 
-use lmdz_blowing_snow_ini, only : coef_eva_bs,RTT,RD,RG,expo_eva_bs, fallv_bs, qbsmin
-use lmdz_blowing_snow_ini, only : RCPD, RLSTT, RLMLT, RLVTT, RVTMP2, tbsmelt, taumeltbs0
+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, tbsmelt, taumeltbs0, rhobs, r_bs
 USE lmdz_lscp_tools, only : calc_qsat_ecmwf
 
@@ -29,5 +29,5 @@
 real, intent(in), dimension(ngrid,nlay) :: temp
 real, intent(in), dimension(ngrid,nlay) :: qv
-real, intent(in), dimension(ngrid,nlay) :: qbs
+real, intent(in), dimension(ngrid,nlay) :: qb
 real, intent(in), dimension(ngrid,nlay) :: pplay
 real, intent(in), dimension(ngrid,nlay+1) :: paprs
@@ -40,6 +40,6 @@
 
 real, intent(out), dimension(ngrid,nlay) :: dtemp_bs
-real, intent(out), dimension(ngrid,nlay) :: dq_bs
-real, intent(out), dimension(ngrid,nlay) :: dqbs_bs
+real, intent(out), dimension(ngrid,nlay) :: dqv_bs
+real, intent(out), dimension(ngrid,nlay) :: dqb_bs
 real, intent(out), dimension(ngrid,nlay+1) :: bsfl
 real, intent(out), dimension(ngrid)      :: precip_bs
@@ -51,10 +51,12 @@
 
 integer                                  :: k,i,n
-real                                     :: cpd, cpw, dqsub, maxdqsub, dqbsmelt
-real                                     :: dqsedim,precbs, dqmelt, zmelt, taumeltbs
-real                                     :: maxdqmelt, rhoair, dz
-real, dimension(ngrid)                   :: zt,zq,zqbs,qsi,dqsi,qsl, dqsl,qzero,sedim,sedimn
-real, dimension(ngrid)                   :: zqbsi, zqbs_up, zmair
-real, dimension(ngrid)                   :: zvelo 
+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
+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
 
 !++++++++++++++++++++++++++++++++++++++++++++++++++
@@ -64,137 +66,229 @@
 qzero(:)=0.
 dtemp_bs(:,:)=0.
-dq_bs(:,:)=0.
-dqbs_bs(:,:)=0.
+dqv_bs(:,:)=0.
+dqb_bs(:,:)=0.
 zvelo(:)=0.
-zt(:)=0.
-zq(:)=0.
-zqbs(:)=0.
 sedim(:)=0.
-sedimn(:)=0.
-
+precip_bs(:)=0.
+bsfl(:,:)=0.
+
+
+Ka=2.4e-2      ! thermal conductivity of the air, SI
+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)
+   ENDDO
+ENDDO
+
+
+! Sedimentation scheme
+!----------------------
+
+IF (iflag_sedim_bs .GT. 0) THEN
 ! begin of top-down loop
 DO k = nlay, 1, -1
     
     DO i=1,ngrid
-        zt(i)=temp(i,k)
-        zq(i)=qv(i,k)
-        zqbs(i)=qbs(i,k)
+        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)
     ENDDO
 
-     ! thermalization of blowing snow precip coming from above     
+    ! thermalization of blowing snow precip coming from above     
     IF (k.LE.nlay-1) THEN
 
         DO i = 1, ngrid
-            zmair(i)=(paprs(i,k)-paprs(i,k+1))/RG
+            zmair=(zpaprsdn(i)-zpaprsup(i))/RG
             ! RVTMP2=rcpv/rcpd-1
-            cpd=RCPD*(1.0+RVTMP2*zq(i))
+            cpd=RCPD*(1.0+RVTMP2*zqv(i))
             cpw=RCPD*RVTMP2
-            ! zqbs_up: blowing snow mass that has to be thermalized with 
+            ! 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
-            zqbs_up(i) = sedim(i)*dtime/zmair(i)
+            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
-            zt(i) = ( (temp(i,k+1)+dtemp_bs(i,k+1))*zqbs_up(i)*cpw + cpd*zt(i) ) &
-             / (cpd + zqbs_up(i)*cpw)
+            ztemp(i) = ( ztemp_up(i)*zqb_up(i)*cpw + cpd*ztemp(i) ) &
+                  / (cpd + zqb_up(i)*cpw)
         ENDDO
-    ENDIF
-
-
-
-    ! calulation saturation specific humidity
-    CALL CALC_QSAT_ECMWF(ngrid,zt(:),qzero(:),pplay(:,k),RTT,2,.false.,qsi(:),dqsi(:))
-    CALL CALC_QSAT_ECMWF(ngrid,zt(:),qzero(:),pplay(:,k),RTT,1,.false.,qsl(:),dqsl(:))
-    
-    
-    ! 3 processes: melting, sublimation and precipitation of blowing snow
+
+     ENDIF
+
+     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) 
+
+     ENDDO ! Loop on ngrid
+
+
+ENDDO ! vertical loop
+
+
+!surface bs flux
+DO i = 1, ngrid
+  precip_bs(i) = sedim(i)
+ENDDO
+
+ENDIF
+
+
+
+
+!+++++++++++++++++++++++++++++++++++++++++++++++
+! 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)
+    ENDDO
+
+   ! calulation saturation specific humidity
+    CALL CALC_QSAT_ECMWF(ngrid,ztemp(:),qzero(:),zpres(:),RTT,2,.false.,qsi(:),dqsi(:))
+
+
     DO i = 1, ngrid
 
-          rhoair  = pplay(i,k) / zt(i) / RD
-          dz      = (paprs(i,k)-paprs(i,k+1)) / (rhoair*RG)
+          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 (zt(i) .GT. RTT) THEN
-             
+          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.,(zt(i)-RTT)/(tbsmelt-RTT)))
-             dqmelt=min(zqbs(i),-1.*zqbs(i)*(exp(-dtime/taumeltbs)-1.))
-             maxdqmelt= max(RCPD*(1.0+RVTMP2*(zq(i)+zqbs(i)))*(zt(i)-RTT)/(RLMLT+RLVTT),0.)
-             dqmelt=min(dqmelt,maxdqmelt)
-             ! q update, melting + evaporation
-             zq(i) = zq(i) + dqmelt
+
+             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
-             zt(i) = zt(i) - dqmelt * (RLMLT+RLVTT)/RCPD/(1.0+RVTMP2*(zq(i)+zqbs(i)))
-             ! qbs update melting + evaporation
-             zqbs(i)=zqbs(i)-dqmelt
-
-          ELSE 
+             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   
-             
-              rhoair=pplay(i,k)/zt(i)/RD 
-              dqsub = 1./rhoair*coef_eva_bs*(1.0-zq(i)/qsi(i))*((rhoair*zvelo(i)*zqbs(i))**expo_eva_bs)*dtime
-              dqsub = MAX(0.0,MIN(dqsub,zqbs(i)))
-
-              ! Sublimation limit: we ensure that the whole mesh does not exceed saturation wrt ice
-              maxdqsub = MAX(0.0, qsi(i)-zq(i))
-              dqsub = MIN(dqsub,maxdqsub)
-
+
+
+              ! 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*r_bs/(Aprime+Bprime)
+              ! assuming Mi=rhobs*4/3*pi*r_bs**3 
+              ! rhoair qb=nc*Mi -> nc=rhoair qb/Mi
+              ! dqb/dt_sub=-coef_sub_bs*(1-qv/qsi)*6*rhoair*qb/(rhobs*pi*r_bs**2)/(Aprime+Bprime)
+              ! dqb/dt_sub=-c_sub(1-qv/qsi)*qb
+              ! c_sub=coef_sub_bs*6*rhoair/(rhobs*pi*r_bs**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
+                 !diffusivity of water vapor
+                 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*6.*rhoair/(rhobs*RPI*r_bs**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.
+              ENDIF
 
               ! vapor, temperature, precip fluxes update following sublimation
-              zq(i) = zq(i) + dqsub     
-              zqbs(i)=zqbs(i)-dqsub           
-              zt(i) = zt(i) - dqsub*RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zqbs(i)))
+              zqv(i) = zqv(i) - dqbsub
+              zqb(i) = zqb(i) + dqbsub
+              ztemp(i) = ztemp(i) + dqbsub*RLSTT/RCPD/(1.0+RVTMP2*(zqv(i)+zqb(i)))
 
           ENDIF
 
-        ! Sedimentation scheme
-        !----------------------
-
-        sedimn(i) = rhoair*zqbs(i)*zvelo(i)
-
-        ! exact numerical resolution of sedimentation
-        ! assuming fall velocity is constant 
-
-        zqbs(i) = zqbs(i)*exp(-zvelo(i)/dz*dtime)+sedim(i)/rhoair/zvelo(i)
-
-        ! flux update
-        sedim(i) = sedimn(i)
-
-
-
-        ! if qbs<qbsmin, sublimate or melt and evaporate qbs
-        ! see Gerber et al. 2023, JGR Atmos for the choice of qbsmin
-        IF (zqbs(i) .LT. qbsmin) THEN                       
-              zq(i) = zq(i)+zqbs(i)
-              IF (zt(i) .LT. RTT) THEN
-                 zt(i) = zt(i) - zqbs(i) * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zqbs(i)))
+          ! 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
-                 zt(i) = zt(i) - zqbs(i) * (RLVTT+RLMLT)/RCPD/(1.0+RVTMP2*(zq(i)+zqbs(i)))
+                 ztemp(i) = ztemp(i) - zqb(i) * (RLVTT+RLMLT)/RCPD/(1.0+RVTMP2*(zqv(i)))
               ENDIF
-              zqbs(i)=0.
-        ENDIF
-
-
-    ! Outputs:
-        bsfl(i,k)=sedim(i)
-        dqbs_bs(i,k) = zqbs(i)-qbs(i,k)
-        dq_bs(i,k) = zq(i) - qv(i,k)
-        dtemp_bs(i,k) = zt(i) - temp(i,k)
-
-    ENDDO ! Loop on ngrid
-
-
-ENDDO ! vertical loop
-
-
-!surface bs flux
-DO i = 1, ngrid
-  precip_bs(i) = sedim(i)
-ENDDO
+              zqb(i)=0.
+          ENDIF
+
+     ! variables update
+     temp_seri(i,k)=ztemp(i)
+     qv_seri(i,k)=zqv(i)
+     qb_seri(i,k)=zqb(i)
+     ENDDO
+ENDDO
+
+ENDIF
+
+
+! 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)
+   ENDDO
+ENDDO
+
 
 
Index: /LMDZ6/trunk/libf/phylmd/phys_local_var_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/phys_local_var_mod.F90	(revision 4834)
+++ /LMDZ6/trunk/libf/phylmd/phys_local_var_mod.F90	(revision 4835)
@@ -335,6 +335,6 @@
       REAL,ALLOCATABLE,SAVE,DIMENSION(:) :: tpot, tpote, ue, uq, uwat, ve, vq, vwat, zxffonte
 !$OMP THREADPRIVATE(tpot, tpote, ue, uq, uwat, ve, vq, vwat, zxffonte)
-      REAL,ALLOCATABLE,SAVE,DIMENSION(:) :: zxustartlic, zxrhoslic
-!$OMP THREADPRIVATE(zxustartlic, zxrhoslic)
+      REAL,ALLOCATABLE,SAVE,DIMENSION(:) :: zxustartlic, zxrhoslic, zxqsaltlic
+!$OMP THREADPRIVATE(zxustartlic, zxrhoslic, zxqsaltlic)
       REAL,ALLOCATABLE,SAVE,DIMENSION(:) :: zxfqcalving
 !$OMP THREADPRIVATE(zxfqcalving)
@@ -842,6 +842,6 @@
       ALLOCATE(zxtsol(klon), snow_lsc(klon), zxfqfonte(klon), zxqsurf(klon))
       ALLOCATE(zxrunofflic(klon))
-      ALLOCATE(zxustartlic(klon), zxrhoslic(klon))
-      zxustartlic(:)=0. ; zxrhoslic(:)=0.
+      ALLOCATE(zxustartlic(klon), zxrhoslic(klon), zxqsaltlic(klon))
+      zxustartlic(:)=0. ; zxrhoslic(:)=0. ; zxqsaltlic(:)=0.
       ALLOCATE(rain_lsc(klon))
       ALLOCATE(rain_num(klon))
@@ -1174,5 +1174,5 @@
       DEALLOCATE(zxfqcalving, zxfluxlat)
       DEALLOCATE(zxrunofflic)
-      DEALLOCATE(zxustartlic, zxrhoslic)
+      DEALLOCATE(zxustartlic, zxrhoslic, zxqsaltlic)
       DEALLOCATE(zxtsol, snow_lsc, zxfqfonte, zxqsurf)
       DEALLOCATE(rain_lsc)
Index: /LMDZ6/trunk/libf/phylmd/phys_output_ctrlout_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/phys_output_ctrlout_mod.F90	(revision 4834)
+++ /LMDZ6/trunk/libf/phylmd/phys_output_ctrlout_mod.F90	(revision 4835)
@@ -388,4 +388,6 @@
   TYPE(ctrl_out), SAVE :: o_rhosnow_lic = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     'rhosnow_lic', 'snow density lic', 'kg/m3', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_qsalt_lic = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'qsalt_lic', 'qb in saltation layer lic', 'kg/kg', (/ ('', i=1, 10) /))
   TYPE(ctrl_out), SAVE :: o_sens_prec_liq_oce = ctrl_out((/ 5, 5, 10, 10, 5, 10, 11, 11, 11, 11/), &
     'sens_rain_oce', 'Sensible heat flux of liquid prec. over ocean', 'W/m2', (/ ('', i=1, 10) /))
Index: /LMDZ6/trunk/libf/phylmd/phys_output_write_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/phys_output_write_mod.F90	(revision 4834)
+++ /LMDZ6/trunk/libf/phylmd/phys_output_write_mod.F90	(revision 4835)
@@ -48,5 +48,5 @@
          o_psol, o_mass, o_qsurf, o_qsol, &
          o_precip, o_rain_fall, o_rain_con, o_ndayrain, o_plul, o_pluc, o_plun, &
-         o_snow, o_msnow, o_fsnow, o_evap, o_snowerosion, o_ustart_lic, o_rhosnow_lic, o_bsfall, & 
+         o_snow, o_msnow, o_fsnow, o_evap, o_snowerosion, o_ustart_lic, o_qsalt_lic, o_rhosnow_lic, o_bsfall, & 
          o_ep,o_epmax_diag, & ! epmax_cape
          o_tops, o_tops0, o_topl, o_topl0, &
@@ -310,5 +310,5 @@
     USE phys_local_var_mod, ONLY: zxfluxlat, slp, ptstar, pt0, zxtsol, zt2m, &
          zn2mout, t2m_min_mon, t2m_max_mon, evap, &
-         snowerosion, zxustartlic, zxrhoslic, &
+         snowerosion, zxustartlic, zxrhoslic, zxqsaltlic, &
          l_mixmin,l_mix, tke_dissip, &
          zu10m, zv10m, zq2m, zustar, zxqsurf, &
@@ -897,4 +897,5 @@
            CALL histwrite_phy(o_ustart_lic, zxustartlic)
            CALL histwrite_phy(o_rhosnow_lic, zxrhoslic)
+           CALL histwrite_phy(o_qsalt_lic, zxqsaltlic)
        ENDIF
 
Index: /LMDZ6/trunk/libf/phylmd/physiq_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/physiq_mod.F90	(revision 4834)
+++ /LMDZ6/trunk/libf/phylmd/physiq_mod.F90	(revision 4835)
@@ -1864,5 +1864,5 @@
        CALL lscp_ini(pdtphys,lunout,prt_level,ok_ice_sursat,iflag_ratqs,fl_cor_ebil,RCPD,RLSTT,RLVTT,RLMLT,RVTMP2,RTT,RD,RG,RPI)
        CALL blowing_snow_ini(RCPD, RLSTT, RLVTT, RLMLT, &
-                             RVTMP2, RTT,RD,RG)
+                             RVTMP2, RTT,RD,RG, RV, RPI)
        ! Test de coherence sur oc_cdnc utilisé uniquement par cloud_optics_prop
        IF (ok_newmicro) then
Index: /LMDZ6/trunk/libf/phylmd/surf_landice_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/surf_landice_mod.F90	(revision 4834)
+++ /LMDZ6/trunk/libf/phylmd/surf_landice_mod.F90	(revision 4835)
@@ -31,10 +31,10 @@
     USE cpl_mod,          ONLY : cpl_send_landice_fields
     USE calcul_fluxs_mod
-    USE phys_local_var_mod, ONLY : zxrhoslic, zxustartlic
+    USE phys_local_var_mod, ONLY : zxrhoslic, zxustartlic, zxqsaltlic
     USE phys_output_var_mod, ONLY : snow_o,zfra_o
 !FC
     USE ioipsl_getin_p_mod, ONLY : getin_p
-    USE lmdz_blowing_snow_ini, ONLY : zeta_bs, pbst_bs, prt_bs, iflag_saltation_bs
-
+    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
 #ifdef CPP_INLANDSIS
     USE surf_inlandsis_mod,  ONLY : surf_inlandsis
@@ -140,9 +140,8 @@
     REAL,DIMENSION(klon) :: precip_totsnow, evap_totsnow
     REAL, DIMENSION (klon,6) :: alb6
-    REAL                   :: rho0, rhoice, ustart0,esalt, rhod
+    REAL                   :: esalt
     REAL                   :: lambdasalt,fluxsalt, csalt, nunu, aa, bb, cc
-    REAL                   :: tau_dens, tau_dens0, tau_densmin, rhomax, rhohard
-    REAL, DIMENSION(klon)  :: ws1, rhos, ustart, qsalt, hsalt, snowini
-    REAL                   :: ustarmin, maxerosion, tau_eq_salt 
+    REAL                   :: tau_dens, maxerosion, fluxbs_2
+    REAL, DIMENSION(klon)  :: ws1, rhod, rhos, ustart0, ustart, qsalt, hsalt
 
 ! End definition
@@ -331,12 +330,5 @@
 !
 !****************************************************************************************
- 
-    CALL albsno(klon,knon,dtime,agesno(:),alb_neig(:), precip_snow(:))  
-
-
-! EV: following lines are obsolete since we set alb1 and alb2 to constant values
-! I therefore comment them
-!    alb1(1:knon) = alb_neig(1:knon)*zfra(1:knon) + &
-!         0.6 * (1.0-zfra(1:knon))
+
 !
 !IM: plusieurs choix/tests sur l'albedo des "glaciers continentaux"
@@ -361,18 +353,55 @@
 
 
-
-  ! Simple blowing snow param
-  if (ok_bs) then
-       ustarmin=1.e-3
-       ustart0 = 0.211
-       rhoice = 920.0
-       rho0 = 300.0
-       rhomax=450.0
-       rhohard=450.0
-       tau_dens0=86400.0*10.  ! 10 days by default, in s
-       tau_densmin=21600.0    ! 1/4 days according to in situ obs by C. Amory during blowing snow + 
-                              ! Marshall et al. 1999 (snow densification during rain)
-       tau_eq_salt=10.
-
+!****************************************************************************************
+! Simple blowing snow param
+!****************************************************************************************
+! we proceed in 2 steps:
+! first we erode - if possible -the accumulated snow during the time step
+! then we update the density of the underlying layer and see if we can also erode
+! this layer
+
+
+   if (ok_bs) then
+       fluxbs(:)=0.
+       do j=1,klon
+          ws1(j)=(u1(j)**2+v1(j)**2)**0.5
+          ustar(j)=(cdragm(j)*(u1(j)**2+v1(j)**2))**0.5
+          rhod(j)=p1lay(j)/RD/temp_air(j)
+          ustart0(j) =(log(2.868)-log(1.625))/0.085*sqrt(cdragm(j))
+       enddo
+
+       ! 1st step: erosion of fresh snow accumulated during the time step
+       do j=1, knon
+
+           rhos(j)=rhofresh_bs
+           ! blowing snow flux formula used in MAR
+           ustart(j)=ustart0(j)*exp(max(rhoice_bs/rhofresh_bs-rhoice_bs/rhos(j),0.))*exp(max(0.,rhos(j)-rhohard_bs))
+           ! we have multiplied by exp to prevent erosion when rhos>rhohard_bs 
+           ! computation of qbs at the top of the saltation layer
+           ! default formulation from MAR model (Amory et al. 2021, Gallee et al. 2001)
+           esalt=1./(c_esalt_bs*max(1.e-6,ustar(j)))
+           hsalt(j)=0.08436*(max(1.e-6,ustar(j))**1.27)
+           qsalt(j)=(max(ustar(j)**2-ustart(j)**2,0.))/(RG*hsalt(j))*esalt
+           ! calculation of erosion (flux positive towards the surface here) 
+           ! consistent with implicit resolution of turbulent mixing equation
+           ! Nemoto and Nishimura 2004 show that steady-state saltation is achieved within a time tau_eqsalt_bs of about 10s
+           ! we thus prevent snowerosion (snow particle transfer from the saltation layer to the first model level)
+           ! integrated over tau_eqsalt_bs to exceed the total mass of snow particle in the saltation layer
+           ! (rho*qsalt*hsalt)
+           ! during this first step we also lower bound the erosion to the amount of fresh snow accumulated during the time step
+           maxerosion=min(precip_snow(j),hsalt(j)*qsalt(j)*rhod(j)/tau_eqsalt_bs)
+
+           fluxbs(j)=rhod(j)*ws1(j)*cdragm(j)*(AcoefQBS(j)-qsalt(j)) &
+                   / (1.-rhod(j)*ws1(j)*cdragm(j)*BcoefQBS(j)*dtime)
+           !fluxbs(j)= rhod*ws1(j)*cdragm(j)*(qbs1(j)-qsalt(j))
+           fluxbs(j)=max(-maxerosion,fluxbs(j))
+       enddo
+
+
+       ! we now compute the snow age of the overlying layer (snow surface after erosion of the fresh snow accumulated during the time step)
+       ! this is done through the routine albsno
+       CALL albsno(klon,knon,dtime,agesno(:),alb_neig(:), precip_snow(:)+fluxbs(:))
+
+       ! 2nd step:
        ! computation of threshold friction velocity
        ! which depends on surface snow density
@@ -381,70 +410,43 @@
            ! snow density increases with snow age and
            ! increases even faster in case of sedimentation of blowing snow or rain
-           tau_dens=max(tau_densmin, tau_dens0*exp(-abs(precip_bs(j))/pbst_bs - & 
-                    abs(precip_rain(j))/prt_bs)*exp(-max(tsurf(j)-272.0,0.)))
-           rhos(j)=rho0+(rhohard-rho0)*(1.-exp(-agesno(j)*86400.0/tau_dens))
+           tau_dens=max(tau_densmin_bs, tau_dens0_bs*exp(-abs(precip_bs(j))/pbst_bs - & 
+                    abs(precip_rain(j))/prt_bs)*exp(-max(tsurf(j)-RTT,0.)))
+           rhos(j)=rhofresh_bs+(rhohard_bs-rhofresh_bs)*(1.-exp(-agesno(j)*86400.0/tau_dens))
            ! blowing snow flux formula used in MAR
-           ws1(j)=(u1(j)**2+v1(j)**2)**0.5
-           ustar(j)=(cdragm(j)*(u1(j)**2+v1(j)**2))**0.5
-           ustart(j)=ustart0*exp(max(rhoice/rho0-rhoice/rhos(j),0.))*exp(max(0.,rhos(j)-rhomax)) 
-           ! we have multiplied by exp to prevent erosion when rhos>rhomax (usefull till
-           ! rhohard<450)
+           ustart(j)=ustart0(j)*exp(max(rhoice_bs/rhofresh_bs-rhoice_bs/rhos(j),0.))*exp(max(0.,rhos(j)-rhohard_bs)) 
+           ! we have multiplied by exp to prevent erosion when rhos>rhohard_bs 
+           ! computation of qbs at the top of the saltation layer
+           ! default formulation from MAR model (Amory et al. 2021, Gallee et al. 2001)
+           esalt=c_esalt_bs*ustar(j)
+           hsalt(j)=0.08436*(max(1.e-6,ustar(j))**1.27)
+           qsalt(j)=(max(ustar(j)**2-ustart(j)**2,0.))/(RG*hsalt(j))*esalt
+           ! calculation of erosion (flux positive towards the surface here) 
+           ! consistent with implicit resolution of turbulent mixing equation
+           ! Nemoto and Nishimura 2004 show that steady-state saltation is achieved within a time tau_eqsalt_bs of about 10s
+           ! we thus prevent snowerosion (snow particle transfer from the saltation layer to the first model level)
+           ! integrated over tau_eqsalt_bs to exceed the total mass of snow particle in the saltation layer
+           ! (rho*qsalt*hsalt)
+           maxerosion=hsalt(j)*qsalt(j)*rhod(j)/tau_eqsalt_bs
+           fluxbs_2=rhod(j)*ws1(j)*cdragm(j)*(AcoefQBS(j)-qsalt(j)) &
+                   / (1.-rhod(j)*ws1(j)*cdragm(j)*BcoefQBS(j)*dtime)
+           !fluxbs_2= rhod*ws1(j)*cdragm(j)*(qbs1(j)-qsalt(j))
+           fluxbs_2=max(-maxerosion,fluxbs_2)
+           fluxbs(j)=fluxbs(j)+fluxbs_2
        enddo
-       
-       ! computation of qbs at the top of the saltation layer
-       ! two formulations possible
-       ! pay attention that qbs is a mixing ratio and has to be converted
-       ! to specific content
-       
-       if (iflag_saltation_bs .eq. 1) then
-       ! expression from CRYOWRF (Sharma et al. 2022)
-          aa=2.6
-          bb=2.5
-          cc=2.0
-          lambdasalt=0.45
-          do j =1, knon
-               rhod=p1lay(j)/RD/temp_air(j)
-               nunu=max(ustar(j)/ustart(j),1.e-3)
-               fluxsalt=rhod/RG*(ustar(j)**3)*(1.-nunu**(-2)) * &
-                        (aa+bb*nunu**(-2)+cc*nunu**(-1))  
-               csalt=fluxsalt/(2.8*ustart(j))
-               hsalt(j)=0.08436*ustar(j)**1.27
-               qsalt(j)=1./rhod*csalt*lambdasalt*RG/(max(ustar(j)**2,ustarmin)) &
-                       * exp(-lambdasalt*RG*hsalt(j)/max(ustar(j)**2,ustarmin))
-               qsalt(j)=max(qsalt(j),0.)
-          enddo
-
-
-       else
-       ! default formulation from MAR model (Amory et al. 2021, Gallee et al. 2001)        
-          do j=1, knon
-              esalt=1./(3.25*max(ustarmin,ustar(j)))
-              hsalt(j)=0.08436*(max(ustarmin,ustar(j))**1.27)
-              qsalt(j)=(max(ustar(j)**2-ustart(j)**2,0.))/(RG*hsalt(j))*esalt
-              !ep=qsalt*cdragm(j)*sqrt(u1(j)**2+v1(j)**2)
-          enddo
-       endif
-
-        ! calculation of erosion (flux positive towards the surface here) 
-        ! consistent with implicit resolution of turbulent mixing equation
-       do j=1, knon
+
+       ! for outputs       
+        do j=1, knon
               i = knindex(j)
-              rhod=p1lay(j)/RD/temp_air(j)
-              ! Nemoto and Nishimura 2004 show that steady-state saltation is achieved within a time tau_eq_salt of about 10s
-              ! we thus prevent snowerosion (snow particle transfer from the saltation layer to the first model level)
-              ! integrated over tau_eq_salt to exceed the total mass of snow particle in the saltation layer
-              ! (rho*qsalt*hsalt)
-              maxerosion=hsalt(j)*qsalt(j)*rhod/tau_eq_salt
-              fluxbs(j)=rhod*ws1(j)*cdragm(j)*(AcoefQBS(j)-qsalt(j)) &
-                       / (1.-rhod*ws1(j)*cdragm(j)*BcoefQBS(j)*dtime)
-              !fluxbs(j)= rhod*ws1(j)*cdragm(j)*(qbs1(j)-qsalt(j))
-              fluxbs(j)=min(0.,max(-maxerosion,fluxbs(j)))
-              ! for outputs
-
               zxustartlic(i) = ustart(j)
               zxrhoslic(i) = rhos(j)
-       enddo
-
-  endif
+              zxqsaltlic(i)=qsalt(j)
+        enddo
+
+
+  else
+  ! those lines are useful to calculate the snow age
+       CALL albsno(klon,knon,dtime,agesno(:),alb_neig(:), precip_snow(:))
+
+  endif ! if ok_bs
 
 
Index: /LMDZ6/trunk/libf/phylmdiso/phys_local_var_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmdiso/phys_local_var_mod.F90	(revision 4834)
+++ /LMDZ6/trunk/libf/phylmdiso/phys_local_var_mod.F90	(revision 4835)
@@ -363,6 +363,6 @@
       REAL,ALLOCATABLE,SAVE,DIMENSION(:) :: tpot, tpote, ue, uq, uwat, ve, vq, vwat, zxffonte
 !$OMP THREADPRIVATE(tpot, tpote, ue, uq, uwat, ve, vq, vwat, zxffonte)
-      REAL,ALLOCATABLE,SAVE,DIMENSION(:) :: zxustartlic, zxrhoslic
-!$OMP THREADPRIVATE(zxustartlic, zxrhoslic)
+      REAL,ALLOCATABLE,SAVE,DIMENSION(:) :: zxustartlic, zxrhoslic, zxqsaltlic
+!$OMP THREADPRIVATE(zxustartlic, zxrhoslic, zxqsaltlic)
       REAL,ALLOCATABLE,SAVE,DIMENSION(:) :: zxfqcalving
 !$OMP THREADPRIVATE(zxfqcalving)
@@ -958,4 +958,6 @@
       ALLOCATE(zxtsol(klon), snow_lsc(klon), zxfqfonte(klon), zxqsurf(klon))
       ALLOCATE(zxrunofflic(klon))
+      ALLOCATE(zxustartlic(klon), zxrhoslic(klon), zxqsaltlic(klon))
+      zxustartlic(:)=0. ; zxrhoslic(:)=0. ; zxqsaltlic(:)=0.
       ALLOCATE(rain_lsc(klon))
       ALLOCATE(rain_num(klon))
Index: /LMDZ6/trunk/libf/phylmdiso/phys_output_ctrlout_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmdiso/phys_output_ctrlout_mod.F90	(revision 4834)
+++ /LMDZ6/trunk/libf/phylmdiso/phys_output_ctrlout_mod.F90	(revision 4835)
@@ -388,4 +388,6 @@
   TYPE(ctrl_out), SAVE :: o_rhosnow_lic = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     'rhosnow_lic', 'snow density lic', 'kg/m3', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_qsalt_lic = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'qsalt_lic', 'qb in saltation layer lic', 'kg/kg', (/ ('', i=1, 10) /)) 
   TYPE(ctrl_out), SAVE :: o_sens_prec_liq_oce = ctrl_out((/ 5, 5, 10, 10, 5, 10, 11, 11, 11, 11/), &
     'sens_rain_oce', 'Sensible heat flux of liquid prec. over ocean', 'W/m2', (/ ('', i=1, 10) /))
Index: /LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90	(revision 4834)
+++ /LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90	(revision 4835)
@@ -1952,5 +1952,5 @@
        CALL lscp_ini(pdtphys,lunout,prt_level,ok_ice_sursat,iflag_ratqs,fl_cor_ebil,RCPD,RLSTT,RLVTT,RLMLT,RVTMP2,RTT,RD,RG,RPI)
        CALL blowing_snow_ini(RCPD, RLSTT, RLVTT, RLMLT, &
-                             RVTMP2, RTT,RD,RG)
+                             RVTMP2, RTT,RD,RG, RV, RPI)
        ! Test de coherence sur oc_cdnc utilisé uniquement par cloud_optics_prop
        IF (ok_newmicro) then
Index: /LMDZ6/trunk/libf/phylmdiso/surf_landice_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmdiso/surf_landice_mod.F90	(revision 4834)
+++ /LMDZ6/trunk/libf/phylmdiso/surf_landice_mod.F90	(revision 4835)
@@ -35,5 +35,5 @@
     USE cpl_mod,          ONLY : cpl_send_landice_fields
     USE calcul_fluxs_mod
-    USE phys_local_var_mod, ONLY : zxrhoslic, zxustartlic
+    USE phys_local_var_mod, ONLY : zxrhoslic, zxustartlic, zxqsaltlic
     USE phys_output_var_mod, ONLY : snow_o,zfra_o
 #ifdef ISO   
@@ -50,6 +50,6 @@
 !FC
     USE ioipsl_getin_p_mod, ONLY : getin_p
-    USE lmdz_blowing_snow_ini, ONLY : zeta_bs, pbst_bs, prt_bs, iflag_saltation_bs
-
+    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
 #ifdef CPP_INLANDSIS
     USE surf_inlandsis_mod,  ONLY : surf_inlandsis
@@ -57,5 +57,4 @@
 
     USE indice_sol_mod
-
 
 !    INCLUDE "indicesol.h"
@@ -185,9 +184,8 @@
     REAL,DIMENSION(klon) :: precip_totsnow, evap_totsnow
     REAL, DIMENSION (klon,6) :: alb6
-    REAL                   :: rho0, rhoice, ustart0,esalt, rhod
+    REAL                   :: esalt
     REAL                   :: lambdasalt,fluxsalt, csalt, nunu, aa, bb, cc
-    REAL                   :: tau_dens, tau_dens0, tau_densmin, rhomax, rhohard
-    REAL, DIMENSION(klon)  :: ws1, rhos, ustart, qsalt, hsalt, snowini
-    REAL                   :: maxerosion ! in kg/s
+    REAL                   :: tau_dens, maxerosion, fluxbs_2
+    REAL, DIMENSION(klon)  :: ws1, rhod, rhos, ustart0, ustart, qsalt, hsalt
 
 ! End definition
@@ -420,12 +418,5 @@
 !
 !****************************************************************************************
- 
-    CALL albsno(klon,knon,dtime,agesno(:),alb_neig(:), precip_snow(:))  
-
-
-! EV: following lines are obsolete since we set alb1 and alb2 to constant values
-! I therefore comment them
-!    alb1(1:knon) = alb_neig(1:knon)*zfra(1:knon) + &
-!         0.6 * (1.0-zfra(1:knon))
+
 !
 !IM: plusieurs choix/tests sur l'albedo des "glaciers continentaux"
@@ -450,16 +441,55 @@
 
 
-
-  ! Simple blowing snow param
-  if (ok_bs) then
-       ustart0 = 0.211
-       rhoice = 920.0
-       rho0 = 300.0
-       rhomax=450.0
-       rhohard=450.0
-       tau_dens0=86400.0*10.  ! 10 days by default, in s
-       tau_densmin=21600.0    ! 1/4 days according to in situ obs by C. Amory during blowing snow +
-                              ! Marshall et al. 1999 (snow densification during rain)
-
+!****************************************************************************************
+! Simple blowing snow param
+!****************************************************************************************
+! we proceed in 2 steps:
+! first we erode - if possible -the accumulated snow during the time step
+! then we update the density of the underlying layer and see if we can also erode
+! this layer
+
+
+   if (ok_bs) then
+       fluxbs(:)=0.
+       do j=1,klon
+          ws1(j)=(u1(j)**2+v1(j)**2)**0.5
+          ustar(j)=(cdragm(j)*(u1(j)**2+v1(j)**2))**0.5
+          rhod(j)=p1lay(j)/RD/temp_air(j)
+          ustart0(j) =(log(2.868)-log(1.625))/0.085*sqrt(cdragm(j))
+       enddo
+
+       ! 1st step: erosion of fresh snow accumulated during the time step
+       do j=1, knon
+
+           rhos(j)=rhofresh_bs
+           ! blowing snow flux formula used in MAR
+           ustart(j)=ustart0(j)*exp(max(rhoice_bs/rhofresh_bs-rhoice_bs/rhos(j),0.))*exp(max(0.,rhos(j)-rhohard_bs))
+           ! we have multiplied by exp to prevent erosion when rhos>rhohard_bs 
+           ! computation of qbs at the top of the saltation layer
+           ! default formulation from MAR model (Amory et al. 2021, Gallee et al. 2001)
+           esalt=1./(c_esalt_bs*max(1.e-6,ustar(j)))
+           hsalt(j)=0.08436*(max(1.e-6,ustar(j))**1.27)
+           qsalt(j)=(max(ustar(j)**2-ustart(j)**2,0.))/(RG*hsalt(j))*esalt
+           ! calculation of erosion (flux positive towards the surface here) 
+           ! consistent with implicit resolution of turbulent mixing equation
+           ! Nemoto and Nishimura 2004 show that steady-state saltation is achieved within a time tau_eqsalt_bs of about 10s
+           ! we thus prevent snowerosion (snow particle transfer from the saltation layer to the first model level)
+           ! integrated over tau_eqsalt_bs to exceed the total mass of snow particle in the saltation layer
+           ! (rho*qsalt*hsalt)
+           ! during this first step we also lower bound the erosion to the amount of fresh snow accumulated during the time step
+           maxerosion=min(precip_snow(j),hsalt(j)*qsalt(j)*rhod(j)/tau_eqsalt_bs)
+
+           fluxbs(j)=rhod(j)*ws1(j)*cdragm(j)*(AcoefQBS(j)-qsalt(j)) &
+                   / (1.-rhod(j)*ws1(j)*cdragm(j)*BcoefQBS(j)*dtime)
+           !fluxbs(j)= rhod*ws1(j)*cdragm(j)*(qbs1(j)-qsalt(j))
+           fluxbs(j)=max(-maxerosion,fluxbs(j))
+       enddo
+
+
+       ! we now compute the snow age of the overlying layer (snow surface after erosion of the fresh snow accumulated during the time step)
+       ! this is done through the routine albsno
+       CALL albsno(klon,knon,dtime,agesno(:),alb_neig(:), precip_snow(:)+fluxbs(:))
+
+       ! 2nd step:
        ! computation of threshold friction velocity
        ! which depends on surface snow density
@@ -468,76 +498,51 @@
            ! snow density increases with snow age and
            ! increases even faster in case of sedimentation of blowing snow or rain
-           tau_dens=max(tau_densmin, tau_dens0*exp(-abs(precip_bs(j))/pbst_bs &
-                   -abs(precip_rain(j))/prt_bs) *exp(-max(tsurf(j)-272.0,0.)))
-           rhos(j)=rho0+(rhohard-rho0)*(1.-exp(-agesno(j)*86400.0/tau_dens))
+           tau_dens=max(tau_densmin_bs, tau_dens0_bs*exp(-abs(precip_bs(j))/pbst_bs - & 
+                    abs(precip_rain(j))/prt_bs)*exp(-max(tsurf(j)-RTT,0.)))
+           rhos(j)=rhofresh_bs+(rhohard_bs-rhofresh_bs)*(1.-exp(-agesno(j)*86400.0/tau_dens))
            ! blowing snow flux formula used in MAR
-           ws1(j)=(u1(j)**2+v1(j)**2)**0.5
-           ustar(j)=(cdragm(j)*(u1(j)**2+v1(j)**2))**0.5
-           ustart(j)=ustart0*exp(max(rhoice/rho0-rhoice/rhos(j),0.))*exp(max(0.,rhos(j)-rhomax)) 
-           ! we have multiplied by exp to prevent erosion when rhos>rhomax (usefull till
-           ! rhohard<450)
+           ustart(j)=ustart0(j)*exp(max(rhoice_bs/rhofresh_bs-rhoice_bs/rhos(j),0.))*exp(max(0.,rhos(j)-rhohard_bs)) 
+           ! we have multiplied by exp to prevent erosion when rhos>rhohard_bs 
+           ! computation of qbs at the top of the saltation layer
+           ! default formulation from MAR model (Amory et al. 2021, Gallee et al. 2001)
+           esalt=c_esalt_bs*ustar(j)
+           hsalt(j)=0.08436*(max(1.e-6,ustar(j))**1.27)
+           qsalt(j)=(max(ustar(j)**2-ustart(j)**2,0.))/(RG*hsalt(j))*esalt
+           ! calculation of erosion (flux positive towards the surface here) 
+           ! consistent with implicit resolution of turbulent mixing equation
+           ! Nemoto and Nishimura 2004 show that steady-state saltation is achieved within a time tau_eqsalt_bs of about 10s
+           ! we thus prevent snowerosion (snow particle transfer from the saltation layer to the first model level)
+           ! integrated over tau_eqsalt_bs to exceed the total mass of snow particle in the saltation layer
+           ! (rho*qsalt*hsalt)
+           maxerosion=hsalt(j)*qsalt(j)*rhod(j)/tau_eqsalt_bs
+           fluxbs_2=rhod(j)*ws1(j)*cdragm(j)*(AcoefQBS(j)-qsalt(j)) &
+                   / (1.-rhod(j)*ws1(j)*cdragm(j)*BcoefQBS(j)*dtime)
+           !fluxbs_2= rhod*ws1(j)*cdragm(j)*(qbs1(j)-qsalt(j))
+           fluxbs_2=max(-maxerosion,fluxbs_2)
+           fluxbs(j)=fluxbs(j)+fluxbs_2
        enddo
-       
-       ! computation of qbs at the top of the saltation layer
-       ! two formulations possible
-       ! pay attention that qbs is a mixing ratio and has to be converted
-       ! to specific content
-       
-       if (iflag_saltation_bs .eq. 1) then
-       ! expression from CRYOWRF (Sharma et al. 2022)
-          aa=2.6
-          bb=2.5
-          cc=2.0
-          lambdasalt=0.45
-          do j =1, knon
-               rhod=p1lay(j)/RD/temp_air(j)
-               nunu=max(ustar(j)/ustart(j),1.e-3)
-               fluxsalt=rhod/RG*(ustar(j)**3)*(1.-nunu**(-2)) * &
-                        (aa+bb*nunu**(-2)+cc*nunu**(-1))  
-               csalt=fluxsalt/(2.8*ustart(j))
-               hsalt(j)=0.08436*ustar(j)**1.27
-               qsalt(j)=1./rhod*csalt*lambdasalt*RG/(max(ustar(j)**2,1E-6)) &
-                       * exp(-lambdasalt*RG*hsalt(j)/max(ustar(j)**2,1E-6))
-               qsalt(j)=max(qsalt(j),0.)
-          enddo
-
-
-       else
-       ! default formulation from MAR model (Amory et al. 2021, Gallee et al. 2001)        
-          do j=1, knon
-              esalt=1./(3.25*max(ustar(j),0.001))
-              hsalt(j)=0.08436*ustar(j)**1.27
-              qsalt(j)=(max(ustar(j)**2-ustart(j)**2,0.))/(RG*hsalt(j))*esalt
-              !ep=qsalt*cdragm(j)*sqrt(u1(j)**2+v1(j)**2)
-          enddo
-       endif
-
-        ! calculation of erosion (flux positive towards the surface here) 
-        ! consistent with implicit resolution of turbulent mixing equation
-       do j=1, knon
+
+       ! for outputs       
+        do j=1, knon
               i = knindex(j)
-              rhod=p1lay(j)/RD/temp_air(j)
-              ! Nemoto and Nishimura 2004 show that steady-state saltation is achieved within an interval of about 10s
-              ! we thus prevent snowerosion (snow particle transfer from the saltation layer to the first model level)
-              ! to exceed (in intensity) integrated over 10s to exceed the total mass of snow particle in the saltation layer
-              ! (rho*qsalt*hsalt)
-              maxerosion=hsalt(j)*qsalt(j)*rhod/10.
-              fluxbs(j)=rhod*ws1(j)*cdragm(j)*(AcoefQBS(j)-qsalt(j)) &
-                       / (1.-rhod*ws1(j)*cdragm(j)*BcoefQBS(j)*dtime)
-              !fluxbs(j)= rhod*ws1(j)*cdragm(j)*(qbs1(j)-qsalt(j))
-              fluxbs(j)=max(-maxerosion,fluxbs(j))
-
-              ! for outputs
-
               zxustartlic(i) = ustart(j)
               zxrhoslic(i) = rhos(j)
-       enddo
-
-  endif
-
-
-
-!****************************************************************************************
-! Calculate surface snow amount
+              zxqsaltlic(i)=qsalt(j)
+        enddo
+
+
+  else
+  ! those lines are useful to calculate the snow age
+       CALL albsno(klon,knon,dtime,agesno(:),alb_neig(:), precip_snow(:))
+
+  endif ! if ok_bs
+
+
+
+
+
+
+!****************************************************************************************
+! Calculate snow amount
 !    
 !****************************************************************************************
