Index: LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmdiso/add_phys_tend_mod.F90
===================================================================
--- LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmdiso/add_phys_tend_mod.F90	(revision 4660)
+++ LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmdiso/add_phys_tend_mod.F90	(revision 4669)
@@ -16,5 +16,5 @@
 CONTAINS
 
-SUBROUTINE add_pbl_tend(zdu, zdv, zdt, zdq, zdql, zdqi, paprs, text,abortphy,flag_inhib_tend, itap &
+SUBROUTINE add_pbl_tend(zdu, zdv, zdt, zdq, zdql, zdqi, zdqbs, paprs, text,abortphy,flag_inhib_tend, itap &
 #ifdef ISO
         ,zdxt,zdxtl,zdxti &
@@ -62,5 +62,5 @@
   ! ------------
   REAL zdu(klon, klev), zdv(klon, klev)
-  REAL zdt(klon, klev), zdq(klon, klev), zdql(klon, klev), zdqi(klon, klev)
+  REAL zdt(klon, klev), zdq(klon, klev), zdql(klon, klev), zdqi(klon, klev), zdqbs(klon,klev)
   CHARACTER *(*) text
   REAL paprs(klon,klev+1)
@@ -104,5 +104,5 @@
     PRINT *, ' add_pbl_tend, zzdt ', zzdt
     PRINT *, ' add_pbl_tend, zzdq ', zzdq
-    CALL add_phys_tend(zdu, zdv, zzdt, zzdq, zdql, zdqi, paprs, text,abortphy,flag_inhib_tend, itap, 0 &
+    CALL add_phys_tend(zdu, zdv, zzdt, zzdq, zdql, zdqi, zdqbs, paprs, text,abortphy,flag_inhib_tend, itap, 0 &
 #ifdef ISO
         ,zzdxt,zdxtl,zdxti &
@@ -110,5 +110,5 @@
         )
   ELSE
-    CALL add_phys_tend(zdu, zdv, zdt, zdq, zdql, zdqi, paprs, text,abortphy,flag_inhib_tend, itap, 0 &
+    CALL add_phys_tend(zdu, zdv, zdt, zdq, zdql, zdqi, zdqbs, paprs, text,abortphy,flag_inhib_tend, itap, 0 &
 #ifdef ISO
         ,zdxt,zdxtl,zdxti &
@@ -123,5 +123,5 @@
 ! $Id$
 !
-SUBROUTINE add_phys_tend (zdu,zdv,zdt,zdq,zdql,zdqi,paprs,text, &
+SUBROUTINE add_phys_tend (zdu,zdv,zdt,zdq,zdql,zdqi,zdqbs,paprs,text, &
                           abortphy,flag_inhib_tend, itap, diag_mode &
 #ifdef ISO
@@ -142,5 +142,5 @@
 USE dimphy, ONLY: klon, klev
 USE phys_state_var_mod, ONLY : phys_tstep
-USE phys_local_var_mod, ONLY: u_seri, v_seri, ql_seri, qs_seri, q_seri, t_seri
+USE phys_local_var_mod, ONLY: u_seri, v_seri, ql_seri, qs_seri, qbs_seri, q_seri, t_seri
 #ifdef ISO
 USE phys_local_var_mod, ONLY: xtl_seri, xts_seri, xt_seri
@@ -150,6 +150,6 @@
 USE print_control_mod, ONLY: prt_level
 USE cmp_seri_mod
-USE phys_output_var_mod, ONLY : d_qw_col, d_ql_col, d_qs_col, d_qt_col, d_ek_col, d_h_dair_col &
-  &           , d_h_qw_col, d_h_ql_col, d_h_qs_col, d_h_col
+USE phys_output_var_mod, ONLY : d_qw_col, d_ql_col, d_qs_col, d_qbs_col, d_qt_col, d_ek_col, d_h_dair_col &
+  &           , d_h_qw_col, d_h_ql_col, d_h_qs_col, d_h_qbs_col, d_h_col
 
 #ifdef ISO
@@ -167,5 +167,5 @@
 !------------
 REAL, DIMENSION(klon,klev),     INTENT(IN)    :: zdu, zdv
-REAL, DIMENSION(klon,klev),     INTENT(IN)    :: zdt, zdql, zdqi
+REAL, DIMENSION(klon,klev),     INTENT(IN)    :: zdt, zdql, zdqi, zdqbs
 REAL, DIMENSION(klon,klev+1),   INTENT(IN)    :: paprs
 CHARACTER*(*),                  INTENT(IN)    :: text
@@ -197,5 +197,5 @@
 ! Save variables, used in diagnostic mode (diag_mode=1).
 REAL, DIMENSION(klon,klev)   :: sav_u_seri, sav_v_seri
-REAL, DIMENSION(klon,klev)   :: sav_ql_seri, sav_qs_seri, sav_q_seri
+REAL, DIMENSION(klon,klev)   :: sav_ql_seri, sav_qs_seri, sav_qbs_seri, sav_q_seri
 REAL, DIMENSION(klon,klev)   :: sav_t_seri
 REAL, DIMENSION(klon,klev)   :: sav_zdq
@@ -228,7 +228,9 @@
 ! zh_ql_col----  total enthalpy of liquid watter (J/m2)
 ! zh_qs_col----  total enthalpy of solid watter  (J/m2)
+! zh_qbs_col----  total enthalpy of blowing snow  (J/m2)
 ! zqw_col------  total mass of watter vapour (kg/m2)
 ! zql_col------  total mass of liquid watter (kg/m2)
-! zqs_col------  total mass of solid watter (kg/m2)
+! zqs_col------  total mass of cloud ice (kg/m2)
+! zqbs_col------  total mass of blowing snow (kg/m2)
 ! zek_col------  total kinetic energy (kg/m2)
 !
@@ -237,7 +239,8 @@
 REAL zql_col(klon,2)
 REAL zqs_col(klon,2)
+REAL zqbs_col(klon,2)
 REAL zek_col(klon,2)
 REAL zh_dair_col(klon,2)
-REAL zh_qw_col(klon,2), zh_ql_col(klon,2), zh_qs_col(klon,2)
+REAL zh_qw_col(klon,2), zh_ql_col(klon,2), zh_qs_col(klon,2), zh_qbs_col(klon,2)
 REAL zh_col(klon,2)
 
@@ -278,4 +281,5 @@
     sav_ql_seri(:,:) = ql_seri(:,:)
     sav_qs_seri(:,:) = qs_seri(:,:)
+    sav_qbs_seri(:,:) = qbs_seri(:,:)
     sav_q_seri(:,:)  = q_seri(:,:)
     sav_t_seri(:,:)  = t_seri(:,:)
@@ -304,7 +308,7 @@
 
     CALL integr_v(klon, klev, zcpvap, &
-                  t_seri, q_seri, ql_seri, qs_seri, u_seri, v_seri, zairm, &
-                    zqw_col(:,n), zql_col(:,n), zqs_col(:,n), zek_col(:,n), zh_dair_col(:,n), &
-                    zh_qw_col(:,n), zh_ql_col(:,n), zh_qs_col(:,n), zh_col(:,n))
+                  t_seri, q_seri, ql_seri, qs_seri, qbs_seri, u_seri, v_seri, zairm, &
+                    zqw_col(:,n), zql_col(:,n), zqs_col(:,n), zqbs_col(:,n), zek_col(:,n), zh_dair_col(:,n), &
+                    zh_qw_col(:,n), zh_ql_col(:,n), zh_qs_col(:,n), zh_qbs_col(:,n), zh_col(:,n))
 
   end if ! end if (fl_ebil .GT. 0)
@@ -337,4 +341,5 @@
      ql_seri(:,:)=ql_seri(:,:)+zdql(:,:)
      qs_seri(:,:)=qs_seri(:,:)+zdqi(:,:)
+     qbs_seri(:,:)=qbs_seri(:,:)+zdqbs(:,:)
 #ifdef ISO 
      xtl_seri(:,:,:)=xtl_seri(:,:,:)+zdxtl(:,:,:) 
@@ -601,4 +606,5 @@
         ql_seri(:,:)=ql_seri(:,:)-zdql(:,:)
         qs_seri(:,:)=qs_seri(:,:)-zdqi(:,:)
+        qbs_seri(:,:)=qbs_seri(:,:)-zdqbs(:,:)
       ENDIF
 
@@ -615,7 +621,7 @@
 
     CALL integr_v(klon, klev, zcpvap, &
-                  t_seri, q_seri, ql_seri, qs_seri, u_seri, v_seri, zairm, &
-                    zqw_col(:,n), zql_col(:,n), zqs_col(:,n), zek_col(:,n), zh_dair_col(:,n), &
-                    zh_qw_col(:,n), zh_ql_col(:,n), zh_qs_col(:,n), zh_col(:,n))
+                  t_seri, q_seri, ql_seri, qs_seri, qbs_seri, u_seri, v_seri, zairm, &
+                    zqw_col(:,n), zql_col(:,n), zqs_col(:,n), zqbs_col(:,n), zek_col(:,n), zh_dair_col(:,n), &
+                    zh_qw_col(:,n), zh_ql_col(:,n), zh_qs_col(:,n), zh_qbs_col(:,n), zh_col(:,n))
 
     ! ------------------------------------------------
@@ -626,5 +632,6 @@
     d_ql_col(:) = (zql_col(:,2)-zql_col(:,1))/phys_tstep
     d_qs_col(:) = (zqs_col(:,2)-zqs_col(:,1))/phys_tstep
-    d_qt_col(:) = d_qw_col(:) + d_ql_col(:) + d_qs_col(:)
+    d_qbs_col(:) = (zqbs_col(:,2)-zqbs_col(:,1))/phys_tstep
+    d_qt_col(:) = d_qw_col(:) + d_ql_col(:) + d_qs_col(:) + d_qbs_col(:)
 
     d_ek_col(:) = (zek_col(:,2)-zek_col(:,1))/phys_tstep
@@ -634,4 +641,5 @@
     d_h_ql_col(:) = (zh_ql_col(:,2)-zh_ql_col(:,1))/phys_tstep
     d_h_qs_col(:) = (zh_qs_col(:,2)-zh_qs_col(:,1))/phys_tstep
+    d_h_qbs_col(:) = (zh_qbs_col(:,2)-zh_qbs_col(:,1))/phys_tstep
 
     d_h_col = (zh_col(:,2)-zh_col(:,1))/phys_tstep 
@@ -645,4 +653,5 @@
     ql_seri(:,:) = sav_ql_seri(:,:)
     qs_seri(:,:) = sav_qs_seri(:,:)
+    qbs_seri(:,:) = sav_qbs_seri(:,:)
     q_seri(:,:)  = sav_q_seri(:,:)
     t_seri(:,:)  = sav_t_seri(:,:)
@@ -659,6 +668,6 @@
 END SUBROUTINE add_phys_tend
 
-SUBROUTINE diag_phys_tend (nlon, nlev, uu, vv, temp, qv, ql, qs, &
-                          zdu,zdv,zdt,zdq,zdql,zdqs,paprs,text)
+SUBROUTINE diag_phys_tend (nlon, nlev, uu, vv, temp, qv, ql, qs, qbs, &
+                          zdu,zdv,zdt,zdq,zdql,zdqs,zdqbs,paprs,text)
 !======================================================================
 ! Ajoute les tendances des variables physiques aux variables 
@@ -676,6 +685,6 @@
 USE print_control_mod, ONLY: prt_level
 USE cmp_seri_mod
-USE phys_output_var_mod, ONLY : d_qw_col, d_ql_col, d_qs_col, d_qt_col, d_ek_col, d_h_dair_col &
-  &           , d_h_qw_col, d_h_ql_col, d_h_qs_col, d_h_col
+USE phys_output_var_mod, ONLY : d_qw_col, d_ql_col, d_qs_col, d_qbs_col, d_qt_col, d_ek_col, d_h_dair_col &
+  &           , d_h_qw_col, d_h_ql_col, d_h_qs_col, d_h_qbs_col, d_h_col
 IMPLICIT none
   include "YOMCST.h"
@@ -686,7 +695,7 @@
 INTEGER, INTENT(IN)                           :: nlon, nlev
 REAL, DIMENSION(nlon,nlev),     INTENT(IN)    :: uu, vv
-REAL, DIMENSION(nlon,nlev),     INTENT(IN)    :: temp, qv, ql, qs
+REAL, DIMENSION(nlon,nlev),     INTENT(IN)    :: temp, qv, ql, qs, qbs
 REAL, DIMENSION(nlon,nlev),     INTENT(IN)    :: zdu, zdv
-REAL, DIMENSION(nlon,nlev),     INTENT(IN)    :: zdt, zdq, zdql, zdqs
+REAL, DIMENSION(nlon,nlev),     INTENT(IN)    :: zdt, zdq, zdql, zdqs, zdqbs
 REAL, DIMENSION(nlon,nlev+1),   INTENT(IN)    :: paprs
 CHARACTER*(*),                  INTENT(IN)    :: text
@@ -695,5 +704,5 @@
 !--------
 REAL, DIMENSION(nlon,nlev)      :: uu_n, vv_n
-REAL, DIMENSION(nlon,nlev)      :: temp_n, qv_n, ql_n, qs_n
+REAL, DIMENSION(nlon,nlev)      :: temp_n, qv_n, ql_n, qs_n, qbs_n
 
 
@@ -718,5 +727,6 @@
 ! zqw_col------  total mass of watter vapour (kg/m2)
 ! zql_col------  total mass of liquid watter (kg/m2)
-! zqs_col------  total mass of solid watter (kg/m2)
+! zqs_col------  total mass of cloud ice (kg/m2)
+! zqbs_col------  total mass of blowing snow (kg/m2)
 ! zek_col------  total kinetic energy (kg/m2)
 !
@@ -725,7 +735,8 @@
 REAL zql_col(nlon,2)
 REAL zqs_col(nlon,2)
+REAL zqbs_col(nlon,2)
 REAL zek_col(nlon,2)
 REAL zh_dair_col(nlon,2)
-REAL zh_qw_col(nlon,2), zh_ql_col(nlon,2), zh_qs_col(nlon,2)
+REAL zh_qw_col(nlon,2), zh_ql_col(nlon,2), zh_qs_col(nlon,2), zh_qbs_col(nlon,2)
 REAL zh_col(nlon,2)
 
@@ -760,7 +771,7 @@
 
     CALL integr_v(nlon, nlev, rcpv, &
-                  temp, qv, ql, qs, uu, vv, zairm, &
-                    zqw_col(:,n), zql_col(:,n), zqs_col(:,n), zek_col(:,n), zh_dair_col(:,n), &
-                    zh_qw_col(:,n), zh_ql_col(:,n), zh_qs_col(:,n), zh_col(:,n))
+                  temp, qv, ql, qs, qbs, uu, vv, zairm, &
+                    zqw_col(:,n), zql_col(:,n), zqs_col(:,n), zqbs_col(:,n), zek_col(:,n), zh_dair_col(:,n), &
+                    zh_qw_col(:,n), zh_ql_col(:,n), zh_qs_col(:,n), zh_qbs_col(:,n), zh_col(:,n))
 
   end if ! end if (fl_ebil .GT. 0)
@@ -775,4 +786,5 @@
      ql_n(:,:)=ql(:,:)+zdql(:,:)
      qs_n(:,:)=qs(:,:)+zdqs(:,:)
+     qbs_n(:,:)=qbs(:,:)+zdqbs(:,:)
      temp_n(:,:)=temp(:,:)+zdt(:,:)
 
@@ -791,7 +803,7 @@
 
     CALL integr_v(nlon, nlev, rcpv, &
-                  temp_n, qv_n, ql_n, qs_n, uu_n, vv_n, zairm, &
-                    zqw_col(:,n), zql_col(:,n), zqs_col(:,n), zek_col(:,n), zh_dair_col(:,n), &
-                    zh_qw_col(:,n), zh_ql_col(:,n), zh_qs_col(:,n), zh_col(:,n))
+                  temp_n, qv_n, ql_n, qs_n, qbs_n, uu_n, vv_n, zairm, &
+                    zqw_col(:,n), zql_col(:,n), zqs_col(:,n), zqbs_col(:,n), zek_col(:,n), zh_dair_col(:,n), &
+                    zh_qw_col(:,n), zh_ql_col(:,n), zh_qs_col(:,n), zh_qbs_col(:,n), zh_col(:,n))
 
     ! ------------------------------------------------
@@ -802,5 +814,6 @@
     d_ql_col(:) = (zql_col(:,2)-zql_col(:,1))/phys_tstep
     d_qs_col(:) = (zqs_col(:,2)-zqs_col(:,1))/phys_tstep
-    d_qt_col(:) = d_qw_col(:) + d_ql_col(:) + d_qs_col(:)
+    d_qbs_col(:) = (zqbs_col(:,2)-zqbs_col(:,1))/phys_tstep
+    d_qt_col(:) = d_qw_col(:) + d_ql_col(:) + d_qs_col(:) + d_qbs_col(:)
 
     d_ek_col(:) = (zek_col(:,2)-zek_col(:,1))/phys_tstep
@@ -814,4 +827,5 @@
     d_h_ql_col(:) = (zh_ql_col(:,2)-zh_ql_col(:,1))/phys_tstep
     d_h_qs_col(:) = (zh_qs_col(:,2)-zh_qs_col(:,1))/phys_tstep
+    d_h_qbs_col(:) = (zh_qbs_col(:,2)-zh_qbs_col(:,1))/phys_tstep
 
     d_h_col = (zh_col(:,2)-zh_col(:,1))/phys_tstep 
@@ -824,7 +838,7 @@
 
 SUBROUTINE integr_v(nlon, nlev, zcpvap, &
-                    temp, qv, ql, qs, uu, vv, zairm,  &
-                    zqw_col, zql_col, zqs_col, zek_col, zh_dair_col, &
-                    zh_qw_col, zh_ql_col, zh_qs_col, zh_col)
+                    temp, qv, ql, qs, qbs, uu, vv, zairm,  &
+                    zqw_col, zql_col, zqs_col, zqbs_col, zek_col, zh_dair_col, &
+                    zh_qw_col, zh_ql_col, zh_qs_col, zh_qbs_col, zh_col)
 
 IMPLICIT none
@@ -833,14 +847,14 @@
 INTEGER,                    INTENT(IN)    :: nlon,nlev
 REAL,                       INTENT(IN)    :: zcpvap
-REAL, DIMENSION(nlon,nlev), INTENT(IN)    :: temp, qv, ql, qs, uu, vv
+REAL, DIMENSION(nlon,nlev), INTENT(IN)    :: temp, qv, ql, qs, qbs, uu, vv
 REAL, DIMENSION(nlon,nlev), INTENT(IN)    :: zairm
 REAL, DIMENSION(nlon),      INTENT(OUT)   :: zqw_col
 REAL, DIMENSION(nlon),      INTENT(OUT)   :: zql_col
-REAL, DIMENSION(nlon),      INTENT(OUT)   :: zqs_col
+REAL, DIMENSION(nlon),      INTENT(OUT)   :: zqs_col, zqbs_col
 REAL, DIMENSION(nlon),      INTENT(OUT)   :: zek_col
 REAL, DIMENSION(nlon),      INTENT(OUT)   :: zh_dair_col
 REAL, DIMENSION(nlon),      INTENT(OUT)   :: zh_qw_col
 REAL, DIMENSION(nlon),      INTENT(OUT)   :: zh_ql_col
-REAL, DIMENSION(nlon),      INTENT(OUT)   :: zh_qs_col
+REAL, DIMENSION(nlon),      INTENT(OUT)   :: zh_qs_col, zh_qbs_col
 REAL, DIMENSION(nlon),      INTENT(OUT)   :: zh_col
 
@@ -852,4 +866,5 @@
     zql_col(:) = 0.
     zqs_col(:) = 0.
+    zqbs_col(:) = 0.
     zek_col(:) = 0.
     zh_dair_col(:) = 0.
@@ -857,4 +872,5 @@
     zh_ql_col(:) = 0.
     zh_qs_col(:) = 0.
+    zh_qbs_col(:) = 0.
 
 !JLD    write (*,*) "rcpd, zcpvap, zcwat, zcice ",rcpd, zcpvap, zcwat, zcice
@@ -869,4 +885,5 @@
         zql_col(i) = zql_col(i) + ql(i, k)*zairm(i, k)
         zqs_col(i) = zqs_col(i) + qs(i, k)*zairm(i, k)
+        zqbs_col(i)= zqbs_col(i) + qbs(i,k)*zairm(i,k)
         ! Kinetic Energy
         zek_col(i) = zek_col(i) + 0.5*(uu(i,k)**2+vv(i,k)**2)*zairm(i, k)
@@ -877,8 +894,9 @@
         zh_ql_col(i) = zh_ql_col(i) + (zcpvap*temp(i, k) - rlvtt)*ql(i, k)*zairm(i, k)   !jyg
         zh_qs_col(i) = zh_qs_col(i) + (zcpvap*temp(i, k) - rlstt)*qs(i, k)*zairm(i, k)   !jyg
+        zh_qbs_col(i) = zh_qbs_col(i) + (zcpvap*temp(i, k) - rlstt)*qbs(i, k)*zairm(i, k)   !jyg
       END DO
     END DO
     ! compute total air enthalpy
-    zh_col(:) = zh_dair_col(:) + zh_qw_col(:) + zh_ql_col(:) + zh_qs_col(:)
+    zh_col(:) = zh_dair_col(:) + zh_qw_col(:) + zh_ql_col(:) + zh_qs_col(:) + zh_qbs_col(:)
 
 END SUBROUTINE integr_v
@@ -895,12 +913,12 @@
 USE dimphy, ONLY: klon, klev
 USE phys_state_var_mod, ONLY : phys_tstep
-USE phys_state_var_mod, ONLY : topsw, toplw, solsw, sollw, rain_con, snow_con
+USE phys_state_var_mod, ONLY : topsw, toplw, solsw, sollw, rain_con, snow_con, bs_fall
 USE geometry_mod, ONLY: longitude_deg, latitude_deg
 USE print_control_mod, ONLY: prt_level
 USE cmp_seri_mod
-USE phys_output_var_mod, ONLY : d_qw_col, d_ql_col, d_qs_col, d_qt_col, d_ek_col, d_h_dair_col &
-  &           , d_h_qw_col, d_h_ql_col, d_h_qs_col, d_h_col
+USE phys_output_var_mod, ONLY : d_qw_col, d_ql_col, d_qs_col, d_qbs_col, d_qt_col, d_ek_col, d_h_dair_col &
+  &           , d_h_qw_col, d_h_ql_col, d_h_qs_col, d_h_qbs_col, d_h_col
 USE phys_local_var_mod, ONLY: evap, sens
-USE phys_local_var_mod, ONLY: u_seri, v_seri, ql_seri, qs_seri, q_seri, t_seri &
+USE phys_local_var_mod, ONLY: u_seri, v_seri, ql_seri, qs_seri, qbs_seri, q_seri, t_seri &
    &  , rain_lsc, snow_lsc
 USE climb_hq_mod, ONLY : d_h_col_vdf, f_h_bnd
@@ -939,4 +957,7 @@
       bilh_bnd = (-(rcw-rcpd)*t_seri(1,1) + rlvtt) * rain_lsc(1) &
     &         + (-(rcs-rcpd)*t_seri(1,1) + rlstt) * snow_lsc(1)
+  CASE("bs") param
+      bilq_bnd = - bs_fall(1)
+      bilh_bnd = (-(rcs-rcpd)*t_seri(1,1) + rlstt) * bs_fall(1)
   CASE("convection") param
       bilq_bnd = - rain_con(1) - snow_con(1)
@@ -975,6 +996,6 @@
   if ( prt_level .GE. 5) then
     write(*,9000) text,"enerbil at boundaries: Q, H",bilq_bnd, bilh_bnd
-    write(*,9000) text,"enerbil: water budget",d_qt_col(1),d_qw_col(1),d_ql_col(1),d_qs_col(1)
-    write(*,9000) text,"enerbil: enthalpy budget",d_h_col(1),d_h_dair_col(1),d_h_qw_col(1),d_h_ql_col(1),d_h_qs_col(1)
+    write(*,9000) text,"enerbil: water budget",d_qt_col(1),d_qw_col(1),d_ql_col(1),d_qs_col(1), d_qbs_col(1)
+    write(*,9000) text,"enerbil: enthalpy budget",d_h_col(1),d_h_dair_col(1),d_h_qw_col(1),d_h_ql_col(1),d_h_qs_col(1),d_h_qbs_col(1)
   end if
 
Index: LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmdiso/blowing_snow_ini_mod.F90
===================================================================
--- LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmdiso/blowing_snow_ini_mod.F90	(revision 4669)
+++ LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmdiso/blowing_snow_ini_mod.F90	(revision 4669)
@@ -0,0 +1,1 @@
+link ../phylmd/blowing_snow_ini_mod.F90
Index: LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmdiso/blowing_snow_sublim_sedim.F90
===================================================================
--- LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmdiso/blowing_snow_sublim_sedim.F90	(revision 4669)
+++ LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmdiso/blowing_snow_sublim_sedim.F90	(revision 4669)
@@ -0,0 +1,1 @@
+link ../phylmd/blowing_snow_sublim_sedim.F90
Index: LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmdiso/call_atke_mod.F90
===================================================================
--- LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmdiso/call_atke_mod.F90	(revision 4669)
+++ LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmdiso/call_atke_mod.F90	(revision 4669)
@@ -0,0 +1,1 @@
+link ../phylmd/call_atke_mod.F90
Index: LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmdiso/call_blowing_snow_sublim_sedim.F90
===================================================================
--- LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmdiso/call_blowing_snow_sublim_sedim.F90	(revision 4669)
+++ LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmdiso/call_blowing_snow_sublim_sedim.F90	(revision 4669)
@@ -0,0 +1,1 @@
+link ../phylmd/call_blowing_snow_sublim_sedim.F90
Index: LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmdiso/climb_qbs_mod.F90
===================================================================
--- LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmdiso/climb_qbs_mod.F90	(revision 4669)
+++ LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmdiso/climb_qbs_mod.F90	(revision 4669)
@@ -0,0 +1,1 @@
+link ../phylmd/climb_qbs_mod.F90
Index: LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmdiso/fisrtilp.F90
===================================================================
--- LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmdiso/fisrtilp.F90	(revision 4660)
+++ LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmdiso/fisrtilp.F90	(revision 4669)
@@ -26,4 +26,12 @@
   ! flag to include modifications to ensure energy conservation (if flag >0)
   USE add_phys_tend_mod, only : fl_cor_ebil 
+  USE lscp_ini_mod, ONLY: iflag_t_glace,t_glace_min, t_glace_max, exposant_glace
+  USE lscp_ini_mod, ONLY: iflag_cloudth_vert, iflag_rain_incloud_vol 
+  USE lscp_ini_mod, ONLY: coef_eva, coef_eva_i, ffallv_lsc, ffallv_con
+  USE lscp_ini_mod, ONLY: cld_tau_lsc, cld_tau_con, cld_lc_lsc, cld_lc_con
+  USE lscp_ini_mod, ONLY: reevap_ice, iflag_bergeron, iflag_fisrtilp_qsat, iflag_pdf 
+
+
+
 #ifdef ISO
   USE infotrac_phy, ONLY: ntraciso=>ntiso,niso,itZonIso
@@ -82,7 +90,4 @@
   !======================================================================
   include "YOMCST.h"
-  include "fisrtilp.h"
-  include "nuage.h" ! JBM (3/14)
-
   !
   ! Principaux inputs:
Index: LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmdiso/fisrtilp.h
===================================================================
--- LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmdiso/fisrtilp.h	(revision 4660)
+++ 	(revision )
@@ -1,1 +1,0 @@
-link ../phylmd/fisrtilp.h
Index: LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmdiso/pbl_surface_mod.F90
===================================================================
--- LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmdiso/pbl_surface_mod.F90	(revision 4660)
+++ LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmdiso/pbl_surface_mod.F90	(revision 4669)
@@ -21,4 +21,5 @@
   USE cpl_mod,             ONLY : gath2cpl
   USE climb_hq_mod,        ONLY : climb_hq_down, climb_hq_up
+  USE climb_qbs_mod,       ONLY : climb_qbs_down, climb_qbs_up
   USE climb_wind_mod,      ONLY : climb_wind_down, climb_wind_up
   USE coef_diff_turb_mod,  ONLY : coef_diff_turb
@@ -261,7 +262,7 @@
        rlon,      rlat,      rugoro,   rmu0,          &
        lwdown_m,  cldt,          &
-       rain_f,    snow_f,    solsw_m,  solswfdiff_m, sollw_m,       &
+       rain_f,    snow_f,    bs_f, solsw_m,  solswfdiff_m, sollw_m,       &
        gustiness,                                     &
-       t,         q,         u,        v,             &
+       t,         q,        qbs,  u,        v,             &
 !!! nrlmd+jyg le 02/05/2011 et le 20/02/2012
 !!       t_x,       q_x,       t_w,      q_w,           &
@@ -275,8 +276,8 @@
        beta, &
 !>jyg
-       alb_dir_m,    alb_dif_m,  zxsens,   zxevap,    &
+       alb_dir_m,    alb_dif_m,  zxsens,   zxevap,  zxsnowerosion,  &
        alb3_lic,  runoff,    snowhgt,   qsnow,     to_ice,    sissnow,  &
        zxtsol,    zxfluxlat, zt2m,     qsat2m, zn2mout, &
-       d_t,       d_q,       d_u,      d_v, d_t_diss, & 
+       d_t,       d_q,    d_qbs,    d_u,      d_v, d_t_diss, & 
 !!! nrlmd+jyg le 02/05/2011 et le 20/02/2012
        d_t_w,     d_q_w,                              &
@@ -301,5 +302,5 @@
        rh2m,      zxfluxu,  zxfluxv,               &
        z0m, z0h,   agesno,  sollw,    solsw,         &
-       d_ts,      evap,    fluxlat,  t2m,           &
+       d_ts,      evap,    fluxlat,   t2m,           &
        wfbils,    wfbilo, wfevap, wfrain, wfsnow,   & 
        flux_t,   flux_u, flux_v,                    &
@@ -307,5 +308,5 @@
 !jyg<
 !!       zxfluxt,   zxfluxq,   q2m,      flux_q, tke,   &
-       zxfluxt,   zxfluxq,   q2m,      flux_q, tke_x,  &
+       zxfluxt,   zxfluxq, zxfluxqbs,   q2m, flux_q, flux_qbs, tke_x,  &
 !>jyg
 !!! nrlmd+jyg le 02/05/2011 et le 20/02/2012
@@ -410,4 +411,5 @@
          dser, dt_ds, zsig, zmea
     use phys_output_var_mod, only: tkt, tks, taur, sss
+    use blowing_snow_ini_mod, only : zeta_bs
 #ifdef CPP_XIOS
     USE wxios, ONLY: missing_val
@@ -444,4 +446,5 @@
     REAL, DIMENSION(klon),        INTENT(IN)        :: rain_f  ! rain fall
     REAL, DIMENSION(klon),        INTENT(IN)        :: snow_f  ! snow fall
+    REAL, DIMENSION(klon),        INTENT(IN)        :: bs_f  ! blowing snow fall
     REAL, DIMENSION(klon),        INTENT(IN)        :: solsw_m ! net shortwave radiation at mean surface
     REAL, DIMENSION(klon),        INTENT(IN)        :: solswfdiff_m ! diffuse fraction fordownward shortwave radiation at mean surface
@@ -449,4 +452,5 @@
     REAL, DIMENSION(klon,klev),   INTENT(IN)        :: t       ! temperature (K)
     REAL, DIMENSION(klon,klev),   INTENT(IN)        :: q       ! water vapour (kg/kg)
+    REAL, DIMENSION(klon,klev),   INTENT(IN)        :: qbs       ! blowing snow specific content (kg/kg)
     REAL, DIMENSION(klon,klev),   INTENT(IN)        :: u       ! u speed
     REAL, DIMENSION(klon,klev),   INTENT(IN)        :: v       ! v speed
@@ -521,4 +525,5 @@
                                                                   ! (=> positive sign upwards)
     REAL, DIMENSION(klon),        INTENT(OUT)       :: zxevap     ! water vapour flux at surface, positiv upwards
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: zxsnowerosion     ! blowing snow flux at surface
     REAL, DIMENSION(klon),        INTENT(OUT)       :: zxtsol     ! temperature at surface, mean for each grid point
 !!! jyg le ???
@@ -537,4 +542,6 @@
     REAL, DIMENSION(klon, klev),  INTENT(OUT)       :: d_u        ! change in u speed
     REAL, DIMENSION(klon, klev),  INTENT(OUT)       :: d_v        ! change in v speed
+    REAL, DIMENSION(klon, klev),  INTENT(OUT)       :: d_qbs        ! change in blowing snow specific content
+
 
     REAL, INTENT(OUT):: zcoefh(:, :, :) ! (klon, klev, nbsrf + 1)
@@ -604,5 +611,5 @@
     REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: sollw      ! net longwave radiation at surface
     REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: d_ts       ! change in temperature at surface
-    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)     :: evap     ! evaporation at surface
+    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)     :: evap       ! evaporation at surface
     REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: fluxlat    ! latent flux
     REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: t2m        ! temperature at 2 meter height
