Index: /LMDZ6/trunk/libf/phylmd/blowing_snow_ini_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/blowing_snow_ini_mod.F90	(revision 4671)
+++ /LMDZ6/trunk/libf/phylmd/blowing_snow_ini_mod.F90	(revision 4672)
@@ -8,5 +8,5 @@
    real coef_eva_bs, expo_eva_bs, fallv_bs, zeta_bs
    real prt_bs, pbst_bs, qbst_bs
-   integer :: niter_bs, iflag_saltation_bs
+   integer :: iflag_saltation_bs
 
    !$OMP THREADPRIVATE(prt_level, lunout)
@@ -14,5 +14,5 @@
    !$OMP THREADPRIVATE(coef_eva_bs, expo_eva_bs, fallv_bs, zeta_bs)
    !$OMP THREADPRIVATE(pbst_bs, prt_bs, qbst_bs)
-   !$OMP THREADPRIVATE(niter_bs, iflag_saltation_bs)
+   !$OMP THREADPRIVATE(iflag_saltation_bs)
 
       contains
@@ -45,5 +45,5 @@
          CALL getin_p('qbst_bs',qbst_bs)
 
-         pbst_bs= 0.01
+         pbst_bs= 0.0003
          CALL getin_p('pbst_bs',pbst_bs)
 
@@ -64,7 +64,4 @@
          CALL getin_p('expo_eva_bs',expo_eva_bs)
 
-         niter_bs = 5
-         CALL getin_p('niter_bs',niter_bs)
-
          iflag_saltation_bs=0
          CALL getin_p('iflag_saltation_bs',iflag_saltation_bs)
Index: /LMDZ6/trunk/libf/phylmd/blowing_snow_sublim_sedim.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/blowing_snow_sublim_sedim.F90	(revision 4671)
+++ /LMDZ6/trunk/libf/phylmd/blowing_snow_sublim_sedim.F90	(revision 4672)
@@ -9,5 +9,5 @@
 
 use blowing_snow_ini_mod, only : coef_eva_bs,RTT,RD,RG,expo_eva_bs, fallv_bs
-use blowing_snow_ini_mod, only : RCPD, RLSTT, RLMLT, RLVTT, RVTMP2, niter_bs
+use blowing_snow_ini_mod, only : RCPD, RLSTT, RLMLT, RLVTT, RVTMP2
 USE lmdz_lscp_tools, only : calc_qsat_ecmwf
 
@@ -125,5 +125,7 @@
               ! vapor, temperature, precip fluxes update
               zq(i) = zq(i) - (sedimn(i)-sedim(i))*(RG/(paprs(i,k)-paprs(i,k+1)))*dtime
-              zqbs(i) = zqbs(i) + (sedimn(i)-sedim(i)) * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
+              zq(i) = max(0., zq(i))
+              !zqbs(i) = zqbs(i) + (sedimn(i)-sedim(i)) * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
+              !zqbs(i) = max(0., zqbs(i))
               ! melting
               zt(i) = zt(i)                             &
@@ -151,7 +153,9 @@
 
 
-              ! vapor, temperature, precip fluxes update
+              ! vapor, temperature, precip fluxes update following sublimation
               zq(i) = zq(i) - (sedimn(i)-sedim(i))*(RG/(paprs(i,k)-paprs(i,k+1)))*dtime                
-              zqbs(i) = zqbs(i) + (sedimn(i)-sedim(i)) * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
+              zq(i) = max(0., zq(i))
+              !zqbs(i) = zqbs(i) + (sedimn(i)-sedim(i)) * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
+              !zqbs(i) = max(0., zqbs(i))
               zt(i) = zt(i)                             &
                 + (sedimn(i)-sedim(i))                      &
@@ -173,20 +177,40 @@
     ENDDO ! loop on ngrid
 
