Index: /LMDZ6/trunk/libf/phylmd/lmdz_blowing_snow_sublim_sedim.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/lmdz_blowing_snow_sublim_sedim.F90	(revision 4827)
+++ /LMDZ6/trunk/libf/phylmd/lmdz_blowing_snow_sublim_sedim.F90	(revision 4828)
@@ -110,4 +110,10 @@
     DO i = 1, ngrid
 
+          rhoair  = pplay(i,k) / zt(i) / RD
+          dz      = (paprs(i,k)-paprs(i,k+1)) / (rhoair*RG)
+          ! BS fall velocity assumed to be constant for now
+          zvelo(i) = fallv_bs
+
+
           IF (zt(i) .GT. RTT) THEN
              
@@ -117,7 +123,7 @@
            
              taumeltbs=taumeltbs0*exp(-max(0.,(zt(i)-RTT)/(tbsmelt-RTT)))
-             dqmelt=min(zqbs(i),zqbs(i)/taumeltbs*dtime)
+             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(max(dqmelt,0.),maxdqmelt)
+             dqmelt=min(dqmelt,maxdqmelt)
              ! q update, melting + evaporation
              zq(i) = zq(i) + dqmelt
@@ -130,5 +136,5 @@
               ! 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
@@ -148,9 +154,5 @@
 
         ! Sedimentation scheme
-           
-        rhoair  = pplay(i,k) / zt(i) / RD
-        dz      = (paprs(i,k)-paprs(i,k+1)) / (rhoair*RG)
-        ! BS fall velocity assumed to be constant for now
-        zvelo(i) = fallv_bs
+        !----------------------
 
         sedimn(i) = rhoair*zqbs(i)*zvelo(i)
@@ -173,5 +175,5 @@
                  zt(i) = zt(i) - zqbs(i) * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zqbs(i)))
               ELSE
-               zt(i) = zt(i) - zqbs(i) * (RLVTT+RLMLT)/RCPD/(1.0+RVTMP2*(zq(i)+zqbs(i)))
+                 zt(i) = zt(i) - zqbs(i) * (RLVTT+RLMLT)/RCPD/(1.0+RVTMP2*(zq(i)+zqbs(i)))
               ENDIF
               zqbs(i)=0.
Index: /LMDZ6/trunk/libf/phylmd/surf_landice_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/surf_landice_mod.F90	(revision 4827)
+++ /LMDZ6/trunk/libf/phylmd/surf_landice_mod.F90	(revision 4828)
@@ -144,5 +144,5 @@
     REAL                   :: tau_dens, tau_dens0, tau_densmin, rhomax, rhohard
     REAL, DIMENSION(klon)  :: ws1, rhos, ustart, qsalt, hsalt, snowini
-    REAL                   :: maxerosion ! in kg/s   
+    REAL                   :: ustarmin, maxerosion, tau_eq_salt 
 
 ! End definition
@@ -364,4 +364,5 @@
   ! Simple blowing snow param
   if (ok_bs) then
+       ustarmin=1.e-3
        ustart0 = 0.211
        rhoice = 920.0
@@ -372,4 +373,5 @@
        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.
 
        ! computation of threshold friction velocity
@@ -408,6 +410,6 @@
                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)=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
@@ -417,6 +419,6 @@
        ! 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
+              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)
@@ -429,14 +431,13 @@
               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
+              ! 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)
-              ! to exceed (in intensity) integrated over 10s to exceed the total mass of snow particle in the saltation layer
+              ! 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/10.
+              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)=max(-maxerosion,fluxbs(j))
-
+              fluxbs(j)=min(0.,max(-maxerosion,fluxbs(j)))
               ! for outputs
 