@@ -631,6 +638,9 @@
     REAL, DIMENSION(klon, klev), INTENT(OUT)        :: zxfluxt    ! sensible heat flux, mean for each grid point
     REAL, DIMENSION(klon, klev), INTENT(OUT)        :: zxfluxq    ! water vapour flux, mean for each grid point
+    REAL, DIMENSION(klon, klev), INTENT(OUT)        :: zxfluxqbs    ! blowing snow flux, mean for each grid point
     REAL, DIMENSION(klon, nbsrf),INTENT(OUT)        :: q2m        ! water vapour at 2 meter height
     REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_q     ! water vapour flux(latent flux) (kg/m**2/s)
+    REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_qbs   ! blowind snow vertical flux (kg/m**2
+
 #ifdef ISO   
     REAL, DIMENSION(ntraciso,klon),       INTENT(OUT)        :: dflux_xt    ! change of water vapour flux
@@ -683,8 +693,8 @@
     REAL, DIMENSION(klon)              :: yalb,yalb_vis
 !albedo SB <<<
-    REAL, DIMENSION(klon)              :: yt1, yq1, yu1, yv1
+    REAL, DIMENSION(klon)              :: yt1, yq1, yu1, yv1, yqbs1
     REAL, DIMENSION(klon)              :: yqa
     REAL, DIMENSION(klon)              :: ysnow, yqsurf, yagesno, yqsol
-    REAL, DIMENSION(klon)              :: yrain_f, ysnow_f
+    REAL, DIMENSION(klon)              :: yrain_f, ysnow_f, ybs_f
 #ifdef ISO
     REAL, DIMENSION(ntraciso,klon)              :: yxt1
@@ -699,4 +709,5 @@
     REAL, DIMENSION(klon)              :: yrugoro
     REAL, DIMENSION(klon)              :: yfluxlat
+    REAL, DIMENSION(klon)              :: yfluxbs
     REAL, DIMENSION(klon)              :: y_d_ts
     REAL, DIMENSION(klon)              :: y_flux_t1, y_flux_q1
@@ -707,4 +718,5 @@
 #endif
     REAL, DIMENSION(klon)              :: y_flux_u1, y_flux_v1
+    REAL, DIMENSION(klon)              :: y_flux_bs, y_flux0
     REAL, DIMENSION(klon)              :: yt2m, yq2m, yu10m
     INTEGER, DIMENSION(klon, nbsrf, 6) :: yn2mout, yn2mout_x, yn2mout_w
@@ -736,4 +748,5 @@
 #endif
     REAL, DIMENSION(klon)              :: AcoefU, AcoefV, BcoefU, BcoefV
+    REAL, DIMENSION(klon)              :: AcoefQBS, BcoefQBS
     REAL, DIMENSION(klon)              :: ypsref
     REAL, DIMENSION(klon)              :: yevap, yevap_pot, ytsurf_new, yalb3_new
@@ -744,12 +757,12 @@
     REAL, DIMENSION(klon)              :: meansqT ! mean square deviation of subsurface temperatures
     REAL, DIMENSION(klon)              :: alb_m  ! mean albedo for whole SW interval
-    REAL, DIMENSION(klon,klev)         :: y_d_t, y_d_q, y_d_t_diss
+    REAL, DIMENSION(klon,klev)         :: y_d_t, y_d_q, y_d_t_diss, y_d_qbs
     REAL, DIMENSION(klon,klev)         :: y_d_u, y_d_v
-    REAL, DIMENSION(klon,klev)         :: y_flux_t, y_flux_q
+    REAL, DIMENSION(klon,klev)         :: y_flux_t, y_flux_q, y_flux_qbs
     REAL, DIMENSION(klon,klev)         :: y_flux_u, y_flux_v
-    REAL, DIMENSION(klon,klev)         :: ycoefh, ycoefm,ycoefq
+    REAL, DIMENSION(klon,klev)         :: ycoefh, ycoefm,ycoefq, ycoefqbs
     REAL, DIMENSION(klon)              :: ycdragh, ycdragq, ycdragm
     REAL, DIMENSION(klon,klev)         :: yu, yv
-    REAL, DIMENSION(klon,klev)         :: yt, yq
+    REAL, DIMENSION(klon,klev)         :: yt, yq, yqbs
 #ifdef ISO
     REAL, DIMENSION(ntraciso,klon)              :: yxtevap
@@ -819,4 +832,5 @@
     REAL, DIMENSION(klon,klev)         :: CcoefH, CcoefQ, DcoefH, DcoefQ
     REAL, DIMENSION(klon,klev)         :: CcoefU, CcoefV, DcoefU, DcoefV
+    REAL, DIMENSION(klon,klev)         :: CcoefQBS, DcoefQBS
     REAL, DIMENSION(klon,klev)         :: CcoefH_x, CcoefQ_x, DcoefH_x, DcoefQ_x
     REAL, DIMENSION(klon,klev)         :: CcoefH_w, CcoefQ_w, DcoefH_w, DcoefQ_w
@@ -824,4 +838,5 @@
     REAL, DIMENSION(klon,klev)         :: CcoefU_w, CcoefV_w, DcoefU_w, DcoefV_w
     REAL, DIMENSION(klon,klev)         :: Kcoef_hq, Kcoef_m, gama_h, gama_q
+    REAL, DIMENSION(klon,klev)         :: gama_qbs, Kcoef_qbs 
     REAL, DIMENSION(klon,klev)         :: Kcoef_hq_x, Kcoef_m_x, gama_h_x, gama_q_x
     REAL, DIMENSION(klon,klev)         :: Kcoef_hq_w, Kcoef_m_w, gama_h_w, gama_q_w
@@ -1016,4 +1031,5 @@
     REAL, DIMENSION(klon,nbsrf)        :: zx_t1
     REAL, DIMENSION(klon, nbsrf)       :: alb          ! mean albedo for whole SW interval
+    REAL, DIMENSION(klon,nbsrf)        :: snowerosion    
     REAL, DIMENSION(klon)              :: ylwdown      ! jg : temporary (ysollwdown)
     REAL, DIMENSION(klon)              :: ygustiness      ! jg : temporary (ysollwdown)
@@ -1189,10 +1205,10 @@
   alb_dir_m=0. ; alb_dif_m=0. ; alb3_lic(:)=0.
 !albedo SB <<<
- zxsens(:)=0. ; zxevap(:)=0. ; zxtsol(:)=0.
+ zxsens(:)=0. ; zxevap(:)=0. ; zxtsol(:)=0. ; zxsnowerosion(:)=0.
  d_t_w(:,:)=0. ; d_q_w(:,:)=0. ; d_t_x(:,:)=0. ; d_q_x(:,:)=0.
  zxfluxlat(:)=0.
  zt2m(:)=0. ; zq2m(:)=0. ; qsat2m(:)=0. ; rh2m(:)=0.
  zn2mout(:,:)=0 ;
- d_t(:,:)=0. ; d_t_diss(:,:)=0. ; d_q(:,:)=0. ; d_u(:,:)=0. ; d_v(:,:)=0.
+ d_t(:,:)=0. ; d_t_diss(:,:)=0. ; d_q(:,:)=0. ; d_qbs(:,:)=0. ; d_u(:,:)=0. ; d_v(:,:)=0.
  zcoefh(:,:,:)=0. ; zcoefm(:,:,:)=0.
  zxsens_x(:)=0. ; zxsens_w(:)=0. ; zxfluxlat_x(:)=0. ; zxfluxlat_w(:)=0.
@@ -1214,11 +1230,13 @@
  d_ts(:,:)=0.
  evap(:,:)=0.
+ snowerosion(:,:)=0. 
  fluxlat(:,:)=0.
  wfbils(:,:)=0. ; wfbilo(:,:)=0.
  wfevap(:,:)=0. ; wfrain(:,:)=0. ; wfsnow(:,:)=0.
  flux_t(:,:,:)=0. ; flux_q(:,:,:)=0. ; flux_u(:,:,:)=0. ; flux_v(:,:,:)=0.
+ flux_qbs(:,:,:)=0.
  dflux_t(:)=0. ; dflux_q(:)=0.
  zxsnow(:)=0.
- zxfluxt(:,:)=0. ; zxfluxq(:,:)=0.
+ zxfluxt(:,:)=0. ; zxfluxq(:,:)=0.; zxfluxqbs(:,:)=0.
  qsnow(:)=0. ; snowhgt(:)=0. ; to_ice(:)=0. ; sissnow(:)=0.
  runoff(:)=0.
@@ -1266,12 +1284,13 @@
     yqsurf = 0.0  ; yalb = 0.0 ; yalb_vis = 0.0
 !albedo SB <<<
-    yrain_f = 0.0 ; ysnow_f = 0.0    ; yfder = 0.0     ; ysolsw = 0.0
+    yrain_f = 0.0 ; ysnow_f = 0.0  ; ybs_f=0.0  ; yfder = 0.0     ; ysolsw = 0.0
     ysollw = 0.0  ; yz0m = 0.0 ; yz0h = 0.0    ; yu1 = 0.0    
-    yv1 = 0.0     ; ypaprs = 0.0     ; ypplay = 0.0
+    yv1 = 0.0     ; ypaprs = 0.0     ; ypplay = 0.0     ; yqbs1 = 0.0
     ydelp = 0.0   ; yu = 0.0         ; yv = 0.0        ; yt = 0.0         
     yq = 0.0      ; y_dflux_t = 0.0  ; y_dflux_q = 0.0 
+    yqbs(:,:)=0.0  
     yrugoro = 0.0 ; ywindsp = 0.0   
 !!    d_ts = 0.0    ; yfluxlat=0.0     ; flux_t = 0.0    ; flux_q = 0.0     
-    yfluxlat=0.0 
+    yfluxlat=0.0 ; y_flux0(:)=0.0
 !!    flux_u = 0.0  ; flux_v = 0.0     ; d_t = 0.0       ; d_q = 0.0      
 !!    d_t_diss= 0.0 ;d_u = 0.0     ; d_v = 0.0 
@@ -1288,4 +1307,5 @@
     ycldt = 0.0      ; yrmu0 = 0.0
     ! Martin
+    y_d_qbs(:,:)=0.0
 
 !!! nrlmd+jyg le 02/05/2011 et le 20/02/2012
@@ -1629,4 +1649,5 @@
           yrain_f(j) = rain_f(i)
           ysnow_f(j) = snow_f(i)
+          ybs_f(j)   = bs_f(i)
           yagesno(j) = agesno(i,nsrf)
           yfder(j)   = fder(i)
@@ -1640,4 +1661,5 @@
           yu1(j)     = u(i,1)
           yv1(j)     = v(i,1)
+          yqbs1(j)   = qbs(i,1)
           ypaprs(j,klev+1) = paprs(i,klev+1)
 !jyg<
@@ -1653,4 +1675,6 @@
 !!! nrlmd le 13/06/2011
           y_delta_tsurf(j)=delta_tsurf(i,nsrf)
+          yfluxbs(j)=0.0
+          y_flux_bs(j) = 0.0
 !!!
 #ifdef ISO
@@ -1721,4 +1745,5 @@
              yt(j,k) = t(i,k)
              yq(j,k) = q(i,k)
+             yqbs(j,k)=qbs(i,k)
 #ifdef ISO
              do ixt=1,ntraciso   
@@ -1951,5 +1976,4 @@
       print *,' args coef_diff_turb: ycdragh ', ycdragh
       print *,' args coef_diff_turb: ytke ', ytke
-
        ENDIF
 
@@ -1960,5 +1984,4 @@
 
         ELSE
-
 
         CALL coef_diff_turb(dtime, nsrf, knon, ni,  &
@@ -1981,6 +2004,8 @@
 
         IF (prt_level >=10) print *,'coef_diff_turb -> ycoefh ',ycoefh
-!
+
+
        ELSE  !(iflag_split .eq.0)
+
       IF (prt_level >=10) THEN
       print *,' args coef_diff_turb: yu_x ',  yu_x  
@@ -2020,5 +2045,5 @@
        ENDIF
 
-       ENDIF ! iflag_pbl >= 50
+        ENDIF ! iflag_pbl >= 50
 
         IF (prt_level >=10) print *,'coef_diff_turb -> ycoefh_x ',ycoefh_x
@@ -2034,5 +2059,5 @@
       print *,' args coef_diff_turb: ycdragh_w ', ycdragh_w
       print *,' args coef_diff_turb: ytke_w ', ytke_w
-       ENDIF
+      ENDIF
 
         IF (iflag_pbl>=50) THEN
@@ -2062,5 +2087,5 @@
 
         IF (prt_level >=10) print *,'coef_diff_turb -> ycoefh_w ',ycoefh_w
-!
+
 !!!jyg le 10/04/2013
 !!   En attendant de traiter le transport des traceurs dans les poches froides, formule
@@ -2072,7 +2097,5 @@
         ENDDO
       ENDDO
-!!!
        ENDIF  ! (iflag_split .eq.0)
-!!!
        
 !****************************************************************************************
@@ -2171,4 +2194,24 @@
        ENDIF  ! (iflag_split .eq.0)
 !!!
+
+ ! For blowing snow:
+    IF (ok_bs) THEN
+     ! following Bintanja et al 2000, part II
+     ! we assume that the eddy diffsivity coefficient for
+     ! suspended particles is larger than Km by a factor zeta_bs
+     ! which is equal to 3 by default
+     do k=1,klev
+        do j=1,knon
+           ycoefqbs(j,k)=ycoefm(j,k)*zeta_bs 
+        enddo
+     enddo
+     CALL climb_qbs_down(knon, ycoefqbs, ypaprs, ypplay, &
+     ydelp, yt, yqbs, dtime, & 
+     CcoefQBS, DcoefQBS, &
+     Kcoef_qbs, gama_qbs, &
+     AcoefQBS, BcoefQBS)
+    ENDIF
+
+
 
 !****************************************************************************************
@@ -2356,5 +2399,5 @@
                debut, lafin, ydelp(:,1), r_co2_ppm, ysolsw, ysollw, yalb, &
 !!jyg               yts, ypplay(:,1), ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
-               yts, ypplay(:,1), ycdragh, ycdragm, yrain_f, ysnow_f, yt1, yq1,&
+               yts, ypplay(:,1), ycdragh, ycdragm, yrain_f, ysnow_f, ybs_f, yt1, yq1,&
                AcoefH, AcoefQ, BcoefH, BcoefQ, & 
                AcoefU, AcoefV, BcoefU, BcoefV, & 
@@ -2362,5 +2405,5 @@
                ylwdown, yq2m, yt2m, &
                ysnow, yqsol, yagesno, ytsoil, &
-               yz0m, yz0h, SFRWL, yalb_dir_new, yalb_dif_new, yevap, yfluxsens,yfluxlat,&
+               yz0m, yz0h, SFRWL, yalb_dir_new, yalb_dif_new, yevap, yfluxsens,yfluxlat,yfluxbs,&
                yqsurf, ytsurf_new, y_dflux_t, y_dflux_q, &
                y_flux_u1, y_flux_v1, &
@@ -2433,12 +2476,12 @@
                   yrmu0, ylwdown, yalb, zgeo1, &
                   ysolsw, ysollw, yts, ypplay(:,1), &
-                  !!jyg               ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
-                  ycdragh, ycdragm, yrain_f, ysnow_f, yt1, yq1,&
+                  ycdragh, ycdragm, yrain_f, ysnow_f, ybs_f, yt1, yq1,&
                   AcoefH, AcoefQ, BcoefH, BcoefQ, &
                   AcoefU, AcoefV, BcoefU, BcoefV, &
+                  AcoefQBS, BcoefQBS, &
                   ypsref, yu1, yv1, ygustiness, yrugoro, pctsrf, &
-                  ysnow, yqsurf, yqsol, yagesno, &
+                  ysnow, yqsurf, yqsol,yqbs1, yagesno, &
                   ytsoil, yz0m, yz0h, SFRWL, yalb_dir_new, yalb_dif_new, yevap,yfluxsens,yfluxlat, &
-                  ytsurf_new, y_dflux_t, y_dflux_q, &
+                  yfluxbs, ytsurf_new, y_dflux_t, y_dflux_q, &
                   yzmea, yzsig, ycldt, &
                   ysnowhgt, yqsnow, ytoice, ysissnow, &
@@ -2499,5 +2542,5 @@
                itap, dtime, jour, knon, ni, &
 !!jyg               ypplay(:,1), zgeo1/RG, ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
-               ypplay(:,1), zgeo1(1:knon)/RG, ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&    ! ym missing init
+               ypplay(:,1), zgeo1(1:knon)/RG, ycdragh, ycdragm, yrain_f, ysnow_f, ybs_f, yt(:,1), yq(:,1),&    ! ym missing init
                AcoefH, AcoefQ, BcoefH, BcoefQ, &
                AcoefU, AcoefV, BcoefU, BcoefV, &
@@ -2664,4 +2707,11 @@
           ENDDO
         ENDIF ! (ok_flux_surf)
+
+        ! flux of blowing snow at the first level 
+        IF (ok_bs) THEN 
+        DO j=1,knon
+        y_flux_bs(j)=yfluxbs(j)
+        ENDDO
+        ENDIF
 !
 ! ------------------------------------------------------------------------------
@@ -2991,4 +3041,14 @@
 !
        ENDIF  ! (iflag_split .eq.0)
+
+       IF (ok_bs) THEN
+            CALL climb_qbs_up(knon, dtime, yqbs, &
+            y_flux_bs, ypaprs, ypplay, &
+            AcoefQBS, BcoefQBS, &
+            CcoefQBS, DcoefQBS, &
+            Kcoef_qbs, gama_qbs, &
+            y_flux_qbs(:,:), y_d_qbs(:,:))
+       ENDIF
+
 !!!
 !!
@@ -3157,12 +3217,20 @@
 !!!
 
-!      print*,'Dans pbl OK1'
-
-!jyg<
-!!       evap(:,nsrf) = - flux_q(:,1,nsrf)
-!>jyg
+       ! tendencies of blowing snow 
+       IF (ok_bs) THEN 
+           DO k = 1, klev   
+            DO j = 1, knon 
+                i = ni(j)
+                y_d_qbs(j,k)=y_d_qbs(j,k) * ypct(j)
+                flux_qbs(i,k,nsrf) = y_flux_qbs(j,k) 
+            ENDDO
+          ENDDO
+       ENDIF
+
+
        DO j = 1, knon
           i = ni(j)
           evap(i,nsrf) = - flux_q(i,1,nsrf)                  !jyg
+          if (ok_bs) then ; snowerosion(i,nsrf)=flux_qbs(i,1,nsrf); endif
           beta(i,nsrf) = ybeta(j)                             !jyg
           d_ts(i,nsrf) = y_d_ts(j)
@@ -3392,4 +3460,15 @@
           ENDDO
        ENDDO
+
+
+       IF (ok_bs) THEN
+         DO k = 1, klev
+         DO j = 1, knon
+         i = ni(j)
+         d_qbs(i,k) = d_qbs(i,k) + y_d_qbs(j,k)
+         ENDDO
+         ENDDO
+        ENDIF
+
 
 #ifdef ISO
@@ -3894,4 +3973,22 @@
        fder_print(i) = fder(i) + dflux_t(i) + dflux_q(i)
     ENDDO
+
+    ! if blowing snow
+    if (ok_bs) then  
+       DO nsrf = 1, nbsrf 
+       DO k = 1, klev
+       DO i = 1, klon 
+         zxfluxqbs(i,k) = zxfluxqbs(i,k) + flux_qbs(i,k,nsrf) * pctsrf(i,nsrf)
+       ENDDO
+       ENDDO
+       ENDDO
+
+       DO i = 1, klon
+        zxsnowerosion(i)     = zxfluxqbs(i,1) ! blowings snow flux at the surface 
+       END DO
+    endif
+
+   
+    
 #ifdef ISO
     DO i = 1, klon
Index: LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmdiso/phyetat0_mod.F90
===================================================================
--- LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmdiso/phyetat0_mod.F90	(revision 4660)
+++ LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmdiso/phyetat0_mod.F90	(revision 4669)
@@ -24,7 +24,7 @@
        qsol, fevap, z0m, z0h, agesno, &
        du_gwd_rando, du_gwd_front, entr_therm, f0, fm_therm, &
-       falb_dir, falb_dif, prw_ancien, prlw_ancien, prsw_ancien, &
-       ftsol, pbl_tke, pctsrf, q_ancien, ql_ancien, qs_ancien, rneb_ancien, radpas, radsol, rain_fall, ratqs, &
-       rnebcon, rugoro, sig1, snow_fall, solaire_etat0, sollw, sollwdown, &
+       falb_dir, falb_dif, prw_ancien, prlw_ancien, prsw_ancien, prbsw_ancien, &
+       ftsol, pbl_tke, pctsrf, q_ancien, ql_ancien, qs_ancien, qbs_ancien, rneb_ancien, radpas, radsol, rain_fall, ratqs, &
+       rnebcon, rugoro, sig1, snow_fall, bs_fall, solaire_etat0, sollw, sollwdown, &
        solsw, solswfdiff, t_ancien, u_ancien, v_ancien, w01, wake_cstar, wake_deltaq, &
        wake_deltat, wake_delta_pbl_TKE, delta_tsurf, beta_aridity, wake_fip, wake_pe, &
@@ -349,4 +349,11 @@
   found=phyetat0_get(rain_fall,"rain_f","rain fall",0.)
 
+  IF (ok_bs) THEN
+     found=phyetat0_get(bs_fall,"bs_f","blowing snow fall",0.)
+  ELSE
+     bs_fall(:)=0.
+  ENDIF
+
+
 !=======================================================================
 ! Radiation
@@ -410,4 +417,14 @@
   ancien_ok=ancien_ok.AND.phyetat0_get(prsw_ancien,"PRSWANCIEN","PRSWANCIEN",0.)
 
+  ! cas specifique des variables de la neige soufflee
+  IF (ok_bs) THEN
+     ancien_ok=ancien_ok.AND.phyetat0_get(qbs_ancien,"QBSANCIEN","QBSANCIEN",0.)
+     ancien_ok=ancien_ok.AND.phyetat0_get(prbsw_ancien,"PRBSWANCIEN","PRBSWANCIEN",0.)
+  ELSE
+     qbs_ancien(:,:)=0.
+     prbsw_ancien(:)=0.
+  ENDIF
+
+
   ! Ehouarn: addtional tests to check if t_ancien, q_ancien contain
   !          dummy values (as is the case when generated by ce0l,
@@ -423,4 +440,12 @@
     ancien_ok=.false.
   ENDIF
+
+  IF (ok_bs) THEN
+    IF ( (maxval(qbs_ancien).EQ.minval(qbs_ancien))       .OR. &
+         (maxval(prbsw_ancien).EQ.minval(prbsw_ancien)) ) THEN
+       ancien_ok=.false.
+    ENDIF
+  ENDIF
+
 
   found=phyetat0_get(clwcon,"CLWCON","CLWCON",0.)
Index: LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmdiso/phyredem.F90
===================================================================
--- LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmdiso/phyredem.F90	(revision 4660)
+++ LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmdiso/phyredem.F90	(revision 4669)
@@ -15,9 +15,9 @@
                                 ftsol, beta_aridity, delta_tsurf, falb_dir,  &
                                 falb_dif, qsol, fevap, radsol, solsw, sollw, &
-                                sollwdown, rain_fall, snow_fall, z0m, z0h,   &
+                                sollwdown, rain_fall, snow_fall, bs_fall, z0m, z0h, &
                                 agesno, zmea, zstd, zsig, zgam, zthe, zpic,  &
                                 zval, rugoro, t_ancien, q_ancien,            &
-                                prw_ancien, prlw_ancien, prsw_ancien,        &
-                                ql_ancien, qs_ancien,  u_ancien,             &
+                                prw_ancien, prlw_ancien, prsw_ancien, prbsw_ancien,      &
+                                ql_ancien, qs_ancien, qbs_ancien,  u_ancien, &
                                 v_ancien, clwcon, rnebcon, ratqs, pbl_tke,   &
                                 wake_delta_pbl_tke, zmax0, f0, sig1, w01,    &
@@ -260,4 +260,10 @@
 
     CALL put_field(pass,"QSANCIEN", "QSANCIEN", qs_ancien)
+
+    IF (ok_bs) THEN
+       CALL put_field(pass,"bs_f", "precipitation neige soufflee", bs_fall)
+       CALL put_field(pass,"QBSANCIEN", "QBSANCIEN", qbs_ancien)
+       CALL put_field(pass,"PRBSWANCIEN", "PRBSWANCIEN", prbsw_ancien)
+    ENDIF
 
     CALL put_field(pass,"PRWANCIEN", "PRWANCIEN", prw_ancien)
Index: LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmdiso/phys_local_var_mod.F90
===================================================================
--- LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmdiso/phys_local_var_mod.F90	(revision 4660)
+++ LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmdiso/phys_local_var_mod.F90	(revision 4669)
@@ -14,4 +14,6 @@
       REAL, SAVE, ALLOCATABLE :: ql_seri(:,:),qs_seri(:,:)
       !$OMP THREADPRIVATE(ql_seri,qs_seri)
+      REAL, SAVE, ALLOCATABLE :: qbs_seri(:,:)
+      !$OMP THREADPRIVATE(qbs_seri)
       REAL, SAVE, ALLOCATABLE :: u_seri(:,:), v_seri(:,:)
       !$OMP THREADPRIVATE(u_seri, v_seri)
@@ -26,8 +28,8 @@
       REAL, SAVE, ALLOCATABLE :: d_t_dyn(:,:), d_q_dyn(:,:)
       !$OMP THREADPRIVATE(d_t_dyn, d_q_dyn)
-      REAL, SAVE, ALLOCATABLE :: d_ql_dyn(:,:), d_qs_dyn(:,:)
-      !$OMP THREADPRIVATE(d_ql_dyn, d_qs_dyn)
-      REAL, SAVE, ALLOCATABLE :: d_q_dyn2d(:), d_ql_dyn2d(:), d_qs_dyn2d(:)
-      !$OMP THREADPRIVATE(d_q_dyn2d, d_ql_dyn2d, d_qs_dyn2d)
+      REAL, SAVE, ALLOCATABLE :: d_ql_dyn(:,:), d_qs_dyn(:,:), d_qbs_dyn(:,:)
+      !$OMP THREADPRIVATE(d_ql_dyn, d_qs_dyn, d_qbs_dyn)
+      REAL, SAVE, ALLOCATABLE :: d_q_dyn2d(:), d_ql_dyn2d(:), d_qs_dyn2d(:), d_qbs_dyn2d(:)
+      !$OMP THREADPRIVATE(d_q_dyn2d, d_ql_dyn2d, d_qs_dyn2d, d_qbs_dyn2d)
       REAL, SAVE, ALLOCATABLE :: d_u_dyn(:,:), d_v_dyn(:,:)
       !$OMP THREADPRIVATE(d_u_dyn, d_v_dyn)
@@ -69,6 +71,6 @@
       REAL, SAVE, ALLOCATABLE :: d_u_oli(:,:), d_v_oli(:,:)
       !$OMP THREADPRIVATE(d_u_oli, d_v_oli)
-      REAL, SAVE, ALLOCATABLE :: d_t_vdf(:,:), d_q_vdf(:,:), d_t_diss(:,:)
-      !$OMP THREADPRIVATE( d_t_vdf, d_q_vdf,d_t_diss)
+      REAL, SAVE, ALLOCATABLE :: d_t_vdf(:,:), d_q_vdf(:,:), d_qbs_vdf(:,:), d_t_diss(:,:)
+      !$OMP THREADPRIVATE( d_t_vdf, d_q_vdf, d_qbs_vdf, d_t_diss)
       REAL, SAVE, ALLOCATABLE :: d_u_vdf(:,:), d_v_vdf(:,:)
       !$OMP THREADPRIVATE(d_u_vdf, d_v_vdf)
@@ -78,4 +80,6 @@
       REAL, SAVE, ALLOCATABLE :: d_t_vdf_x(:,:), d_q_vdf_x(:,:)
       !$OMP THREADPRIVATE( d_t_vdf_x, d_q_vdf_x)
+      REAL, SAVE, ALLOCATABLE :: d_t_bs(:,:), d_q_bs(:,:), d_qbs_bs(:,:)
+      !$OMP THREADPRIVATE( d_t_bs,d_q_bs, d_qbs_bs)
 !>nrlmd+jyg
       REAL, SAVE, ALLOCATABLE :: d_t_oro(:,:)
@@ -313,4 +317,8 @@
 !$OMP THREADPRIVATE(sens, flwp, fiwp)
 !!
+!FC
+      REAL,ALLOCATABLE,SAVE,DIMENSION(:,:) :: zxfluxt, zxfluxq
+!$OMP THREADPRIVATE(zxfluxt, zxfluxq)
+!FC
 !!         Wake variables
       REAL,ALLOCATABLE,SAVE,DIMENSION(:)            :: alp_wake
@@ -361,6 +369,6 @@
       REAL,ALLOCATABLE,SAVE,DIMENSION(:) :: JrNt
 !$OMP THREADPRIVATE(JrNt)
-      REAL,ALLOCATABLE,SAVE,DIMENSION(:) :: dthmin, evap, fder, plcl, plfc, prw, prlw, prsw
-!$OMP THREADPRIVATE(dthmin, evap, fder, plcl, plfc, prw, prlw, prsw)
+      REAL,ALLOCATABLE,SAVE,DIMENSION(:) :: dthmin, evap, snowerosion, fder, plcl, plfc, prw, prlw, prsw, prbsw
+!$OMP THREADPRIVATE(dthmin, evap, snowerosion, fder, plcl, plfc, prw, prlw, prsw, prbsw)
       REAL,ALLOCATABLE,SAVE,DIMENSION(:) :: zustar, zu10m, zv10m, rh2m
 !$OMP THREADPRIVATE(zustar, zu10m, zv10m, rh2m)
@@ -379,4 +387,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(:) :: zxfqcalving
 !$OMP THREADPRIVATE(zxfqcalving)
@@ -567,6 +577,6 @@
       REAL,ALLOCATABLE,SAVE,DIMENSION(:,:) :: zx_rh, zx_rhl, zx_rhi
 !$OMP THREADPRIVATE(zx_rh, zx_rhl, zx_rhi)
-      REAL,ALLOCATABLE,SAVE,DIMENSION(:,:) :: prfl, psfl, fraca
-!$OMP THREADPRIVATE(prfl, psfl, fraca)
+      REAL,ALLOCATABLE,SAVE,DIMENSION(:,:) :: prfl, psfl, fraca, bsfl
+!$OMP THREADPRIVATE(prfl, psfl, fraca, bsfl)
       REAL,ALLOCATABLE,SAVE,DIMENSION(:,:) :: Vprecip, zw2
 !$OMP THREADPRIVATE(Vprecip, zw2)
@@ -585,4 +595,8 @@
       REAL, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rneb,rnebjn,rneblsvol
 !$OMP THREADPRIVATE(rneb,rnebjn,rneblsvol)
+      REAL, ALLOCATABLE, SAVE, DIMENSION(:,:) :: pfraclr,pfracld
+!$OMP THREADPRIVATE(pfraclr,pfracld)
+      REAL, ALLOCATABLE, SAVE, DIMENSION(:,:) :: distcltop
+!$OMP THREADPRIVATE(distcltop)
 
 ! variables de sorties MM
@@ -734,5 +748,5 @@
 
 IMPLICIT NONE
-      ALLOCATE(t_seri(klon,klev),q_seri(klon,klev),ql_seri(klon,klev),qs_seri(klon,klev))
+      ALLOCATE(t_seri(klon,klev),q_seri(klon,klev),ql_seri(klon,klev),qs_seri(klon,klev), qbs_seri(klon,klev))
       ALLOCATE(u_seri(klon,klev),v_seri(klon,klev))
       ALLOCATE(l_mixmin(klon,klev+1,nbsrf),l_mix(klon,klev+1,nbsrf),tke_dissip(klon,klev+1,nbsrf),wprime(klon,klev+1,nbsrf))
@@ -741,6 +755,6 @@
       ALLOCATE(tr_seri(klon,klev,nbtr))
       ALLOCATE(d_t_dyn(klon,klev),d_q_dyn(klon,klev))
-      ALLOCATE(d_ql_dyn(klon,klev),d_qs_dyn(klon,klev))
-      ALLOCATE(d_q_dyn2d(klon),d_ql_dyn2d(klon),d_qs_dyn2d(klon))
+      ALLOCATE(d_ql_dyn(klon,klev),d_qs_dyn(klon,klev), d_qbs_dyn(klon,klev))
+      ALLOCATE(d_q_dyn2d(klon),d_ql_dyn2d(klon),d_qs_dyn2d(klon), d_qbs_dyn2d(klon))
       ALLOCATE(d_u_dyn(klon,klev),d_v_dyn(klon,klev))
       ALLOCATE(d_tr_dyn(klon,klev,nbtr))                   !RomP
@@ -765,5 +779,6 @@
       ALLOCATE(plul_st(klon),plul_th(klon))
       ALLOCATE(d_t_vdf(klon,klev),d_q_vdf(klon,klev),d_t_diss(klon,klev))
-
+      ALLOCATE (d_qbs_vdf(klon,klev))
+      ALLOCATE(d_t_bs(klon,klev),d_q_bs(klon,klev),d_qbs_bs(klon,klev))
       ALLOCATE(d_t_vdf_w(klon,klev),d_q_vdf_w(klon,klev))
       ALLOCATE(d_t_vdf_x(klon,klev),d_q_vdf_x(klon,klev))
@@ -925,6 +940,6 @@
       ALLOCATE(cldm(klon), cldq(klon), cldt(klon), qsat2m(klon))
       ALLOCATE(JrNt(klon))
-      ALLOCATE(dthmin(klon), evap(klon), fder(klon), plcl(klon), plfc(klon))
-      ALLOCATE(prw(klon), prlw(klon), prsw(klon), zustar(klon), zu10m(klon), zv10m(klon), rh2m(klon))
+      ALLOCATE(dthmin(klon), evap(klon), snowerosion(klon), fder(klon), plcl(klon), plfc(klon))
+      ALLOCATE(prw(klon), prlw(klon), prsw(klon), prbsw(klon), zustar(klon), zu10m(klon), zv10m(klon), rh2m(klon))
       ALLOCATE(s_lcl(klon))
       ALLOCATE(s_pblh(klon), s_pblt(klon), s_therm(klon))
@@ -1012,4 +1027,7 @@
       ALLOCATE(wfevap(klon, nbsrf), wfrain(klon,nbsrf), wfsnow(klon, nbsrf))
       ALLOCATE(evap_pot(klon, nbsrf))
+! FC
+      ALLOCATE(zxfluxq(klon,klev),zxfluxt(klon,klev))
+!
 !
 !  Deep convective variables used in phytrac
@@ -1055,4 +1073,5 @@
       ALLOCATE(prfl(klon, klev+1))
       ALLOCATE(psfl(klon, klev+1), fraca(klon, klev+1), Vprecip(klon, klev+1))
+      ALLOCATE(bsfl(klon,klev+1))
       ALLOCATE(zw2(klon, klev+1))
 
@@ -1068,5 +1087,7 @@
       ALLOCATE(beta_prec(klon,klev))
       ALLOCATE(rneb(klon,klev),rnebjn(klon,klev),rneblsvol(klon,klev))
-
+      ALLOCATE(pfraclr(klon,klev),pfracld(klon,klev))
+      pfraclr(:,:)=0. ; pfracld(:,:)=0. ! because not always defined
+      ALLOCATE(distcltop(klon,klev))
 
       ALLOCATE (zxsnow(klon),snowhgt(klon),qsnow(klon),to_ice(klon))
@@ -1141,5 +1162,5 @@
 USE indice_sol_mod
 IMPLICIT NONE
-      DEALLOCATE(t_seri,q_seri,ql_seri,qs_seri)
+      DEALLOCATE(t_seri,q_seri,ql_seri,qs_seri, qbs_seri)
       DEALLOCATE(u_seri,v_seri)
       DEALLOCATE(l_mixmin,l_mix, tke_dissip,wprime)
@@ -1147,6 +1168,6 @@
       DEALLOCATE(tr_seri)
       DEALLOCATE(d_t_dyn,d_q_dyn)
-      DEALLOCATE(d_ql_dyn,d_qs_dyn)
-      DEALLOCATE(d_q_dyn2d,d_ql_dyn2d,d_qs_dyn2d)
+      DEALLOCATE(d_ql_dyn,d_qs_dyn, d_qbs_dyn)
+      DEALLOCATE(d_q_dyn2d,d_ql_dyn2d,d_qs_dyn2d, d_qbs_dyn2d)
       DEALLOCATE(d_u_dyn,d_v_dyn)
       DEALLOCATE(d_tr_dyn)                      !RomP
@@ -1171,4 +1192,6 @@
       DEALLOCATE(plul_st,plul_th)
       DEALLOCATE(d_t_vdf,d_q_vdf,d_t_diss)
+      DEALLOCATE(d_qbs_vdf)
+      DEALLOCATE(d_t_bs,d_q_bs,d_qbs_bs)
 #ifdef ISO
       deallocate(xt_seri,xtl_seri,xts_seri)
@@ -1308,6 +1331,6 @@
       DEALLOCATE(cldm, cldq, cldt, qsat2m)
       DEALLOCATE(JrNt)
-      DEALLOCATE(dthmin, evap, fder, plcl, plfc)
-      DEALLOCATE(prw, prlw, prsw, zustar, zu10m, zv10m, rh2m, s_lcl)
+      DEALLOCATE(dthmin, evap, snowerosion, fder, plcl, plfc)
+      DEALLOCATE(prw, prlw, prsw, prbsw, zustar, zu10m, zv10m, rh2m, s_lcl)
       DEALLOCATE(s_pblh, s_pblt, s_therm)
 !
@@ -1322,4 +1345,5 @@
       DEALLOCATE(zxfqcalving, zxfluxlat)
       DEALLOCATE(zxrunofflic)
+      DEALLOCATE(zxustartlic, zxrhoslic)
       DEALLOCATE(zxtsol, snow_lsc, zxfqfonte, zxqsurf)
       DEALLOCATE(rain_lsc)
@@ -1364,4 +1388,6 @@
       DEALLOCATE(alp_bl_stat, n2, s2)
       DEALLOCATE(proba_notrig, random_notrig)
+!FC
+      DEALLOCATE(zxfluxq,zxfluxt)
 
       DEALLOCATE(dnwd0)
@@ -1423,5 +1449,5 @@
 
 
-      DEALLOCATE(prfl, psfl, fraca, Vprecip)
+      DEALLOCATE(prfl, psfl, bsfl, fraca, Vprecip)
       DEALLOCATE(zw2)
 
@@ -1436,4 +1462,6 @@
       DEALLOCATE(beta_prec)
       DEALLOCATE(rneb)
+      DEALLOCATE(pfraclr,pfracld)
+      DEALLOCATE(distcltop)
       DEALLOCATE (zxsnow,snowhgt,qsnow,to_ice,sissnow,runoff,albsol3_lic)
 #ifdef ISO
Index: LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmdiso/phys_output_ctrlout_mod.F90
===================================================================
--- LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmdiso/phys_output_ctrlout_mod.F90	(revision 4660)
+++ LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmdiso/phys_output_ctrlout_mod.F90	(revision 4669)
@@ -378,7 +378,14 @@
   TYPE(ctrl_out), SAVE :: o_snow = ctrl_out((/ 1, 1, 10, 10, 5, 10, 11, 11, 11, 11/), &
     'snow', 'Snow fall', 'kg/(s*m2)', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_bsfall = ctrl_out((/ 10, 10, 10, 10, 5, 10, 11, 11, 11, 11/), &
+    'bsfall', 'Blowing Snow fall', 'kg/(s*m2)', (/ ('', i=1, 10) /))
   TYPE(ctrl_out), SAVE :: o_evap = ctrl_out((/ 1, 1, 10, 10, 10, 10, 11, 11, 11, 11/), &
     'evap', 'Evaporat', 'kg/(s*m2)', (/ ('', i=1, 10) /))
-
+  TYPE(ctrl_out), SAVE :: o_snowerosion = ctrl_out((/ 10, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
+   'snowerosion', 'blowing snow flux', 'kg/(s*m2)', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_ustart_lic = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), & 
+    'ustart_lic', 'threshold velocity', 'm/s', (/ ('', i=1, 10) /))
+  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_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) /))
@@ -743,4 +750,6 @@
   TYPE(ctrl_out), SAVE :: o_prsw = ctrl_out((/ 1, 1, 10, 10, 10, 10, 11, 11, 11, 11/), &
     'prsw', 'Precipitable solid water', 'kg/m2', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_prbsw = ctrl_out((/ 10, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
+    'prbsw', 'Precipitable blowing snow', 'kg/m2', (/ ('', i=1, 10) /))
   TYPE(ctrl_out), SAVE :: o_s_pblh = ctrl_out((/ 1, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
     's_pblh', 'Boundary Layer Height', 'm', (/ ('', i=1, 10) /))
@@ -1413,7 +1422,7 @@
     'ec550aer', 'Extinction at 550nm', 'm^-1', (/ ('', i=1, 10) /))
   TYPE(ctrl_out), SAVE :: o_lwcon = ctrl_out((/ 2, 5, 10, 10, 10, 10, 11, 11, 11, 11/), &
-    'lwcon', 'Cloud liquid water content', 'kg/kg', (/ ('', i=1, 10) /))
+    'lwcon', 'Cloud liquid water content seen by radiation', 'kg/kg', (/ ('', i=1, 10) /))
   TYPE(ctrl_out), SAVE :: o_iwcon = ctrl_out((/ 2, 5, 10, 10, 10, 10, 11, 11, 11, 11/), &
-    'iwcon', 'Cloud ice water content', 'kg/kg', (/ ('', i=1, 10) /))
+    'iwcon', 'Cloud ice water content seen by radiation', 'kg/kg', (/ ('', i=1, 10) /))
   TYPE(ctrl_out), SAVE :: o_temp = ctrl_out((/ 2, 3, 4, 10, 10, 10, 11, 11, 11, 11/), &
     'temp', 'Air temperature', 'K', (/ ('', i=1, 10) /))
@@ -1432,4 +1441,6 @@
   TYPE(ctrl_out), SAVE :: o_ocond = ctrl_out((/ 2, 3, 4, 10, 10, 10, 11, 11, 11, 11/), &
     'ocond', 'Condensed water', 'kg/kg', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_qbs = ctrl_out((/ 10, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
+    'qbs', 'Specific content of blowing snow', 'kg/kg', (/ ('', i=1, 10) /))
   TYPE(ctrl_out), SAVE :: o_wvapp = ctrl_out((/ 2, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
     'wvapp', '', '', (/ ('', i=1, 10) /))
@@ -1462,4 +1473,8 @@
   TYPE(ctrl_out), SAVE :: o_rneblsvol = ctrl_out((/ 2, 5, 10, 10, 10, 10, 11, 11, 11, 11/), &
     'rneblsvol', 'LS Cloud fraction by volume', '-', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_pfraclr = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'pfraclr', 'LS precipitation fraction clear-sky part', '-', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_pfracld = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'pfracld', 'LS precipitation fraction cloudy part', '-', (/ ('', i=1, 10) /))
   TYPE(ctrl_out), SAVE :: o_rhum = ctrl_out((/ 2, 5, 10, 10, 10, 10, 11, 11, 11, 11/), &
     'rhum', 'Relative humidity', '-', (/ ('', i=1, 10) /))
@@ -1494,4 +1509,8 @@
   TYPE(ctrl_out), SAVE :: o_dqsphy2d = ctrl_out((/ 2, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
     'dqsphy2d', 'Physics dQS', '(kg/m2)/s', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_dqbsphy = ctrl_out((/ 10, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
+    'dqbsphy', 'Physics dQBS', '(kg/kg)/s', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_dqbsphy2d = ctrl_out((/ 10, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
+    'dqbsphy2d', 'Physics dQBS', '(kg/m2)/s', (/ ('', i=1, 10) /)) 
   TYPE(ctrl_out), SAVE :: o_pr_con_l = ctrl_out((/ 2, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
     'pr_con_l', 'Convective precipitation lic', ' ', (/ ('', i=1, 10) /))
@@ -1502,4 +1521,6 @@
   TYPE(ctrl_out), SAVE :: o_pr_lsc_i = ctrl_out((/ 2, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
     'pr_lsc_i', 'Large scale precipitation ice', ' ', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_pr_bs = ctrl_out((/ 10, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
+    'pr_bs', 'profile of blowing snow flux', ' ', (/ ('', i=1, 10) /))
   TYPE(ctrl_out), SAVE :: o_re = ctrl_out((/ 5, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
     're', 'Cloud droplet effective radius', 'um', (/ ('', i=1, 10) /))
@@ -1530,4 +1551,9 @@
   TYPE(ctrl_out), SAVE :: o_stratomask = ctrl_out((/ 2,  6, 10, 10, 10, 10, 11, 11, 11, 11/), &
     'stratomask', 'Stratospheric fraction', '1', (/ ('', i=1, 10) /))
+!FC
+  TYPE(ctrl_out), SAVE :: o_zxfluxt = ctrl_out((/ 2,  6, 10, 10, 10, 10, 11, 11, 11, 11/), &
+    'fluxt', 'flux h ', 'W/m2', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_zxfluxq = ctrl_out((/ 2,  6, 10, 10, 10, 10, 11, 11, 11, 11/), &
+    'fluxq', 'flux q ', 'kg/(s*m2)', (/ ('', i=1, 10) /))
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
@@ -1599,4 +1625,8 @@
   TYPE(ctrl_out), SAVE :: o_dqsdyn2d = ctrl_out((/ 4, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
     'dqsdyn2d', 'Dynamics dQS', '(kg/m2)/s', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_dqbsdyn = ctrl_out((/ 10, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
+    'dqbsdyn', 'Dynamics dQBS', '(kg/kg)/s', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_dqbsdyn2d = ctrl_out((/ 10, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
+    'dqbsdyn2d', 'Dynamics dQBS', '(kg/m2)/s', (/ ('', i=1, 10) /))
   TYPE(ctrl_out), SAVE :: o_dudyn = ctrl_out((/ 4, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
     'dudyn', 'Dynamics dU', 'm/s2', (/ ('', i=1, 10) /))
@@ -1673,4 +1703,12 @@
   TYPE(ctrl_out), SAVE :: o_dqeva2d = ctrl_out((/ 4, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
     'dqeva2d', 'Reevaporation dQ', '(kg/m2)/s', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_dqbsvdf = ctrl_out((/ 10, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
+    'dqbsvdf', 'Boundary-layer dQBS', '(kg/kg)/s', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_dqbsbs = ctrl_out((/ 10, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
+    'dqbsbs', 'Blowing snow dQBS', '(kg/kg)/s', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_dtbs = ctrl_out((/ 10, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
+    'dtbs', 'Blowing snow dT', '(K)/s', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_dqbs = ctrl_out((/ 10, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
+    'dqbs', 'Blowing snow dQ', '(kg/kg)/s', (/ ('', i=1, 10) /))
 
 !!!!!!!!!!!!!!!! Specifique thermiques
@@ -1835,4 +1873,5 @@
   TYPE(ctrl_out), SAVE, ALLOCATABLE :: o_xtprecip(:)
   TYPE(ctrl_out), SAVE, ALLOCATABLE :: o_xtevap(:)
+  TYPE(ctrl_out), SAVE, ALLOCATABLE :: o_xtevap_srf(:,:) ! ajout Camille 8 mai 2023
   TYPE(ctrl_out), SAVE, ALLOCATABLE :: o_xtplul(:)
   TYPE(ctrl_out), SAVE, ALLOCATABLE :: o_xtpluc(:)
Index: LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmdiso/phys_output_mod.F90
===================================================================
--- LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmdiso/phys_output_mod.F90	(revision 4660)
+++ LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmdiso/phys_output_mod.F90	(revision 4669)
@@ -178,4 +178,5 @@
     ALLOCATE(o_xtpluc(ntraciso))
     ALLOCATE(o_xtevap(ntraciso))
+    ALLOCATE(o_xtevap_srf(ntraciso,4))
     ALLOCATE(o_xtovap(ntraciso))
     ALLOCATE(o_xtoliq(ntraciso))
@@ -558,4 +559,11 @@
       o_xtevap  (ixt)=ctrl_out(flag,   'evap'//TRIM(outiso),             'Evaporat.', unit, [('',i=1,nfiles)])
 
+      ! ajout Camille 8 mai 2023
+      flag = [1, 6, 10, 10, 10, 10, 11, 11, 11, 11]
+      o_xtevap_srf (ixt,1)=ctrl_out(flag,   'evap_ter'//TRIM(outiso), 'Evap sfc'//clnsurf(1), unit, [('',i=1,nfiles)])
+      o_xtevap_srf (ixt,2)=ctrl_out(flag,   'evap_lic'//TRIM(outiso), 'Evap sfc'//clnsurf(2), unit, [('',i=1,nfiles)])
+      o_xtevap_srf (ixt,3)=ctrl_out(flag,   'evap_oce'//TRIM(outiso), 'Evap sfc'//clnsurf(3), unit, [('',i=1,nfiles)])
+      o_xtevap_srf (ixt,4)=ctrl_out(flag,   'evap_sic'//TRIM(outiso), 'Evap sfc'//clnsurf(4), unit, [('',i=1,nfiles)])
+
       flag = [2,  3,  4, 10, 10, 10, 11, 11, 11, 11]; unit = 'kg/kg'
       o_xtovap  (ixt)=ctrl_out(flag,   'ovap'//TRIM(outiso),     'Specific humidity', unit, [('',i=1,nfiles)])
Index: LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmdiso/phys_output_var_mod.F90
===================================================================
--- LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmdiso/phys_output_var_mod.F90	(revision 4660)
+++ LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmdiso/phys_output_var_mod.F90	(revision 4669)
@@ -29,5 +29,6 @@
   REAL, SAVE, ALLOCATABLE :: d_qw_col(:)      ! watter vapour mass budget for each column (kg/m2/s)
   REAL, SAVE, ALLOCATABLE :: d_ql_col(:)      ! liquid watter mass budget for each column (kg/m2/s)
-  REAL, SAVE, ALLOCATABLE :: d_qs_col(:)      ! solid watter mass budget for each column (kg/m2/s)
+  REAL, SAVE, ALLOCATABLE :: d_qs_col(:)      ! cloud ice mass budget for each column (kg/m2/s)
+  REAL, SAVE, ALLOCATABLE :: d_qbs_col(:)     ! blowing snow mass budget for each column (kg/m2/s)
   REAL, SAVE, ALLOCATABLE :: d_qt_col(:)      ! total watter mass budget for each column (kg/m2/s)
   REAL, SAVE, ALLOCATABLE :: d_ek_col(:)      ! kinetic energy budget for each column (W/m2)
@@ -35,8 +36,9 @@
   REAL, SAVE, ALLOCATABLE :: d_h_qw_col(:)    ! enthalpy budget of watter vapour for each column (W/m2)
   REAL, SAVE, ALLOCATABLE :: d_h_ql_col(:)    ! enthalpy budget of liquid watter for each column (W/m2)
-  REAL, SAVE, ALLOCATABLE :: d_h_qs_col(:)    ! enthalpy budget of solid watter  for each column (W/m2)
+  REAL, SAVE, ALLOCATABLE :: d_h_qs_col(:)    ! enthalpy budget of cloud ice  for each column (W/m2)
+  REAL, SAVE, ALLOCATABLE :: d_h_qbs_col(:)    ! enthalpy budget of blowing snow for each column (W/m2)
   REAL, SAVE, ALLOCATABLE :: d_h_col(:)       ! total enthalpy budget for each column (W/m2)
-  !$OMP THREADPRIVATE(d_qw_col, d_ql_col, d_qs_col, d_qt_col, d_ek_col, d_h_dair_col)
-  !$OMP THREADPRIVATE(d_h_qw_col, d_h_ql_col, d_h_qs_col, d_h_col)
+  !$OMP THREADPRIVATE(d_qw_col, d_ql_col, d_qs_col, d_qbs_col, d_qt_col, d_ek_col, d_h_dair_col)
+  !$OMP THREADPRIVATE(d_h_qw_col, d_h_ql_col, d_h_qs_col, d_h_qbs_col, d_h_col)
 
   ! Outputs used in cloudth_vert to extract the moments of the horizontal and 
@@ -173,8 +175,8 @@
 
     allocate (bils_ec(klon),bils_ech(klon),bils_tke(klon),bils_diss(klon),bils_kinetic(klon),bils_enthalp(klon),bils_latent(klon))
-    allocate (d_qw_col(klon), d_ql_col(klon), d_qs_col(klon), d_qt_col(klon), d_ek_col(klon), d_h_dair_col(klon) &
-  &         , d_h_qw_col(klon), d_h_ql_col(klon), d_h_qs_col(klon), d_h_col(klon))
-    d_qw_col=0. ; d_ql_col=0. ; d_qs_col=0. ; d_qt_col=0. ; d_ek_col=0. ; d_h_dair_col =0.
-    d_h_qw_col=0. ; d_h_ql_col=0. ; d_h_qs_col=0. ; d_h_col=0.
+    allocate (d_qw_col(klon), d_ql_col(klon), d_qs_col(klon), d_qbs_col(klon), d_qt_col(klon), d_ek_col(klon), d_h_dair_col(klon) &
+  &         , d_h_qw_col(klon), d_h_ql_col(klon), d_h_qs_col(klon), d_h_qbs_col(klon), d_h_col(klon))
+    d_qw_col=0. ; d_ql_col=0. ; d_qs_col=0. ; d_qbs_col=0. ; d_qt_col=0. ; d_ek_col=0. ; d_h_dair_col =0.
+    d_h_qw_col=0. ; d_h_ql_col=0. ; d_h_qs_col=0. ; d_h_qbs_col=0. ; d_h_col=0.
 
     ! Outputs used in cloudth_vert
@@ -223,6 +225,6 @@
     deallocate(sza_o)
     deallocate (bils_ec,bils_ech,bils_tke,bils_diss,bils_kinetic,bils_enthalp,bils_latent)
-    deallocate (d_qw_col, d_ql_col, d_qs_col, d_qt_col, d_ek_col, d_h_dair_col &
-  &           , d_h_qw_col, d_h_ql_col, d_h_qs_col, d_h_col)
+    deallocate (d_qw_col, d_ql_col, d_qs_col, d_qbs_col, d_qt_col, d_ek_col, d_h_dair_col &
+  &           , d_h_qw_col, d_h_ql_col, d_h_qs_col, d_h_qbs_col, d_h_col)
 
     ! Outputs used in cloudth_vert
Index: LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmdiso/physiq_mod.F90
===================================================================
--- LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmdiso/physiq_mod.F90	(revision 4660)
+++ LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmdiso/physiq_mod.F90	(revision 4669)
@@ -79,10 +79,11 @@
     USE wxios, ONLY: g_ctx, wxios_set_context
 #endif
-    USE atke_turbulence_ini_mod, ONLY : atke_ini
-    USE lscp_ini_mod, ONLY : lscp_ini
     USE lscp_mod, ONLY : lscp
     USE wake_ini_mod, ONLY : wake_ini
     USE yamada_ini_mod, ONLY : yamada_ini
+    USE atke_turbulence_ini_mod, ONLY : atke_ini
     USE thermcell_ini_mod, ONLY : thermcell_ini
+    USE blowing_snow_ini_mod, ONLY : blowing_snow_ini , qbst_bs 
+    USE lscp_ini_mod, ONLY : lscp_ini
 
     !USE cmp_seri_mod
@@ -183,8 +184,8 @@
        ! [Variables internes non sauvegardees de la physique]
        ! Variables locales pour effectuer les appels en serie
-       t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,tr_seri,rneb_seri, &
+       t_seri,q_seri,ql_seri,qs_seri,qbs_seri,u_seri,v_seri,tr_seri,rneb_seri, &
        ! Dynamic tendencies (diagnostics)
-       d_t_dyn,d_q_dyn,d_ql_dyn,d_qs_dyn,d_u_dyn,d_v_dyn,d_tr_dyn,d_rneb_dyn, &
-       d_q_dyn2d,d_ql_dyn2d,d_qs_dyn2d, &
+       d_t_dyn,d_q_dyn,d_ql_dyn,d_qs_dyn,d_qbs_dyn,d_u_dyn,d_v_dyn,d_tr_dyn,d_rneb_dyn, &
+       d_q_dyn2d,d_ql_dyn2d,d_qs_dyn2d,d_qbs_dyn2d, &
        ! Physic tendencies
        d_t_con,d_q_con,d_u_con,d_v_con, &
@@ -203,8 +204,10 @@
        plul_st,plul_th, &
        !
-       d_t_vdf,d_q_vdf,d_u_vdf,d_v_vdf,d_t_diss, &
+       d_t_vdf,d_q_vdf, d_qbs_vdf, d_u_vdf,d_v_vdf,d_t_diss, &
        d_t_vdf_x, d_t_vdf_w, &
        d_q_vdf_x, d_q_vdf_w, &
        d_ts, &
+       !
+       d_t_bs,d_q_bs,d_qbs_bs, &
        !
 !       d_t_oli,d_u_oli,d_v_oli, &
@@ -254,6 +257,6 @@
        cldh, cldl,cldm, cldq, cldt,      &
        JrNt,                             &
-       dthmin, evap, fder, plcl, plfc,   &
-       prw, prlw, prsw,                  &
+       dthmin, evap, snowerosion,fder, plcl, plfc,   &
+       prw, prlw, prsw, prbsw,                  &
        s_lcl, s_pblh, s_pblt, s_therm,   &
        cdragm, cdragh,                   &
@@ -328,4 +331,6 @@
        !
        rneblsvol, &
+       pfraclr,pfracld, &
+       distcltop, &
        zqsatl, zqsats, &
        qclr, qcld, qss, qvc, rnebclr, rnebss, gamma_ss, &
@@ -343,5 +348,5 @@
        fsolsw, wfbils, wfbilo,  &
        wfevap, wfrain, wfsnow,  &  
-       prfl, psfl, fraca, Vprecip,  &
+       prfl, psfl,bsfl, fraca, Vprecip,  &
        zw2,  &
        !
@@ -355,5 +360,6 @@
        beta_prec,  &
        rneb,  &
-       zxsnow,snowhgt,qsnow,to_ice,sissnow,runoff,albsol3_lic
+       zxsnow,snowhgt,qsnow,to_ice,sissnow,runoff,albsol3_lic, &
+       zxfluxt,zxfluxq
 
 
@@ -526,7 +532,7 @@
     !======================================================================
     !
-    ! indices de traceurs eau vapeur, liquide, glace, fraction nuageuse LS (optional)
-    INTEGER,SAVE :: ivap, iliq, isol, irneb
-!$OMP THREADPRIVATE(ivap, iliq, isol, irneb)
+    ! indices de traceurs eau vapeur, liquide, glace, fraction nuageuse LS (optional), blowing snow (optional)
+    INTEGER,SAVE :: ivap, iliq, isol, irneb, ibs
+!$OMP THREADPRIVATE(ivap, iliq, isol, irneb, ibs)
     !
     !
@@ -905,11 +911,13 @@
     REAL dialiq(klon,klev)  ! eau liquide nuageuse
     REAL diafra(klon,klev)  ! fraction nuageuse
-    REAL cldliq(klon,klev)  ! eau liquide nuageuse
+    REAL radocond(klon,klev)  ! eau condensee nuageuse
     !
     !XXX PB 
     REAL fluxq(klon,klev, nbsrf)   ! flux turbulent d'humidite
-    !
-    REAL zxfluxt(klon, klev)
-    REAL zxfluxq(klon, klev)
+    REAL fluxqbs(klon,klev, nbsrf)   ! flux turbulent de neige soufflee
+    !
+    !FC    REAL zxfluxt(klon, klev)
+    !FC    REAL zxfluxq(klon, klev)
+    REAL zxfluxqbs(klon,klev)
     REAL zxfluxu(klon, klev)
     REAL zxfluxv(klon, klev)
@@ -1009,5 +1017,5 @@
     !
     ! tendance nulles
-    REAL, dimension(klon,klev):: du0, dv0, dt0, dq0, dql0, dqi0
+    REAL, dimension(klon,klev):: du0, dv0, dt0, dq0, dql0, dqi0, dqbs0
 #ifdef ISO
     REAL, dimension(ntraciso,klon,klev):: dxt0, dxtl0, dxti0
@@ -1156,4 +1164,6 @@
     REAL ztsol(klon)
     REAL q2m(klon,nbsrf)  ! humidite a 2m
+    REAL fsnowerosion(klon,nbsrf) ! blowing snow flux at surface
+    REAL qbsfra  ! blowing snow fraction
 #ifdef ISO
     REAL d_xtw(ntraciso),d_xtl(ntraciso), d_xts(ntraciso)
@@ -1234,5 +1244,4 @@
     !IM 100106 BEG : pouvoir sortir les ctes de la physique
     include "conema3.h"
-    include "fisrtilp.h"
     include "nuage.h"
     include "compbl.h"
@@ -1372,4 +1381,5 @@
        isol = strIdx(tracers(:)%name, addPhase('H2O', 's'))
        irneb= strIdx(tracers(:)%name, addPhase('H2O', 'r'))
+       ibs  = strIdx(tracers(:)%name, addPhase('H2O', 'b'))
        CALL init_etat0_limit_unstruct
        IF (.NOT. create_etat0_limit) CALL init_limit_read(days_elapsed)
@@ -1421,5 +1431,5 @@
        ENDIF
 
-       IF (ok_ice_sursat.AND.(nqo.NE.4)) THEN
+       IF (ok_ice_sursat.AND.(nqo.LT.4)) THEN
           WRITE (lunout, *) ' ok_ice_sursat=y requires 4 H2O tracers ', &
                '(H2O_g, H2O_l, H2O_s, H2O_r) but nqo=', nqo, '. Might as well stop here.'
@@ -1440,4 +1450,14 @@
        ENDIF
 
+        IF (ok_bs) THEN
+          abort_message='blowing snow cannot be activated with water isotopes yet'
+          CALL abort_physic(modname,abort_message, 1)
+         IF ((ok_ice_sursat.AND.nqo .LT.5).OR.(.NOT.ok_ice_sursat.AND.nqo.LT.4)) THEN
+             WRITE (lunout, *) 'activation of blowing snow needs a specific H2O tracer', &
+                               'but nqo=', nqo
+             abort_message='see above'
+             CALL abort_physic(modname,abort_message, 1)
+         ENDIF
+        ENDIF
        Ncvpaseq1 = 0
        dnwd0=0.0
@@ -1890,7 +1910,8 @@
        CALL thermcell_ini(iflag_thermals,prt_level,tau_thermals,lunout, &
    &    RG,RD,RCPD,RKAPPA,RLVTT,RETV)
-       IF (ok_new_lscp) then
-           CALL lscp_ini(pdtphys,ok_ice_sursat)
-       endif
+       CALL lscp_ini(pdtphys,ok_ice_sursat,RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT,RD,RG)
+       CALL blowing_snow_ini(prt_level,lunout, &
+                             RCPD, RLSTT, RLVTT, RLMLT, &
+                             RVTMP2, RTT,RD,RG)
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -1919,5 +1940,5 @@
        CALL phys_output_write(itap, pdtphys, paprs, pphis,                    &
                               pplay, lmax_th, aerosol_couple,                 &
-                              ok_ade, ok_aie, ok_volcan, ivap, iliq, isol, ok_sync,&
+                              ok_ade, ok_aie, ok_volcan, ivap, iliq, isol, ibs,  ok_sync,&
                               ptconv, read_climoz, clevSTD,                   &
                               ptconvth, d_u, d_t, qx, d_qx, zmasse,           &
@@ -2406,4 +2427,5 @@
     dql0(:,:)=0.
     dqi0(:,:)=0.
+    dqbs0(:,:)=0.
 #ifdef ISO
       dxt0(:,:,:)=0.
@@ -2463,4 +2485,5 @@
           q_seri(i,k)  = qx(i,k,ivap)
           ql_seri(i,k) = qx(i,k,iliq)
+          qbs_seri(i,k) = 0.
           !CR: ATTENTION, on rajoute la variable glace
           IF (nqo.EQ.2) THEN             !--vapour and liquid only
@@ -2470,8 +2493,16 @@
              qs_seri(i,k) = qx(i,k,isol)
              rneb_seri(i,k) = 0.
-          ELSE IF (nqo.EQ.4) THEN        !--vapour, liquid, ice and rneb
+          ELSE IF (nqo.GE.4) THEN        !--vapour, liquid, ice and rneb and blowing snow
              qs_seri(i,k) = qx(i,k,isol)
+             IF (ok_ice_sursat) THEN
              rneb_seri(i,k) = qx(i,k,irneb)
+             ENDIF
+             IF (ok_bs) THEN
+             qbs_seri(i,k)= qx(i,k,ibs)
+             ENDIF
+
           ENDIF
+
+
        ENDDO
     ENDDO
@@ -2515,5 +2546,5 @@
     qql1(:)=0.0
     DO k = 1, klev
-      qql1(:)=qql1(:)+(q_seri(:,k)+ql_seri(:,k)+qs_seri(:,k))*zmasse(:,k)
+      qql1(:)=qql1(:)+(q_seri(:,k)+ql_seri(:,k)+qs_seri(:,k)+qbs_seri(:,k))*zmasse(:,k)
     ENDDO
 #ifdef ISO
@@ -2631,4 +2662,5 @@
        d_ql_dyn(:,:) = (ql_seri(:,:)-ql_ancien(:,:))/phys_tstep
        d_qs_dyn(:,:) = (qs_seri(:,:)-qs_ancien(:,:))/phys_tstep
+       d_qbs_dyn(:,:) = (qbs_seri(:,:)-qbs_ancien(:,:))/phys_tstep
        CALL water_int(klon,klev,q_seri,zmasse,zx_tmp_fi2d)
        d_q_dyn2d(:)=(zx_tmp_fi2d(:)-prw_ancien(:))/phys_tstep
@@ -2637,4 +2669,6 @@
        CALL water_int(klon,klev,qs_seri,zmasse,zx_tmp_fi2d)
        d_qs_dyn2d(:)=(zx_tmp_fi2d(:)-prsw_ancien(:))/phys_tstep
+       CALL water_int(klon,klev,qbs_seri,zmasse,zx_tmp_fi2d)
+       d_qbs_dyn2d(:)=(zx_tmp_fi2d(:)-prbsw_ancien(:))/phys_tstep
        ! !! RomP >>>   td dyn traceur
        IF (nqtot > nqo) d_tr_dyn(:,:,:)=(tr_seri(:,:,:)-tr_ancien(:,:,:))/phys_tstep
@@ -2723,4 +2757,5 @@
        d_ql_dyn2d(:) = 0.0
        d_qs_dyn2d(:) = 0.0
+       d_qbs_dyn2d(:)= 0.0
 
 #ifdef ISO
@@ -2746,4 +2781,5 @@
        ! !! RomP <<<
        d_rneb_dyn(:,:)=0.0
+       d_qbs_dyn(:,:)=0.0
        ancien_ok = .TRUE.
     ENDIF
@@ -2878,5 +2914,5 @@
 
      CALL add_phys_tend &
-            (du0,dv0,d_t_eva,d_q_eva,d_ql_eva,d_qi_eva,paprs,&
+            (du0,dv0,d_t_eva,d_q_eva,d_ql_eva,d_qi_eva,dqbs0,paprs,&
                'eva',abortphy,flag_inhib_tend,itap,0 &
 #ifdef ISO
@@ -3047,7 +3083,7 @@
             longitude_deg, latitude_deg, rugoro,  zrmu0,      &
             sollwdown,    cldt,      &
-            rain_fall, snow_fall, solsw,   solswfdiff, sollw,     &
+            rain_fall, snow_fall, bs_fall, solsw,   solswfdiff, sollw,     &
             gustiness,                                &
-            t_seri,    q_seri,    u_seri,  v_seri,    &
+            t_seri,    q_seri,   qbs_seri,  u_seri,  v_seri,    &
                                 !nrlmd+jyg<
             wake_deltat, wake_deltaq, wake_cstar, wake_s, &
@@ -3060,9 +3096,9 @@
                                 !albedo SB >>>
                                 ! albsol1,   albsol2,   sens,    evap,      &
-            albsol_dir,   albsol_dif,   sens,    evap,   &  
+            albsol_dir,   albsol_dif,   sens,    evap, snowerosion, &  
                                 !albedo SB <<<
             albsol3_lic,runoff,   snowhgt,   qsnow, to_ice, sissnow, &
             zxtsol,    zxfluxlat, zt2m,    qsat2m,  zn2mout, &
-            d_t_vdf,   d_q_vdf,   d_u_vdf, d_v_vdf, d_t_diss, &
+            d_t_vdf,   d_q_vdf, d_qbs_vdf,  d_u_vdf, d_v_vdf, d_t_diss, &
                                 !nrlmd<
                                 !jyg<
@@ -3090,5 +3126,5 @@
             fluxt,   fluxu,  fluxv, &
             dsens,     devap,     zxsnow, &
-            zxfluxt,   zxfluxq,   q2m,     fluxq, pbl_tke, &
+            zxfluxt,   zxfluxq,  zxfluxqbs,  q2m, fluxq, fluxqbs, pbl_tke, &
                                 !nrlmd+jyg<
             wake_delta_pbl_TKE, &
@@ -3202,5 +3238,5 @@
        IF (klon_glo==1) THEN
           CALL add_pbl_tend &
-               (d_u_vdf,d_v_vdf,d_t_vdf+d_t_diss,d_q_vdf,dql0,dqi0,paprs,&
+               (d_u_vdf,d_v_vdf,d_t_vdf+d_t_diss,d_q_vdf,dql0,dqi0,d_qbs_vdf,paprs,&
                'vdf',abortphy,flag_inhib_tend,itap &
 #ifdef ISO
@@ -3210,5 +3246,5 @@
        ELSE
           CALL add_phys_tend &
-               (d_u_vdf,d_v_vdf,d_t_vdf+d_t_diss,d_q_vdf,dql0,dqi0,paprs,&
+               (d_u_vdf,d_v_vdf,d_t_vdf+d_t_diss,d_q_vdf,dql0,dqi0,d_qbs_vdf,paprs,&
                'vdf',abortphy,flag_inhib_tend,itap,0 &
 #ifdef ISO
@@ -3272,4 +3308,28 @@
 
     ENDIF
+
+    ! ==================================================================
+    ! Blowing snow sublimation and sedimentation
+
+    d_t_bs(:,:)=0.
+    d_q_bs(:,:)=0.
+    d_qbs_bs(:,:)=0.
+    bsfl(:,:)=0.
+    bs_fall(:)=0.
+    IF (ok_bs) THEN
+
+     CALL call_blowing_snow_sublim_sedim(klon,klev,phys_tstep,t_seri,q_seri,qbs_seri,pplay,paprs, &
+                                        d_t_bs,d_q_bs,d_qbs_bs,bsfl,bs_fall)
+
+     CALL add_phys_tend &
+               (du0,dv0,d_t_bs,d_q_bs,dql0,dqi0,d_qbs_bs,paprs,&
+               'bs',abortphy,flag_inhib_tend,itap,0  &
+#ifdef ISO                                       
+       &    ,dxt0,dxtl0,dxti0 &                    
+#endif                                       
+       &   )
+
+    ENDIF
+
     ! =================================================================== c
     !   Calcul de Qsat
@@ -3959,5 +4019,5 @@
 !!
 !!
-    CALL add_phys_tend(d_u_con, d_v_con, d_t_con, d_q_con, dql0, dqi0, paprs, &
+    CALL add_phys_tend(d_u_con, d_v_con, d_t_con, d_q_con, dql0, dqi0, dqbs0, paprs, &
          'convection',abortphy,flag_inhib_tend,itap,0 &
 #ifdef ISO
@@ -4191,5 +4251,5 @@
        !-----------------------------------------------------------------------
        ! ajout des tendances des poches froides
-       CALL add_phys_tend(du0,dv0,d_t_wake,d_q_wake,dql0,dqi0,paprs,'wake', &
+       CALL add_phys_tend(du0,dv0,d_t_wake,d_q_wake,dql0,dqi0,dqbs0,paprs,'wake', &
             abortphy,flag_inhib_tend,itap,0 &
 #ifdef ISO
@@ -4470,5 +4530,5 @@
           !
           CALL add_phys_tend(d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs,  &
-                             dql0,dqi0,paprs,'thermals', abortphy,flag_inhib_tend,itap,0 &
+                             dql0,dqi0,dqbs0,paprs,'thermals', abortphy,flag_inhib_tend,itap,0 &
 #ifdef ISO
      &    ,d_xt_ajs,dxtl0,dxti0 &
@@ -4586,5 +4646,5 @@
           !--------------------------------------------------------------------
           ! ajout des tendances de l'ajustement sec ou des thermiques
-          CALL add_phys_tend(du0,dv0,d_t_ajsb,d_q_ajsb,dql0,dqi0,paprs, &
+          CALL add_phys_tend(du0,dv0,d_t_ajsb,d_q_ajsb,dql0,dqi0,dqbs0,paprs, &
                'ajsb',abortphy,flag_inhib_tend,itap,0 &
 #ifdef ISO
@@ -4627,5 +4687,5 @@
          ptconv,ptconvth,clwcon0th, rnebcon0th,     &
          paprs,pplay,t_seri,q_seri, qtc_cv, sigt_cv, zqsat, &
-         pbl_tke(:,:,is_ave),tke_dissip_ave,l_mix_ave,wprime_ave,t2m,q2m,fm_therm, &
+         pbl_tke(:,:,is_ave),tke_dissip_ave,l_mix_ave,wprime_ave,t2m,q2m,fm_therm,cell_area, &
          ratqs,ratqsc,ratqs_inter)
 
@@ -4731,9 +4791,10 @@
          t_seri, q_seri,ptconv,ratqs, &
          d_t_lsc, d_q_lsc, d_ql_lsc, d_qi_lsc, rneb, rneblsvol, rneb_seri, &
-         cldliq, picefra, rain_lsc, snow_lsc, &
+         pfraclr,pfracld, &
+         radocond, picefra, rain_lsc, snow_lsc, &
          frac_impa, frac_nucl, beta_prec_fisrt, &
          prfl, psfl, rhcl,  &
          zqasc, fraca,ztv,zpspsk,ztla,zthl,iflag_cld_th, &
-         iflag_ice_thermo, ok_ice_sursat, zqsatl, zqsats, &
+         iflag_ice_thermo, ok_ice_sursat, zqsatl, zqsats, distcltop, &
          qclr, qcld, qss, qvc, rnebclr, rnebss, gamma_ss, &
          Tcontr, qcontr, qcontr2, fcontrN, fcontrP )
@@ -4743,5 +4804,5 @@
     CALL fisrtilp(phys_tstep,paprs,pplay, &
          t_seri, q_seri,ptconv,ratqs, &
-         d_t_lsc, d_q_lsc, d_ql_lsc, d_qi_lsc, rneb, cldliq, &
+         d_t_lsc, d_q_lsc, d_ql_lsc, d_qi_lsc, rneb, radocond, &
          rain_lsc, snow_lsc, &
          pfrac_impa, pfrac_nucl, pfrac_1nucl, &
@@ -4791,5 +4852,5 @@
 !    write(*,9000) "rcpv","rcw",rcpv,rcw,rcs,t_seri(1,1)
 !-JLD 
-    CALL add_phys_tend(du0,dv0,d_t_lsc,d_q_lsc,d_ql_lsc,d_qi_lsc,paprs, &
+    CALL add_phys_tend(du0,dv0,d_t_lsc,d_q_lsc,d_ql_lsc,d_qi_lsc,dqbs0,paprs, &
          'lsc',abortphy,flag_inhib_tend,itap,0 &
 #ifdef ISO
@@ -4828,15 +4889,39 @@
     ENDIF
 
-    !---------------------------------------------------------------------------
+
+!---------------------------------------------------------------------------
     DO k = 1, klev
        DO i = 1, klon
           cldfra(i,k) = rneb(i,k)
           !CR: a quoi ca sert? Faut-il ajouter qs_seri?
-          IF (.NOT.new_oliq) cldliq(i,k) = ql_seri(i,k)
+          !EV: en effet etrange, j'ajouterais aussi qs_seri
+          !    plus largement, je nettoierais (enleverrais) ces lignes
+          IF (.NOT.new_oliq) radocond(i,k) = ql_seri(i,k)
        ENDDO
     ENDDO
 
 
-
+    ! Option to activate the radiative effect of blowing snow (ok_rad_bs)
+    ! makes sense only if the new large scale condensation scheme is active
+    ! with the ok_icefra_lscp flag active as well
+
+    IF (ok_bs .AND. ok_rad_bs) THEN
+       IF (ok_new_lscp .AND. ok_icefra_lscp) THEN
+           DO k=1,klev
+             DO i=1,klon
+                radocond(i,k)=radocond(i,k)+qbs_seri(i,k)
+                picefra(i,k)=(radocond(i,k)*picefra(i,k)+qbs_seri(i,k))/(radocond(i,k))
+                qbsfra=min(qbs_seri(i,k)/qbst_bs,1.0) 
+                cldfra(i,k)=max(cldfra(i,k),qbsfra)
+             ENDDO
+           ENDDO
+       ELSE
+          WRITE(lunout,*)"PAY ATTENTION, you try to activate the radiative effect of blowing snow"
+          WRITE(lunout,*)"with ok_new_lscp=false and/or ok_icefra_lscp=false"
+          abort_message='inconsistency in cloud phase for blowing snow'
+          CALL abort_physic(modname,abort_message,1)
+       ENDIF
+
+    ENDIF
 #ifdef ISO      
 !#ifdef ISOVERIF
@@ -5001,5 +5086,5 @@
           DO i = 1, klon
              IF (diafra(i,k).GT.cldfra(i,k)) THEN
-                cldliq(i,k) = dialiq(i,k)
+                radocond(i,k) = dialiq(i,k)
                 cldfra(i,k) = diafra(i,k)
              ENDIF
@@ -5038,9 +5123,9 @@
                 DO i=1,klon
                    IF (ptconv(i,k).AND.ptconvth(i,k)) THEN
-                      cldliq(i,k)=cldliq(i,k)+rnebcon(i,k)*clwcon(i,k)
+                      radocond(i,k)=radocond(i,k)+rnebcon(i,k)*clwcon(i,k)
                       cldfra(i,k)=min(cldfra(i,k)+rnebcon(i,k),1.)
                    ELSE IF (ptconv(i,k)) THEN
                       cldfra(i,k)=rnebcon(i,k)
-                      cldliq(i,k)=rnebcon(i,k)*clwcon(i,k)
+                      radocond(i,k)=rnebcon(i,k)*clwcon(i,k)
                    ENDIF
                 ENDDO
@@ -5051,5 +5136,5 @@
                 DO i=1,klon
                    cldfra(i,k)=min(cldfra(i,k)+rnebcon(i,k),1.)
-                   cldliq(i,k)=cldliq(i,k)+rnebcon(i,k)*clwcon(i,k)
+                   radocond(i,k)=radocond(i,k)+rnebcon(i,k)*clwcon(i,k)
                 ENDDO
              ENDDO
@@ -5069,5 +5154,5 @@
                    IF (ptconv(i,k).AND. .NOT.ptconvth(i,k)) THEN
                       cldfra(i,k)=rnebcon(i,k)
-                      cldliq(i,k)=rnebcon(i,k)*clwcon(i,k)
+                      radocond(i,k)=rnebcon(i,k)*clwcon(i,k)
                    ENDIF
                 ENDDO
@@ -5080,5 +5165,5 @@
           ! Ancienne version
           cldfra(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
-          cldliq(:,:)=cldliq(:,:)+rnebcon(:,:)*clwcon(:,:)
+          radocond(:,:)=radocond(:,:)+rnebcon(:,:)*clwcon(:,:)
        ENDIF
 
@@ -5100,5 +5185,5 @@
           DO i = 1, klon
              IF (diafra(i,k).GT.cldfra(i,k)) THEN
-                cldliq(i,k) = dialiq(i,k)
+                radocond(i,k) = dialiq(i,k)
                 cldfra(i,k) = diafra(i,k)
              ENDIF
@@ -5475,17 +5560,17 @@
           ENDIF
           CALL newmicro (flag_aerosol, ok_cdnc, bl95_b0, bl95_b1, &
-               paprs, pplay, t_seri, cldliq, picefra, cldfra, &
+               paprs, pplay, t_seri, radocond, picefra, cldfra, &
                cldtau, cldemi, cldh, cldl, cldm, cldt, cldq, &
                flwp, fiwp, flwc, fiwc, &
                mass_solu_aero, mass_solu_aero_pi, &
-               cldtaupi, latitude_deg, re, fl, ref_liq, ref_ice, &
+               cldtaupi, latitude_deg, distcltop, re, fl, ref_liq, ref_ice, &
                ref_liq_pi, ref_ice_pi)
        ELSE
           CALL nuage (paprs, pplay, &
-               t_seri, cldliq, picefra, cldfra, cldtau, cldemi, &
+               t_seri, radocond, picefra, cldfra, cldtau, cldemi, &
                cldh, cldl, cldm, cldt, cldq, &
                ok_aie, &
                mass_solu_aero, mass_solu_aero_pi, &
-               bl95_b0, bl95_b1, &
+               bl95_b0, bl95_b1, distcltop, &
                cldtaupi, re, fl)
        ENDIF
@@ -5792,5 +5877,5 @@
     ENDDO
 
-    CALL add_phys_tend(du0,dv0,d_t_swr,dq0,dql0,dqi0,paprs,'SW',abortphy,flag_inhib_tend,itap,0 &
+    CALL add_phys_tend(du0,dv0,d_t_swr,dq0,dql0,dqi0,dqbs0,paprs,'SW',abortphy,flag_inhib_tend,itap,0 &
 #ifdef ISO
      &    ,dxt0,dxtl0,dxti0 &
@@ -5798,5 +5883,5 @@
      &  )
     CALL prt_enerbil('SW',itap)
-    CALL add_phys_tend(du0,dv0,d_t_lwr,dq0,dql0,dqi0,paprs,'LW',abortphy,flag_inhib_tend,itap,0 &
+    CALL add_phys_tend(du0,dv0,d_t_lwr,dq0,dql0,dqi0,dqbs0,paprs,'LW',abortphy,flag_inhib_tend,itap,0 &
 #ifdef ISO
      &    ,dxt0,dxtl0,dxti0 &
@@ -5846,5 +5931,5 @@
           ! -> condition on zrel_oro can deactivate the drag on tilted planar terrains
           !    such as ice sheets (work by V. Wiener)
-          ! zpmm_orodr_t and zstd_orodr_t are activation thresholds set by F. Lott to 
+          ! zpmm_orodr_t and zstd_orodr_t are activation thresholds set by F. Lott to
           ! earn computation time but they are not physical.
           IF (((zpic(i)-zmea(i)).GT.zpmm_orodr_t).AND.(zstd(i).GT.zstd_orodr_t).AND.(zrel_oro(i).LE.zrel_oro_t)) THEN
@@ -5877,5 +5962,5 @@
        !-----------------------------------------------------------------------
        ! ajout des tendances de la trainee de l'orographie
-       CALL add_phys_tend(d_u_oro,d_v_oro,d_t_oro,dq0,dql0,dqi0,paprs,'oro', &
+       CALL add_phys_tend(d_u_oro,d_v_oro,d_t_oro,dq0,dql0,dqi0,dqbs0,paprs,'oro', &
             abortphy,flag_inhib_tend,itap,0 &
 #ifdef ISO
@@ -5932,5 +6017,5 @@
 
        ! ajout des tendances de la portance de l'orographie
-       CALL add_phys_tend(d_u_lif, d_v_lif, d_t_lif, dq0, dql0, dqi0, paprs, &
+       CALL add_phys_tend(d_u_lif, d_v_lif, d_t_lif, dq0, dql0, dqi0, dqbs0,paprs, &
             'lif', abortphy,flag_inhib_tend,itap,0 &
 #ifdef ISO
@@ -5961,5 +6046,5 @@
        d_t_hin(:, :)=0.
        CALL add_phys_tend(du_gwd_hines, dv_gwd_hines, d_t_hin, dq0, dql0, &
-            dqi0, paprs, 'hin', abortphy,flag_inhib_tend,itap,0 &
+            dqi0, dqbs0,paprs, 'hin', abortphy,flag_inhib_tend,itap,0 &
 #ifdef ISO
      &    ,dxt0,dxtl0,dxti0 &
@@ -5983,5 +6068,5 @@
        ENDDO
 
-       CALL add_phys_tend(du_gwd_front, dv_gwd_front, dt0, dq0, dql0, dqi0, &
+       CALL add_phys_tend(du_gwd_front, dv_gwd_front, dt0, dq0, dql0, dqi0, dqbs0, &
             paprs, 'front_gwd_rando', abortphy,flag_inhib_tend,itap,0 &
 #ifdef ISO
@@ -5996,5 +6081,5 @@
             rain_fall + snow_fall, zustr_gwd_rando, zvstr_gwd_rando, &
             du_gwd_rando, dv_gwd_rando, east_gwstress, west_gwstress)
-       CALL add_phys_tend(du_gwd_rando, dv_gwd_rando, dt0, dq0, dql0, dqi0, &
+       CALL add_phys_tend(du_gwd_rando, dv_gwd_rando, dt0, dq0, dql0, dqi0, dqbs0, &
             paprs, 'flott_gwd_rando', abortphy,flag_inhib_tend,itap,0 &
 #ifdef ISO
@@ -6055,5 +6140,5 @@
        d_xt_ch4_dtime(:,:,:) = d_xt_ch4(:,:,:)*phys_tstep
 #endif
-       CALL add_phys_tend(du0, dv0, dt0, d_q_ch4_dtime, dql0, dqi0, paprs, &
+       CALL add_phys_tend(du0, dv0, dt0, d_q_ch4_dtime, dql0, dqi0, dqbs0, paprs, &
             'q_ch4', abortphy,flag_inhib_tend,itap,0 &
 #ifdef ISO
@@ -6121,5 +6206,5 @@
 ! car on peut s'attendre a ce que les petites echelles produisent aussi de la TKE
 ! Mais attention, cela ne va pas dans le sens de la conservation de l'energie! 
-          IF ((zstd(i).GT.1.0).AND.(zrel_oro(i).LE.zrel_oro_t)) THEN
+          IF ((zstd(i).GT.1.0) .AND.(zrel_oro(i).LE.zrel_oro_t)) THEN
              itest(i)=1
              igwd=igwd+1
@@ -6381,5 +6466,5 @@
     !
 
-    IF (type_trac=='repr') THEN
+    IF (type_trac == 'repr') THEN
 !MM pas d'impact, car on recupere q_seri,tr_seri,t_seri via phys_local_var_mod
 !MM                               dans Reprobus
@@ -6431,5 +6516,5 @@
          presnivs, pphis,     pphi,     albsol1, &
          sh_in,   ch_in,    rhcl,      cldfra,   rneb, &
-         diafra,   cldliq,    itop_con, ibas_con, &
+         diafra,   radocond,    itop_con, ibas_con, &
          pmflxr,   pmflxs,    prfl,     psfl, &
          da,       phi,       mp,       upwd, &
@@ -6448,5 +6533,5 @@
 
      CALL add_phys_tend &
-            (du0,dv0,dt0,d_q_rep,d_ql_rep,d_qi_rep,paprs,&
+            (du0,dv0,dt0,d_q_rep,d_ql_rep,d_qi_rep,dqbs0,paprs,&
              'rep',abortphy,flag_inhib_tend,itap,0)
         IF (abortphy==1) Print*,'ERROR ABORT REP'
@@ -6525,11 +6610,14 @@
     !   prlw = colonne eau liquide 
     !   prlw = colonne eau solide
+    !   prbsw = colonne neige soufflee
     prw(:) = 0.
     prlw(:) = 0.
     prsw(:) = 0.
+    prbsw(:) = 0.
     DO k = 1, klev
        prw(:)  = prw(:)  + q_seri(:,k)*zmasse(:,k)
        prlw(:) = prlw(:) + ql_seri(:,k)*zmasse(:,k)
        prsw(:) = prsw(:) + qs_seri(:,k)*zmasse(:,k)
+       prbsw(:)= prbsw(:) + qbs_seri(:,k)*zmasse(:,k)
     ENDDO
 
@@ -6602,7 +6690,13 @@
           ENDIF
           !--ice_sursat: nqo=4, on ajoute rneb
-          IF (nqo == 4) THEN
+          IF (nqo.ge.4 .and. ok_ice_sursat) THEN
              d_qx(i,k,irneb) = ( rneb_seri(i,k) - qx(i,k,irneb) ) / phys_tstep
           ENDIF
+
+           IF (nqo.ge.4 .and. ok_bs) THEN
+             d_qx(i,k,ibs) = ( qbs_seri(i,k) - qx(i,k,ibs) ) / phys_tstep
+          ENDIF
+
+
        ENDDO
     ENDDO
@@ -6684,4 +6778,5 @@
     ql_ancien(:,:) = ql_seri(:,:)
     qs_ancien(:,:) = qs_seri(:,:)
+    qbs_ancien(:,:) = qbs_seri(:,:)
     rneb_ancien(:,:) = rneb_seri(:,:)
 #ifdef ISO
@@ -6693,4 +6788,5 @@
     CALL water_int(klon,klev,ql_ancien,zmasse,prlw_ancien)
     CALL water_int(klon,klev,qs_ancien,zmasse,prsw_ancien)
+    CALL water_int(klon,klev,qbs_ancien,zmasse,prbsw_ancien)
     ! !! RomP >>>
     IF (nqtot > nqo) tr_ancien(:,:,:) = tr_seri(:,:,:)
@@ -6821,5 +6917,5 @@
     CALL phys_output_write(itap, pdtphys, paprs, pphis,  &
          pplay, lmax_th, aerosol_couple,                 &
-         ok_ade, ok_aie, ok_volcan, ivap, iliq, isol,    & 
+         ok_ade, ok_aie, ok_volcan, ivap, iliq, isol, ibs,   & 
          ok_sync, ptconv, read_climoz, clevSTD,          &
          ptconvth, d_u, d_t, qx, d_qx, zmasse,           &
Index: LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmdiso/surf_land_mod.F90
===================================================================
--- LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmdiso/surf_land_mod.F90	(revision 4660)
+++ LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmdiso/surf_land_mod.F90	(revision 4669)
@@ -11,5 +11,5 @@
        rlon, rlat, yrmu0, &
        debut, lafin, zlev, ccanopy, swnet, lwnet, albedo, &
-       tsurf, p1lay, cdragh, cdragm, precip_rain, precip_snow, temp_air, spechum, &
+       tsurf, p1lay, cdragh, cdragm, precip_rain, precip_snow, precip_bs, temp_air, spechum, &
        AcoefH, AcoefQ, BcoefH, BcoefQ, &
        AcoefU, AcoefV, BcoefU, BcoefV, & 
@@ -17,5 +17,5 @@
        lwdown_m, q2m, t2m, &
        snow, qsol, agesno, tsoil, &
-       z0m, z0h, SFRWL, alb_dir_new, alb_dif_new, evap, fluxsens, fluxlat, &   
+       z0m, z0h, SFRWL, alb_dir_new, alb_dif_new, evap, fluxsens, fluxlat, fluxbs, &   
        qsurf, tsurf_new, dflux_s, dflux_l, &
        flux_u1, flux_v1 , & 
@@ -94,5 +94,5 @@
     REAL, DIMENSION(klon), INTENT(IN)       :: p1lay
     REAL, DIMENSION(klon), INTENT(IN)       :: cdragh, cdragm
-    REAL, DIMENSION(klon), INTENT(IN)       :: precip_rain, precip_snow
+    REAL, DIMENSION(klon), INTENT(IN)       :: precip_rain, precip_snow, precip_bs
     REAL, DIMENSION(klon), INTENT(IN)       :: temp_air, spechum
     REAL, DIMENSION(klon), INTENT(IN)       :: AcoefH, AcoefQ, BcoefH, BcoefQ
@@ -129,5 +129,5 @@
 !albedo SB <<<
     REAL, DIMENSION(klon), INTENT(OUT)       :: evap
-    REAL, DIMENSION(klon), INTENT(OUT)       :: fluxsens, fluxlat
+    REAL, DIMENSION(klon), INTENT(OUT)       :: fluxsens, fluxlat, fluxbs
     REAL, DIMENSION(klon), INTENT(OUT)       :: qsurf
     REAL, DIMENSION(klon), INTENT(OUT)       :: tsurf_new
@@ -152,4 +152,5 @@
     REAL, DIMENSION(klon) :: tsol_rad, emis_new ! output from interfsol not used
     REAL, DIMENSION(klon) :: u0, v0     ! surface speed
+    REAL, DIMENSION(klon) :: precip_totsnow     ! total solid precip
     INTEGER               :: i
 
@@ -166,4 +167,13 @@
 #endif  
 
+!****************************************************************************************
+!Total solid precip
+
+IF (ok_bs) THEN
+precip_totsnow(:)=precip_snow(:)+precip_bs(:)
+ELSE
+precip_totsnow(:)=precip_snow(:)
+ENDIF
+!****************************************************************************************
 #ifdef ISO
 #ifdef ISOVERIF
@@ -228,7 +238,7 @@
             zlev,  u1, v1, gustiness, temp_air, spechum, epot_air, ccanopy, & 
             cdragh, AcoefH, AcoefQ, BcoefH, BcoefQ, &
-            precip_rain, precip_snow, lwdown_m, swnet, swdown, &
+            precip_rain, precip_totsnow, lwdown_m, swnet, swdown, &
             pref_tmp, q2m, t2m, &
-            evap, fluxsens, fluxlat, &              
+            evap, fluxsens, fluxlat,  &              
             tsol_rad, tsurf_new, alb1_new, alb2_new, &
             emis_new, z0m, z0h, qsurf, &
@@ -277,5 +287,5 @@
 #endif
        CALL surf_land_bucket(itime, jour, knon, knindex, debut, dtime,&
-            tsurf, p1lay, cdragh, precip_rain, precip_snow, temp_air, &
+            tsurf, p1lay, cdragh, precip_rain, precip_totsnow, temp_air, &
             spechum, AcoefH, AcoefQ, BcoefH, BcoefQ, pref, &
             u1, v1, gustiness, rugoro, swnet, lwnet, &
@@ -293,4 +303,6 @@
     ENDIF ! ok_veget
 
+        ! blowing snow not treated yet over land
+        fluxbs(:)=0.
 !****************************************************************************************
 ! Calculation for all land models
Index: LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmdiso/surf_landice_mod.F90
===================================================================
--- LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmdiso/surf_landice_mod.F90	(revision 4660)
+++ LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmdiso/surf_landice_mod.F90	(revision 4669)
@@ -12,10 +12,11 @@
        rmu0, lwdownm, albedo, pphi1, &
        swnet, lwnet, tsurf, p1lay, &
-       cdragh, cdragm, precip_rain, precip_snow, temp_air, spechum, &
+       cdragh, cdragm, precip_rain, precip_snow, precip_bs, temp_air, spechum, &
        AcoefH, AcoefQ, BcoefH, BcoefQ, &
        AcoefU, AcoefV, BcoefU, BcoefV, &
+       AcoefQBS, BcoefQBS, &
        ps, u1, v1, gustiness, rugoro, pctsrf, &
-       snow, qsurf, qsol, agesno, &
-       tsoil, z0m, z0h, SFRWL, alb_dir, alb_dif, evap, fluxsens, fluxlat, &
+       snow, qsurf, qsol, qbs1, agesno, &
+       tsoil, z0m, z0h, SFRWL, alb_dir, alb_dif, evap, fluxsens, fluxlat, fluxbs, &
        tsurf_new, dflux_s, dflux_l, &
        alt, slope, cloudf, &
@@ -34,5 +35,6 @@
     USE cpl_mod,          ONLY : cpl_send_landice_fields
     USE calcul_fluxs_mod
-    USE phys_output_var_mod
+    USE phys_local_var_mod, ONLY : zxrhoslic, zxustartlic
+    USE phys_output_var_mod, ONLY : snow_o,zfra_o
 #ifdef ISO   
     USE fonte_neige_mod,  ONLY : xtrun_off_lic
@@ -48,5 +50,5 @@
 !FC
     USE ioipsl_getin_p_mod, ONLY : getin_p
-
+    USE blowing_snow_ini_mod, ONLY : zeta_bs, pbst_bs, prt_bs, iflag_saltation_bs
 
 #ifdef CPP_INLANDSIS
@@ -72,11 +74,12 @@
     REAL, DIMENSION(klon), INTENT(IN)             :: p1lay
     REAL, DIMENSION(klon), INTENT(IN)             :: cdragh, cdragm
-    REAL, DIMENSION(klon), INTENT(IN)             :: precip_rain, precip_snow
+    REAL, DIMENSION(klon), INTENT(IN)             :: precip_rain, precip_snow, precip_bs
     REAL, DIMENSION(klon), INTENT(IN)             :: temp_air, spechum
     REAL, DIMENSION(klon), INTENT(IN)             :: AcoefH, AcoefQ
     REAL, DIMENSION(klon), INTENT(IN)             :: BcoefH, BcoefQ
     REAL, DIMENSION(klon), INTENT(IN)             :: AcoefU, AcoefV, BcoefU, BcoefV
+    REAL, DIMENSION(klon), INTENT(IN)             :: AcoefQBS, BcoefQBS
     REAL, DIMENSION(klon), INTENT(IN)             :: ps
-    REAL, DIMENSION(klon), INTENT(IN)             :: u1, v1, gustiness
+    REAL, DIMENSION(klon), INTENT(IN)             :: u1, v1, gustiness, qbs1
     REAL, DIMENSION(klon), INTENT(IN)             :: rugoro
     REAL, DIMENSION(klon,nbsrf), INTENT(IN)       :: pctsrf
@@ -118,4 +121,5 @@
 !albedo SB <<<
     REAL, DIMENSION(klon), INTENT(OUT)            :: evap, fluxsens, fluxlat
+    REAL, DIMENSION(klon), INTENT(OUT)            :: fluxbs
     REAL, DIMENSION(klon), INTENT(OUT)            :: tsurf_new
     REAL, DIMENSION(klon), INTENT(OUT)            :: dflux_s, dflux_l      
@@ -179,5 +183,10 @@
 
     REAL,DIMENSION(klon) :: alb1,alb2
+    REAL,DIMENSION(klon) :: precip_totsnow, evap_totsnow
     REAL, DIMENSION (klon,6) :: alb6
+    REAL                   :: rho0, rhoice, ustart0, hsalt, 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
 ! End definition
 !****************************************************************************************
@@ -214,15 +223,13 @@
            PRINT*, 'alb_nir_sno_lic',alb_nir_sno_lic
   
-!  z0m=1.e-3
-!  z0h = z0m
   firstcall=.false.
   ENDIF
 !******************************************************************************************
-!
+
 ! Initialize output variables
     alb3(:) = 999999.
     alb2(:) = 999999.
     alb1(:) = 999999.
-    
+    fluxbs(:)=0.  
     runoff(:) = 0.
 !****************************************************************************************
@@ -232,4 +239,6 @@
     radsol(:) = 0.0
     radsol(1:knon) = swnet(1:knon) + lwnet(1:knon)
+
+!****************************************************************************************
 
 !****************************************************************************************
@@ -327,5 +336,4 @@
 
 
-
     ELSE 
 
@@ -342,6 +350,5 @@
     IF (soil_model) THEN 
        CALL soil(dtime, is_lic, knon, snow, tsurf, qsol, &
-    & longitude(knindex(1:knon)), latitude(knindex(1:knon)), tsoil, soilcap, soilflux)
-
+        & longitude(knindex(1:knon)), latitude(knindex(1:knon)), tsoil, soilcap, soilflux)
        cal(1:knon) = RCPD / soilcap(1:knon)
        radsol(1:knon)  = radsol(1:knon) + soilflux(1:knon)
@@ -406,11 +413,135 @@
          flux_u1, flux_v1)
 
-!****************************************************************************************
-! Calculate snow height, age, run-off,..
+
+!****************************************************************************************
+! Calculate albedo
+!
+!****************************************************************************************
+ 
+    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"
+!       alb1(1 : knon)  = 0.6 !IM cf FH/GK 
+!       alb1(1 : knon)  = 0.82
+!       alb1(1 : knon)  = 0.77 !211003 Ksta0.77
+!       alb1(1 : knon)  = 0.8 !KstaTER0.8 & LMD_ARMIP5
+!IM: KstaTER0.77 & LMD_ARMIP6    
+
+! Attantion: alb1 and alb2 are not the same!
+    alb1(1:knon)  = alb_vis_sno_lic
+    alb2(1:knon)  = alb_nir_sno_lic
+
+
+!****************************************************************************************
+! Rugosity
+!
+!****************************************************************************************
+    z0m = z0m_landice
+    z0h = z0h_landice
+    !z0m = SQRT(z0m**2+rugoro**2)
+
+
+
+  ! Simple blowing snow param
+  if (ok_bs) then
+       ustart0 = 0.211
+       rhoice = 920.0
+       rho0 = 200.0
+       rhomax=450.0
+       rhohard=400.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
+
+       ! computation of threshold friction velocity
+       ! which depends on surface snow density
+       do i = 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))
+           ! 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)) 
+           ! we have multiplied by exp to prevent erosion when rhos>rhomax (usefull till
+           ! rhohard<450)
+       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 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)) * &
+                        (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.)
+          enddo
+
+
+       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)
+          enddo
+       endif
+
+        ! calculation of erosion (emission flux towards the first atmospheric level)
+        ! 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)
+       enddo
+
+  endif
+
+
+
+!****************************************************************************************
+! Calculate surface snow amount
 !    
 !****************************************************************************************
+    IF (ok_bs) THEN
+      precip_totsnow(:)=precip_snow(:)+precip_bs(:)
+      evap_totsnow(:)=evap(:)-fluxbs(:) ! flux bs is positive towards the surface (snow erosion)
+    ELSE
+      precip_totsnow(:)=precip_snow(:)
+      evap_totsnow(:)=evap(:)
+    ENDIF
+
     CALL fonte_neige(knon, is_lic, knindex, dtime, &
-         tsurf, precip_rain, precip_snow, &
-         snow, qsol, tsurf_new, evap &
+         tsurf, precip_rain, precip_totsnow, &
+         snow, qsol, tsurf_new, evap_totsnow &
 #ifdef ISO    
      & ,fq_fonte_diag,fqfonte_diag,snow_evap_diag,fqcalving_diag   &
@@ -446,35 +577,6 @@
 
 
-!****************************************************************************************
-! Calculate albedo
-!
-!****************************************************************************************
-    CALL albsno(klon,knon,dtime,agesno(:),alb_neig(:), precip_snow(:))  
-
-    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)))
-    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"
-!       alb1(1 : knon)  = 0.6 !IM cf FH/GK 
-!       alb1(1 : knon)  = 0.82
-!       alb1(1 : knon)  = 0.77 !211003 Ksta0.77
-!       alb1(1 : knon)  = 0.8 !KstaTER0.8 & LMD_ARMIP5
-!IM: KstaTER0.77 & LMD_ARMIP6    
-
-! Attantion: alb1 and alb2 are not the same!
-    alb1(1:knon)  = alb_vis_sno_lic
-    alb2(1:knon)  = alb_nir_sno_lic
-
-
-!****************************************************************************************
-! Rugosity
-!
-!****************************************************************************************
-    z0m=1.e-3
-    z0h = z0m
-    z0m = SQRT(z0m**2+rugoro**2)
-
+    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)))  
 
 
@@ -494,5 +596,5 @@
           run_off_lic_frac(j) = pctsrf(i,is_lic)
        ENDDO
-       
+
        CALL cpl_send_landice_fields(itime, knon, knindex, run_off_lic, run_off_lic_frac)
     ENDIF
@@ -501,16 +603,4 @@
     runoff(1:knon)=run_off_lic(1:knon)/dtime
 
-  
-!****************************************************************************************
-! Etienne: comment these lines because of duplication just below
-!       snow_o=0.
-!       zfra_o = 0.
-!       DO j = 1, knon
-!           i = knindex(j)
-!           snow_o(i) = snow(j) 
-!           zfra_o(i) = zfra(j)
-!       ENDDO
-!
-!****************************************************************************************
        snow_o=0.
        zfra_o = 0.
@@ -553,6 +643,4 @@
 !albedo SB <<<
 
-  
-  
 
   END SUBROUTINE surf_landice
Index: LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmdiso/surf_ocean_mod.F90
===================================================================
--- LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmdiso/surf_ocean_mod.F90	(revision 4660)
+++ LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmdiso/surf_ocean_mod.F90	(revision 4669)
@@ -13,5 +13,5 @@
        windsp, rmu0, fder, tsurf_in, &
        itime, dtime, jour, knon, knindex, &
-       p1lay, z1lay, cdragh, cdragm, precip_rain, precip_snow, temp_air, spechum, &
+       p1lay, z1lay, cdragh, cdragm, precip_rain, precip_snow, precip_bs, temp_air, spechum, &
        AcoefH, AcoefQ, BcoefH, BcoefQ, &
        AcoefU, AcoefV, BcoefU, BcoefV, &
@@ -72,5 +72,5 @@
     REAL, DIMENSION(klon), INTENT(IN)        :: cdragh
     REAL, DIMENSION(klon), INTENT(IN)        :: cdragm
-    REAL, DIMENSION(klon), INTENT(IN)        :: precip_rain, precip_snow
+    REAL, DIMENSION(klon), INTENT(IN)        :: precip_rain, precip_snow, precip_bs
     REAL, DIMENSION(klon), INTENT(IN)        :: temp_air, spechum
     REAL, DIMENSION(klon), INTENT(IN)        :: AcoefH, AcoefQ, BcoefH, BcoefQ
@@ -169,4 +169,5 @@
     REAL, DIMENSION(klon) :: radsol
     REAL, DIMENSION(klon) :: cdragq ! Cdrag pour l'evaporation
+    REAL, DIMENSION(klon) :: precip_totsnow
     CHARACTER(len=20),PARAMETER :: modname="surf_ocean"
     real rhoa(knon) ! density of moist air  (kg / m3)
@@ -200,4 +201,14 @@
     radsol(1:knon) = swnet(1:knon) + lwnet(1:knon)
 
+
+    !****************************************************************************************
+    !Total solid precip
+
+    IF (ok_bs) THEN
+       precip_totsnow(:)=precip_snow(:)+precip_bs(:)
+    ELSE
+       precip_totsnow(:)=precip_snow(:)
+    ENDIF
+
     !******************************************************************************
     ! Cdragq computed from cdrag
@@ -227,5 +238,5 @@
             windsp, fder, & 
             itime, dtime, knon, knindex, &
-            p1lay, cdragh, cdragq, cdragm, precip_rain, precip_snow,temp_air,spechum,& 
+            p1lay, cdragh, cdragq, cdragm, precip_rain, precip_totsnow,temp_air,spechum,& 
             AcoefH, AcoefQ, BcoefH, BcoefQ, &
             AcoefU, AcoefV, BcoefU, BcoefV, &
@@ -239,5 +250,5 @@
        CALL ocean_slab_noice( &
             itime, dtime, jour, knon, knindex, &
-            p1lay, cdragh, cdragq, cdragm, precip_rain, precip_snow, temp_air, spechum,&
+            p1lay, cdragh, cdragq, cdragm, precip_rain, precip_totsnow, temp_air, spechum,&
             AcoefH, AcoefQ, BcoefH, BcoefQ, &
             AcoefU, AcoefV, BcoefU, BcoefV, &
@@ -250,5 +261,5 @@
        CALL ocean_forced_noice( &
             itime, dtime, jour, knon, knindex, &
-            p1lay, cdragh, cdragq, cdragm, precip_rain, precip_snow, &
+            p1lay, cdragh, cdragq, cdragm, precip_rain, precip_totsnow, &
             temp_air, spechum, &
             AcoefH, AcoefQ, BcoefH, BcoefQ, &
@@ -370,5 +381,5 @@
        call bulk_flux(tkt, tks, taur, dter, dser, t_int, s_int, ds_ns, dt_ns, &
             u = windsp(:knon), t_ocean_1 = tsurf_new(:knon), s1 = sss(:knon), &
-            rain = precip_rain(:knon) + precip_snow(:knon), &
+            rain = precip_rain(:knon) + precip_totsnow(:knon), &
             hf = - fluxsens(:knon), hlb = - fluxlat(:knon), &
             rnl = - lwnet(:knon), &