-    ! Now sedimention scheme (alike that in lscp)
-    DO n = 1, niter_bs
-       DO i = 1, ngrid
-         zrho(i,k)  = pplay(i,k) / zt(i) / RD
-         zdz(i)   = (paprs(i,k)-paprs(i,k+1)) / (zrho(i,k)*RG)
-         velo(i,k) = fallv_bs 
-         dqsedim = dtime/REAL(niter_bs)/zdz(i)*zqbs(i)*velo(i,k)   ! dqice/dt=1/rho*d(rho*wice*qice)/dz
-         precbs   = MIN(MAX(dqsedim,0.0),zqbs(i))
-         zqbs(i) = MAX(zqbs(i)-1.*precbs  , 0.0)
-       ENDDO !loop on ngrid
-    ENDDO ! loop on niter_bs
-
-    ! add to non sublimated precip
-    DO i=1,ngrid
-    sedim(i) = sedim(i)+max(0.,zqbsi(i)-zqbs(i))*(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
+    ! Now sedimention scheme 
+
+    ! exact resolution of the conservation equation for qbs with the updated flux from the top (after evap)
+    ! valid only if the fall velocity is constant
+
+    DO i = 1, ngrid
+   
+        zrho(i,k)  = pplay(i,k) / zt(i) / RD
+        zdz(i)   = (paprs(i,k)-paprs(i,k+1)) / (zrho(i,k)*RG)
+        velo(i,k) = fallv_bs
+        zqbs(i) = zqbsi(i)*exp(-velo(i,k)/zdz(i)*dtime)+sedim(i)/zrho(i,k)/velo(i,k)
+        zqbs(i) = max(zqbs(i),0.)
+        ! flux update
+        sedim(i) = sedim(i) + zrho(i,k)*zdz(i)/dtime*(zqbsi(i)-zqbs(i))
+        sedim(i) = max(0.,sedim(i))
+ 
     ENDDO
+
+! old version with bugs
+!    DO n = 1, niter_bs
+!       DO i = 1, ngrid
+!         zrho(i,k)  = pplay(i,k) / zt(i) / RD
+!         zdz(i)   = (paprs(i,k)-paprs(i,k+1)) / (zrho(i,k)*RG)
+!         velo(i,k) = fallv_bs 
+!         dqsedim = dtime/REAL(niter_bs)/zdz(i)*zqbs(i)*velo(i,k)   ! dqice/dt=1/rho*d(rho*wice*qice)/dz
+!         precbs   = MIN(MAX(dqsedim,0.0),zqbs(i))
+!         zqbs(i) = MAX(zqbs(i)-1.*precbs  , 0.0)
+!       ENDDO !loop on ngrid
+!    ENDDO ! loop on niter_bs
+!    ! add to non sublimated precip
+!    DO i=1,ngrid
+!       sedim(i) = sedim(i)+max(0.,zqbsi(i)-zqbs(i))*(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
+!    ENDDO
+!
+
+
 
     ! Outputs:
Index: /LMDZ6/trunk/libf/phylmd/surf_landice_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/surf_landice_mod.F90	(revision 4671)
+++ /LMDZ6/trunk/libf/phylmd/surf_landice_mod.F90	(revision 4672)
@@ -140,8 +140,10 @@
     REAL,DIMENSION(klon) :: precip_totsnow, evap_totsnow
     REAL, DIMENSION (klon,6) :: alb6
-    REAL                   :: rho0, rhoice, ustart0, hsalt, esalt, rhod
+    REAL                   :: rho0, rhoice, ustart0,esalt, rhod
     REAL                   :: lambdasalt,fluxsalt, csalt, nunu, aa, bb, cc
     REAL                   :: tau_dens, tau_dens0, tau_densmin, rhomax, rhohard
-    REAL, DIMENSION(klon)  :: ws1, rhos, ustart, qsalt
+    REAL, DIMENSION(klon)  :: ws1, rhos, ustart, qsalt, hsalt, snowini
+    REAL                   :: maxerosion ! in kg/s   
+
 ! End definition
 !****************************************************************************************
@@ -364,22 +366,24 @@
        ustart0 = 0.211
        rhoice = 920.0
-       rho0 = 200.0
+       rho0 = 300.0
        rhomax=450.0
-       rhohard=400.0
+       rhohard=450.0
        tau_dens0=86400.0*10.  ! 10 days by default, in s
-       tau_densmin=86400.0 ! 1 days according to in situ obs by C. Amory
+       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)
 
        ! computation of threshold friction velocity
        ! which depends on surface snow density
-       do i = 1, knon
+       do j = 1, knon
            ! estimation of snow density
            ! snow density increases with snow age and
-           ! increases even faster in case of sedimentation of blowing snow
-           tau_dens=max(tau_densmin, tau_dens0*exp(-abs(precip_bs(i))/pbst_bs-abs(precip_rain(i))/prt_bs))
-           rhos(i)=rho0+(rhohard-rho0)*(1.-exp(-agesno(i)*86400.0/tau_dens))
+           ! 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))
            ! blowing snow flux formula used in MAR
-           ws1(i)=(u1(i)**2+v1(i)**2)**0.5
-           ustar(i)=(cdragm(i)*(u1(i)**2+v1(i)**2))**0.5
-           ustart(i)=ustart0*exp(max(rhoice/rho0-rhoice/rhos(i),0.))*exp(max(0.,rhos(i)-rhomax)) 
+           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)
@@ -397,14 +401,14 @@
           cc=2.0
           lambdasalt=0.45
-          do i =1, knon
-               rhod=p1lay(i)/RD/temp_air(i)
-               nunu=max(ustar(i)/ustart(i),1.e-3)
-               fluxsalt=rhod/RG*(ustar(i)**3)*(1.-nunu**(-2)) * &
+          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(i))
-               hsalt=0.08436*ustar(i)**1.27
-               qsalt(i)=1./rhod*csalt*lambdasalt*RG/(max(ustar(i)**2,1E-6)) &
-                       * exp(-lambdasalt*RG*hsalt/max(ustar(i)**2,1E-6))
-               qsalt(i)=max(qsalt(i),0.)
+               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
 
@@ -412,28 +416,33 @@
        else
        ! default formulation from MAR model (Amory et al. 2021, Gallee et al. 2001)        
-          do i=1, knon
-              esalt=1./(3.25*max(ustar(i),0.001))
-              hsalt=0.08436*ustar(i)**1.27
-              qsalt(i)=(max(ustar(i)**2-ustart(i)**2,0.))/(RG*hsalt)*esalt
-              !ep=qsalt*cdragm(i)*sqrt(u1(i)**2+v1(i)**2)
+          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 (emission flux towards the first atmospheric level)
+        ! calculation of erosion (flux positive towards the surface here) 
         ! consistent with implicit resolution of turbulent mixing equation
-       do i=1, knon
-              rhod=p1lay(i)/RD/temp_air(i)
-              fluxbs(i)=rhod*ws1(i)*cdragm(i)*zeta_bs*(AcoefQBS(i)-qsalt(i)) &
-                       / (1.-rhod*ws1(i)*zeta_bs*cdragm(i)*BcoefQBS(i)*dtime)
-              !fluxbs(i)= zeta_bs*rhod*ws1(i)*cdragm(i)*(qbs1(i)-qsalt(i))
+       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
 
-       ! for outputs
-       do j = 1, knon
-          i = knindex(j)
-          zxustartlic(i) = ustart(j)
-          zxrhoslic(i) = rhos(j)
-       enddo
-
   endif
 
@@ -441,5 +450,5 @@
 
 !****************************************************************************************
-! Calculate surface snow amount
+! Calculate snow amount
 !    
 !****************************************************************************************
@@ -451,9 +460,11 @@
       evap_totsnow(:)=evap(:)
     ENDIF
-
+    
+ 
     CALL fonte_neige(knon, is_lic, knindex, dtime, &
          tsurf, precip_rain, precip_totsnow,  &
          snow, qsol, tsurf_new, evap_totsnow)
-
+   
+    
     WHERE (snow(1 : knon) .LT. 0.0001) agesno(1 : knon) = 0.                                         
     zfra(1:knon) = MAX(0.0,MIN(1.0,snow(1:knon)/(snow(1:knon)+10.0)))  
Index: /LMDZ6/trunk/libf/phylmdiso/surf_landice_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmdiso/surf_landice_mod.F90	(revision 4671)
+++ /LMDZ6/trunk/libf/phylmdiso/surf_landice_mod.F90	(revision 4672)
@@ -185,8 +185,10 @@
     REAL,DIMENSION(klon) :: precip_totsnow, evap_totsnow
     REAL, DIMENSION (klon,6) :: alb6
-    REAL                   :: rho0, rhoice, ustart0, hsalt, esalt, rhod
+    REAL                   :: rho0, rhoice, ustart0,esalt, rhod
     REAL                   :: lambdasalt,fluxsalt, csalt, nunu, aa, bb, cc
     REAL                   :: tau_dens, tau_dens0, tau_densmin, rhomax, rhohard
-    REAL, DIMENSION(klon)  :: ws1, rhos, ustart, qsalt
+    REAL, DIMENSION(klon)  :: ws1, rhos, ustart, qsalt, hsalt, snowini
+    REAL                   :: maxerosion ! in kg/s
+
 ! End definition
 !****************************************************************************************
@@ -453,22 +455,24 @@
        ustart0 = 0.211
        rhoice = 920.0
-       rho0 = 200.0
+       rho0 = 300.0
        rhomax=450.0
-       rhohard=400.0
+       rhohard=450.0
        tau_dens0=86400.0*10.  ! 10 days by default, in s
-       tau_densmin=86400.0 ! 1 days according to in situ obs by C. Amory
+       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)
 
        ! computation of threshold friction velocity
        ! which depends on surface snow density
-       do i = 1, knon
+       do j = 1, knon
            ! estimation of snow density
            ! snow density increases with snow age and
-           ! increases even faster in case of sedimentation of blowing snow
-           tau_dens=max(tau_densmin, tau_dens0*exp(-abs(precip_bs(i))/pbst_bs-abs(precip_rain(i))/prt_bs))
-           rhos(i)=rho0+(rhohard-rho0)*(1.-exp(-agesno(i)*86400.0/tau_dens))
+           ! 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))
            ! blowing snow flux formula used in MAR
-           ws1(i)=(u1(i)**2+v1(i)**2)**0.5
-           ustar(i)=(cdragm(i)*(u1(i)**2+v1(i)**2))**0.5
-           ustart(i)=ustart0*exp(max(rhoice/rho0-rhoice/rhos(i),0.))*exp(max(0.,rhos(i)-rhomax)) 
+           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)
@@ -486,14 +490,14 @@
           cc=2.0
           lambdasalt=0.45
-          do i =1, knon
-               rhod=p1lay(i)/RD/temp_air(i)
-               nunu=max(ustar(i)/ustart(i),1.e-3)
-               fluxsalt=rhod/RG*(ustar(i)**3)*(1.-nunu**(-2)) * &
+          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(i))
-               hsalt=0.08436*ustar(i)**1.27
-               qsalt(i)=1./rhod*csalt*lambdasalt*RG/(max(ustar(i)**2,1E-6)) &
-                       * exp(-lambdasalt*RG*hsalt/max(ustar(i)**2,1E-6))
-               qsalt(i)=max(qsalt(i),0.)
+               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
 
@@ -501,26 +505,31 @@
        else
        ! default formulation from MAR model (Amory et al. 2021, Gallee et al. 2001)        
-          do i=1, knon
-              esalt=1./(3.25*max(ustar(i),0.001))
-              hsalt=0.08436*ustar(i)**1.27
-              qsalt(i)=(max(ustar(i)**2-ustart(i)**2,0.))/(RG*hsalt)*esalt
-              !ep=qsalt*cdragm(i)*sqrt(u1(i)**2+v1(i)**2)
+          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 (emission flux towards the first atmospheric level)
+        ! calculation of erosion (flux positive towards the surface here) 
         ! consistent with implicit resolution of turbulent mixing equation
-       do i=1, knon
-              rhod=p1lay(i)/RD/temp_air(i)
-              fluxbs(i)=rhod*ws1(i)*cdragm(i)*zeta_bs*(AcoefQBS(i)-qsalt(i)) &
-                       / (1.-rhod*ws1(i)*zeta_bs*cdragm(i)*BcoefQBS(i)*dtime)
-              !fluxbs(i)= zeta_bs*rhod*ws1(i)*cdragm(i)*(qbs1(i)-qsalt(i))
-       enddo
-
-       ! for outputs
-       do j = 1, knon
-          i = knindex(j)
-          zxustartlic(i) = ustart(j)
-          zxrhoslic(i) = rhos(j)
+       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
 
