Index: /LMDZ6/branches/blowing_snow/DefLists/field_def_lmdz.xml
===================================================================
--- /LMDZ6/branches/blowing_snow/DefLists/field_def_lmdz.xml	(revision 4484)
+++ /LMDZ6/branches/blowing_snow/DefLists/field_def_lmdz.xml	(revision 4485)
@@ -112,6 +112,10 @@
         <field id="pluc"    long_name="Convective Precip."      unit="kg/(s*m2)" />
         <field id="snow"    long_name="Snow fall"               unit="kg/(s*m2)" />
-        <field id="evap"    long_name="Evaporat"                unit="kg/(s*m2)" />
-        <field id="evap_ter"    long_name="evaporation at surface ter"    unit="kg/(s*m2)" />
+	<field id="evap"    long_name="Evaporat"                unit="kg/(s*m2)" />
+        <field id="snowerosion"    long_name="blowing snow erosion"                unit="kg/(s*m2)" />
+	<field id="bsfall"    long_name="blowing snow precipitation"                unit="kg/(s*m2)" />
+	<field id="rhosnow_lic"    long_name="snow density lic"                unit="kg/m3" />
+        <field id="ustart_lic"    long_name="ustar threshold for snow erosion"                unit="m/s" />
+	<field id="evap_ter"    long_name="evaporation at surface ter"    unit="kg/(s*m2)" />
         <field id="evap_lic"    long_name="evaporation at surface lic"    unit="kg/(s*m2)" />
         <field id="evap_oce"    long_name="evaporation at surface oce"    unit="kg/(s*m2)" />
@@ -275,5 +279,6 @@
         <field id="prw"    long_name="Precipitable water"    unit="kg/m2" />
         <field id="prlw"   long_name="Precipitable liquid water"  unit="kg/m2" />
-        <field id="prsw"   long_name="Precipitable solid water"  unit="kg/m2" />
+	<field id="prsw"   long_name="Precipitable solid water"  unit="kg/m2" />
+	<field id="prbsw"   long_name="Precipitable blowing snow"  unit="kg/m2" />
         <field id="s_pblh"    long_name="Boundary Layer Height"    unit="m" />
         <field id="s_pblt"    long_name="t at Boundary Layer Height"    unit="K" />
@@ -505,9 +510,11 @@
         <field id="dqphy2d"    long_name="Physics dQ"    unit="(kg/m2)/s" />
         <field id="dqlphy2d"    long_name="Physics dQL"    unit="(kg/m2)/s" />
-        <field id="dqsphy2d"    long_name="Physics dQS"    unit="(kg/m2)/s" />
+	<field id="dqsphy2d"    long_name="Physics dQS"    unit="(kg/m2)/s" />
+	<field id="dqbsphy2d"    long_name="Physics dQBS"    unit="(kg/m2)/s" />
         <field id="dqdyn2d"    long_name="Dynamics dQ"    unit="(kg/m2)/s" />
         <field id="dqldyn2d"   long_name="Dynamics dQL"    unit="(kg/m2)/s" />
         <field id="dqsdyn2d"   long_name="Dynamics dQS"    unit="(kg/m2)/s" />
-        <field id="dqcon2d"    long_name="Convection dQ"    unit="(kg/m2)/s" />
+        <field id="dqbsdyn2d"   long_name="Dynamics dQBS"    unit="(kg/m2)/s" />
+	<field id="dqcon2d"    long_name="Convection dQ"    unit="(kg/m2)/s" />
         <field id="dqwak2d"    long_name="Wake dQ"    unit="(kg/m2)/s" />
         <field id="dqlsc2d"    long_name="Condensation dQ"    unit="(kg/m2)/s" />
@@ -590,5 +597,6 @@
         <field id="ovapinit"    long_name="Specific humidity (begin of timestep)"    unit="kg/kg" />
         <field id="oliq"    long_name="Liquid water"    unit="kg/kg" />
-        <field id="ocond"    long_name="Condensed water"    unit="kg/kg" />
+	<field id="ocond"    long_name="Condensed water"    unit="kg/kg" />
+	<field id="qbs"    long_name="specific content of blowing snow"    unit="kg/kg" />
         <field id="geop"    long_name="Geopotential height"    unit="m2/s2" />
         <field id="geoph" field_ref="geop" long_name="Geopotential height" unit="m" > geop/9.80665 </field>
@@ -616,9 +624,11 @@
         <field id="dqphy"    long_name="Physics dQ"    unit="(kg/kg)/s" />
         <field id="dqlphy"   long_name="Physics dQL"    unit="(kg/kg)/s" />
-        <field id="dqsphy"   long_name="Physics dQS"    unit="(kg/kg)/s" />
+	<field id="dqsphy"   long_name="Physics dQS"    unit="(kg/kg)/s" />
+	<field id="dqbsphy"   long_name="Physics dQBS"    unit="(kg/kg)/s" />
         <field id="pr_con_l"    long_name="Convective precipitation lic"    unit="-" />
         <field id="pr_con_i"    long_name="Convective precipitation ice"    unit="-" />
         <field id="pr_lsc_l"    long_name="Large scale precipitation lic"    unit="-" />
-        <field id="pr_lsc_i"    long_name="Large scale precipitation ice"    unit="-" />
+	<field id="pr_lsc_i"    long_name="Large scale precipitation ice"    unit="-" />
+	<field id="pr_bs"    long_name="precipitation blowing snow"    unit="-" />
         <field id="re"    long_name="Cloud droplet effective radius"    unit="um" />
         <field id="fl"    long_name="Denominator of Cloud droplet effective radius"    unit="-" />
@@ -643,5 +653,6 @@
         <field id="dqldyn"    long_name="Dynamics dQL"    unit="(kg/kg)/s" />
         <field id="dqsdyn"    long_name="Dynamics dQS"    unit="(kg/kg)/s" />
-        <field id="dudyn"    long_name="Dynamics dU"    unit="m/s2" />
+        <field id="dqbsdyn"    long_name="Dynamics dQBS"    unit="(kg/kg)/s" />       
+	<field id="dudyn"    long_name="Dynamics dU"    unit="m/s2" />
         <field id="dvdyn"    long_name="Dynamics dV"    unit="m/s2" />
         <field id="dtcon"    long_name="Convection dT"    unit="K/s" />
@@ -673,5 +684,9 @@
         <field id="dqvdf_w" long_name="dqvdf within_wake" unit="(kg/kg)/s" />
         <field id="dteva"    long_name="Reevaporation dT"    unit="K/s" />
-        <field id="dqeva"    long_name="Reevaporation dQ"    unit="(kg/kg)/s" />
+	<field id="dqeva"    long_name="Reevaporation dQ"    unit="(kg/kg)/s" />
+        <field id="dqbsvdf"    long_name="Boundary-layer dQBS"    unit="(kg/kg)/s" />
+        <field id="dqbsbs"    long_name="Blowing snow dQBS"    unit="(kg/kg)/s" />
+        <field id="dqbs"    long_name="Blowing snow dQ"    unit="(kg/kg)/s" />
+        <field id="dtbs"    long_name="Blowing snow dT"    unit="(K)/s" />
         <field id="dqlscth"    long_name="dQ therm."    unit="(kg/kg)/s" />
         <field id="dqlscst"    long_name="dQ strat."    unit="(kg/kg)/s" />
Index: /LMDZ6/branches/blowing_snow/DefLists/file_def_histhf_lmdz.xml
===================================================================
--- /LMDZ6/branches/blowing_snow/DefLists/file_def_histhf_lmdz.xml	(revision 4484)
+++ /LMDZ6/branches/blowing_snow/DefLists/file_def_histhf_lmdz.xml	(revision 4485)
@@ -2,5 +2,5 @@
 <file_definition>
     <file_group id="defile">
-        <file id="histhf" name="histhf" output_freq="3h" output_level="_AUTO_" enabled="_AUTO_" compression_level="2" sync_freq="3h">
+        <file id="histhf" name="histhf" output_freq="1h" output_level="_AUTO_" enabled="_AUTO_" compression_level="2" sync_freq="1h">
             
             <!-- VARS 1D -->
@@ -12,16 +12,16 @@
         
             <!-- VARS 2D -->
-            <field_group group_ref="remap_3h" >
+            <field_group group_ref="remap_1h" >
               <field_group grid_ref="grid_out" > 
-                <field field_ref="phis" level="10" />
+                <field field_ref="phis" level="5" />
                 <field field_ref="aire" level="10" operation="once" />
                 <field field_ref="contfracATM" level="5" operation="once" />
                 <field field_ref="contfracOR" level="5" />
-                <field field_ref="flat" level="10" />
-                <field field_ref="slp" level="10" />
+                <field field_ref="flat" level="5" />
+                <field field_ref="slp" level="5" />
                 <field field_ref="sza" level="10" />
                 <field field_ref="ptstar" level="10" />
                 <field field_ref="pt0" level="10" />
-                <field field_ref="tsol" level="10" />
+                <field field_ref="tsol" level="5" />
                 <field field_ref="t2m" level="5" />
                 <field field_ref="t2m_min" level="10" operation="minimum" />
@@ -32,10 +32,10 @@
                 <field field_ref="t2m_sic" level="10" />
                 <field field_ref="gusts" level="10" />
-                <field field_ref="wind10m" level="10" />
+                <field field_ref="wind10m" level="5" />
                 <field field_ref="wind100m" level="10" />
                 <field field_ref="wind10max" level="10" operation="maximum" />
                 <field field_ref="sicf" level="10" />
                 <field field_ref="q2m" level="5" />
-                <field field_ref="ustar" level="10" />
+                <field field_ref="ustar" level="5" />
                 <field field_ref="u10m" level="5" />
                 <field field_ref="v10m" level="5" />
@@ -67,11 +67,15 @@
                 <field field_ref="plun" level="10" />
                 <field field_ref="pluc" level="10" />
-                <field field_ref="snow" level="5" />
-                <field field_ref="evap" level="10" />
-                <field field_ref="evap_ter" level="10" />
+		<field field_ref="snow" level="5" />
+		<field field_ref="bsfall" level="5" />
+		<field field_ref="evap" level="5" />
+		<field field_ref="snowerosion" level="5" />
+		<field field_ref="rhosnow_lic" level="5" />
+                <field field_ref="ustart_lic" level="5" />
+		<field field_ref="evap_ter" level="10" />
                 <field field_ref="evap_lic" level="10" />
                 <field field_ref="evap_oce" level="10" />
                 <field field_ref="evap_sic" level="10" />
-                <field field_ref="msnow" level="10" />
+                <field field_ref="msnow" level="5" />
                 <field field_ref="fsnow" level="10" />
                 <field field_ref="tops" level="10" />
@@ -99,9 +103,9 @@
                 <field field_ref="SWupSFC" level="10" />
                 <field field_ref="SWupSFCclr" level="10" />
-                <field field_ref="SWdnSFC" level="10" />
+                <field field_ref="SWdnSFC" level="5" />
                 <field field_ref="SWdnSFCclr" level="10" />
                 <field field_ref="LWupSFC" level="10" />
                 <field field_ref="LWupSFCclr" level="10" />
-                <field field_ref="LWdnSFC" level="10" />
+                <field field_ref="LWdnSFC" level="5" />
                 <field field_ref="LWdnSFCclr" level="10" />
                 <field field_ref="bils" level="10" />
@@ -112,5 +116,5 @@
                 <field field_ref="bils_enthalp" level="10" />
                 <field field_ref="bils_latent" level="10" />
-                <field field_ref="sens" level="10" />
+                <field field_ref="sens" level="5" />
                 <field field_ref="sens_x" level="11" />
                 <field field_ref="sens_w" level="11" />
@@ -198,5 +202,5 @@
                 <field field_ref="cldm" level="10" />
                 <field field_ref="cldh" level="10" />
-                <field field_ref="cldt" level="10" />
+                <field field_ref="cldt" level="5" />
                 <field field_ref="cldq" level="10" />
                 <field field_ref="lwp" level="10" />
@@ -214,7 +218,8 @@
                 <field field_ref="wbeff" level="10" />
                 <field field_ref="convoccur" level="10" />
-                <field field_ref="prw" level="10" />
+                <field field_ref="prw" level="5" />
                 <field field_ref="prlw" level="10" />
-                <field field_ref="prsw" level="10" />
+		<field field_ref="prsw" level="10" />
+                <field field_ref="prbsw" level="5" />
                 <field field_ref="s_pblh" level="10" />
                 <field field_ref="s_pblt" level="10" />
@@ -246,5 +251,5 @@
                     <field field_ref="u850" level="10" />
                     <field field_ref="u700" level="10" />
-                    <field field_ref="u500" level="10" />
+                    <field field_ref="u500" level="5" />
                     <field field_ref="u200" level="10" />
                     <field field_ref="u100" level="10" />
@@ -253,5 +258,5 @@
                     <field field_ref="v850" level="10" />
                     <field field_ref="v700" level="10" />
-                    <field field_ref="v500" level="10" />
+                    <field field_ref="v500" level="5" />
                     <field field_ref="v200" level="10" />
                     <field field_ref="v100" level="10" />
@@ -267,5 +272,5 @@
                     <field field_ref="t850" level="10" />
                     <field field_ref="t700" level="10" />
-                    <field field_ref="t500" level="10" />
+                    <field field_ref="t500" level="5" />
                     <field field_ref="t200" level="10" />
                     <field field_ref="t100" level="10" />
@@ -281,5 +286,5 @@
                     <field field_ref="z850" level="10" />
                     <field field_ref="z700" level="10" />
-                    <field field_ref="z500" level="10" />
+                    <field field_ref="z500" level="5" />
                     <field field_ref="z200" level="10" />
                     <field field_ref="z100" level="10" />
@@ -299,5 +304,5 @@
                 <field field_ref="v10_oce" level="10" />
                 <field field_ref="v10_sic" level="10" />
-                <field field_ref="rh2m" level="10" />
+                <field field_ref="rh2m" level="5" />
                 <field field_ref="rh2m_min" level="10" operation="minimum" />
                 <field field_ref="rh2m_max" level="10" operation="maximum" />
@@ -503,21 +508,22 @@
                 <field field_ref="concno3" level="10" />
                 <field field_ref="ec550aer" level="10" />
-                <field field_ref="lwcon" level="10" />
-                <field field_ref="iwcon" level="10" />
+                <field field_ref="lwcon" level="5" />
+                <field field_ref="iwcon" level="5" />
                 <field field_ref="temp" level="4" />
                 <field field_ref="theta" level="10" />
-                <field field_ref="ovap" level="10" />
+                <field field_ref="ovap" level="5" />
                 <field field_ref="ovapinit" level="10" />
-                <field field_ref="oliq" level="10" />
-                <field field_ref="ocond" level="10" />
+                <field field_ref="oliq" level="5" />
+		<field field_ref="ocond" level="5" />
+		<field field_ref="qbs" level="5" />
                 <field field_ref="geop" level="10" />
                 <field field_ref="geoph" level="10" />
-                <field field_ref="vitu" level="10" />
-                <field field_ref="vitv" level="10" />
+                <field field_ref="vitu" level="5" />
+                <field field_ref="vitv" level="5" />
                 <field field_ref="vitw" level="6" />
-                <field field_ref="pres" level="10" />
+                <field field_ref="pres" level="5" />
                 <field field_ref="paprs" level="10" />
                 <field field_ref="mass" level="10" />
-                <field field_ref="zfull" level="10" />
+                <field field_ref="zfull" level="5" />
                 <field field_ref="zhalf" level="10" />
                 <field field_ref="rneb" level="10" />
@@ -525,11 +531,12 @@
                 <field field_ref="rnebls" level="10" />
                 <field field_ref="rneblsvol" level="10" />
-                <field field_ref="rhum" level="10" />
+		<field field_ref="rhum" level="10" />
+		<field field_ref="rhi" level="5" />
                 <field field_ref="ozone" level="10" />
                 <field field_ref="upwd" level="10" />
                 <field field_ref="ep" level="10" />
                 <field field_ref="duphy" level="10" />
-                <field field_ref="dtphy" level="10" />
-                <field field_ref="dqphy" level="10" />
+                <field field_ref="dtphy" level="5" />
+                <field field_ref="dqphy" level="5" />
                 <field field_ref="dqlphy" level="10" />
                 <field field_ref="dqsphy" level="10" />
@@ -537,5 +544,6 @@
                 <field field_ref="pr_con_i" level="10" />
                 <field field_ref="pr_lsc_l" level="10" />
-                <field field_ref="pr_lsc_i" level="10" />
+		<field field_ref="pr_lsc_i" level="10" />
+		<field field_ref="pr_bs" level="5" />
                 <field field_ref="re" level="10" />
                 <field field_ref="fl" level="10" />
@@ -555,6 +563,6 @@
                 <field field_ref="dnwd0" level="10" />
                 <field field_ref="mc" level="10" />
-                <field field_ref="dtdyn" level="10" />
-                <field field_ref="dqdyn" level="10" />
+                <field field_ref="dtdyn" level="5" />
+                <field field_ref="dqdyn" level="5" />
                 <field field_ref="dqldyn" level="10" />
                 <field field_ref="dqsdyn" level="10" />
@@ -576,5 +584,7 @@
                 <field field_ref="dtlsc" level="10" />
                 <field field_ref="dtlschr" level="10" />
-                <field field_ref="dqlsc" level="10" />
+		<field field_ref="dqlsc" level="10" />
+		<field field_ref="dtbs" level="5" />
+                <field field_ref="dqbs" level="5" />
                 <field field_ref="beta_prec" level="10" />
                 <field field_ref="dtvdf" level="10" />
Index: /LMDZ6/branches/blowing_snow/libf/dyn3d_common/infotrac.F90
===================================================================
--- /LMDZ6/branches/blowing_snow/libf/dyn3d_common/infotrac.F90	(revision 4484)
+++ /LMDZ6/branches/blowing_snow/libf/dyn3d_common/infotrac.F90	(revision 4485)
@@ -37,5 +37,5 @@
 !  |--------------------+-----------------------+-----------------+---------------+----------------------------|
 !  | water in different |    water tagging      |  water isotopes | other tracers | additional tracers moments |
-!  | phases: H2O_[gls]  |      isotopes         |                 |               |  for higher order schemes  |
+!  | phases: H2O_[glsb] |      isotopes         |                 |               |  for higher order schemes  |
 !  |--------------------+-----------------------+-----------------+---------------+----------------------------|
 !  |                    |                       |                 |               |                            |
@@ -66,5 +66,5 @@
 !  | longName    | Long name (with adv. scheme suffix) for outputs      | ttext       |                        |
 !  | type        | Type (so far: tracer or tag)                         | /           | tracer,tag             |
-!  | phase       | Phases list ("g"as / "l"iquid / "s"olid)             | /           | [g][l][s]              |
+!  | phase       | Phases list ("g"as / "l"iquid / "s"olid / "b"lowing) | /           | [g][l][s][b]           |
 !  | component   | Name(s) of the merged/cumulated section(s)           | /           | coma-separated names   |
 !  | iGeneration | Generation (>=1)                                     | /           |                        |
@@ -92,5 +92,5 @@
 !  | trac   | ntiso  | Isotopes + tagging tracers list + number         | / | ntraciso       |                 |
 !  | zone   | nzone  | Geographic tagging zones   list + number         | / | ntraceurs_zone |                 |
-!  | phase  | nphas  | Phases                     list + number         |                    | [g][l][s], 1:3  |
+!  | phase  | nphas  | Phases                     list + number         |                    |[g][l][s][b],1:4 |
 !  | iqIsoPha        | Index in "qx"           = f(name(1:ntiso)),phas) | iqiso              | 1:nqtot         |
 !  | itZonIso        | Index in "trac(1:ntiso)"= f(zone, name(1:niso))  | index_trac         | 1:ntiso         |
@@ -227,5 +227,5 @@
       nbtr = nqINCA + nqCO2                                          !--- Number of tracers passed to phytrac
       nqtrue = nbtr + nqo                                            !--- Total number of "true" tracers
-      IF(ALL([2,3] /= nqo)) CALL abort_gcm(modname, 'Only 2 or 3 water phases allowed ; found nqo='//TRIM(int2str(nqo)), 1)
+      IF(ALL([2,3,4,5] /= nqo)) CALL abort_gcm(modname, 'Only 2, 3, 4 , 5 water phases allowed ; found nqo='//TRIM(int2str(nqo)), 1)
       ALLOCATE(hadv(nqtrue), hadv_inca(nqINCA), conv_flg_inca(nqINCA), solsym_inca(nqINCA))
       ALLOCATE(vadv(nqtrue), vadv_inca(nqINCA),  pbl_flg_inca(nqINCA))
Index: /LMDZ6/branches/blowing_snow/libf/misc/readTracFiles_mod.f90
===================================================================
--- /LMDZ6/branches/blowing_snow/libf/misc/readTracFiles_mod.f90	(revision 4484)
+++ /LMDZ6/branches/blowing_snow/libf/misc/readTracFiles_mod.f90	(revision 4485)
@@ -114,9 +114,9 @@
   !--- SOME PARAMETERS THAT ARE NOT LIKELY TO CHANGE OFTEN
   CHARACTER(LEN=maxlen), SAVE      :: tran0        = 'air'      !--- Default transporting fluid
-  CHARACTER(LEN=maxlen), PARAMETER :: old_phases   = 'vlir'     !--- Old phases for water (no separator)
-  CHARACTER(LEN=maxlen), PARAMETER :: known_phases = 'glsr'     !--- Known phases initials
+  CHARACTER(LEN=maxlen), PARAMETER :: old_phases   = 'vlirb'     !--- Old phases for water (no separator)
+  CHARACTER(LEN=maxlen), PARAMETER :: known_phases = 'glsrb'     !--- Known phases initials
   INTEGER, PARAMETER :: nphases = LEN_TRIM(known_phases)        !--- Number of phases
   CHARACTER(LEN=maxlen), SAVE      :: phases_names(nphases) &   !--- Known phases names
-                                             = ['gaseous', 'liquid ', 'solid  ', 'cloud  ']
+                                = ['gaseous', 'liquid ', 'solid  ', 'cloud  ','blosno ']
   CHARACTER(LEN=1), SAVE :: phases_sep  =  '_'                  !--- Phase separator
   LOGICAL,          SAVE :: tracs_merge = .TRUE.                !--- Merge/stack tracers lists
Index: /LMDZ6/branches/blowing_snow/libf/phylmd/add_phys_tend_mod.F90
===================================================================
--- /LMDZ6/branches/blowing_snow/libf/phylmd/add_phys_tend_mod.F90	(revision 4484)
+++ /LMDZ6/branches/blowing_snow/libf/phylmd/add_phys_tend_mod.F90	(revision 4485)
@@ -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)
   ! ======================================================================
   ! Ajoute les tendances de couche limite, soit determinees par la
@@ -49,5 +49,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)
@@ -76,7 +76,7 @@
     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)
   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)
   END IF
 
@@ -87,5 +87,5 @@
 ! $Id: add_phys_tend.F90 2611 2016-08-03 15:41:26Z jyg $
 !
-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)
 !======================================================================
@@ -102,11 +102,11 @@
 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
 USE phys_state_var_mod, ONLY: ftsol
 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
 IMPLICIT none
   include "YOMCST.h"
@@ -116,5 +116,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
@@ -136,5 +136,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
@@ -163,7 +163,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)
 !
@@ -172,7 +174,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)
 
@@ -213,4 +216,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(:,:)
@@ -236,7 +240,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)
@@ -250,4 +254,5 @@
      ql_seri(:,:)=ql_seri(:,:)+zdql(:,:)
      qs_seri(:,:)=qs_seri(:,:)+zdqi(:,:)
+     qbs_seri(:,:)=qbs_seri(:,:)+zdqbs(:,:)
 
 !======================================================================
@@ -432,4 +437,5 @@
         ql_seri(:,:)=ql_seri(:,:)-zdql(:,:)
         qs_seri(:,:)=qs_seri(:,:)-zdqi(:,:)
+        qbs_seri(:,:)=qbs_seri(:,:)-zdqbs(:,:)
       ENDIF
 
@@ -446,7 +452,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))
 
     ! ------------------------------------------------
@@ -457,5 +463,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
@@ -465,4 +472,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 
@@ -476,4 +484,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(:,:)
@@ -484,6 +493,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 
@@ -501,6 +510,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"
@@ -511,7 +520,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
@@ -520,5 +529,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
 
 
@@ -543,5 +552,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)
 !
@@ -550,7 +560,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)
 
@@ -585,7 +596,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)
@@ -600,4 +611,5 @@
      ql_n(:,:)=ql(:,:)+zdql(:,:)
      qs_n(:,:)=qs(:,:)+zdqs(:,:)
+     qbs_n(:,:)=qbs(:,:)+zdqbs(:,:)
      temp_n(:,:)=temp(:,:)+zdt(:,:)
 
@@ -616,7 +628,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))
 
     ! ------------------------------------------------
@@ -627,5 +639,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
@@ -639,4 +652,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 
@@ -649,7 +663,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
@@ -658,14 +672,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
 
@@ -677,4 +691,5 @@
     zql_col(:) = 0.
     zqs_col(:) = 0.
+    zqbs_col(:) = 0.
     zek_col(:) = 0.
     zh_dair_col(:) = 0.
@@ -682,4 +697,5 @@
     zh_ql_col(:) = 0.
     zh_qs_col(:) = 0.
+    zh_qbs_col(:) = 0.
 
 !JLD    write (*,*) "rcpd, zcpvap, zcwat, zcice ",rcpd, zcpvap, zcwat, zcice
@@ -694,4 +710,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)
@@ -702,8 +719,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
@@ -720,12 +738,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
@@ -764,4 +782,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)
@@ -800,6 +821,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/blowing_snow/libf/phylmd/blowing_snow_ini_mod.F90
===================================================================
--- /LMDZ6/branches/blowing_snow/libf/phylmd/blowing_snow_ini_mod.F90	(revision 4485)
+++ /LMDZ6/branches/blowing_snow/libf/phylmd/blowing_snow_ini_mod.F90	(revision 4485)
@@ -0,0 +1,72 @@
+module blowing_snow_ini_mod
+
+implicit none
+
+save
+   integer :: prt_level=0,lunout
+   real RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT,RD,RG
+   real coef_eva_bs, expo_eva_bs, fallv_bs, zeta_bs
+   real prt_bs, pbst_bs, qbst_bs
+   integer :: niter_bs
+
+   !$OMP THREADPRIVATE(prt_level, lunout)
+   !$OMP THREADPRIVATE(RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT,RD,RG)
+   !$OMP THREADPRIVATE(coef_eva_bs, expo_eva_bs, fallv_bs, zeta_bs)
+   !$OMP THREADPRIVATE(pbst_bs, prt_bs, qbst_bs)
+   !$OMP THREADPRIVATE(niter_bs)
+
+      contains
+
+      subroutine blowing_snow_ini(prt_level_in,lunout_in, &
+                                  RCPD_in, RLSTT_in, RLVTT_in, RLMLT_in,&
+                                  RVTMP2_in, RTT_in,RD_in,RG_in)
+
+         USE ioipsl_getin_p_mod, ONLY : getin_p
+
+         integer, intent(in) :: prt_level_in,lunout_in
+         real, intent(in) :: RCPD_in, RLSTT_in, RLVTT_in, RLMLT_in
+         real, intent(in) ::  RVTMP2_in, RTT_in, RD_in, RG_in
+
+
+         prt_level=prt_level_in
+         RG=RG_in
+         RD=RD_in
+         RCPD=RCPD_in
+         RLVTT=RLVTT_in
+         RLSTT=RLSTT_in
+         RLMLT=RLMLT_in
+         RTT=RTT_in
+         RG=RG_in
+         RVTMP2=RVTMP2_in
+         lunout=lunout_in
+
+
+         qbst_bs= 0.001
+         CALL getin_p('qbst_bs',qbst_bs)
+
+         pbst_bs= 0.01
+         CALL getin_p('pbst_bs',pbst_bs)
+
+         prt_bs= 0.0003
+         CALL getin_p('prt_bs',prt_bs)
+
+
+         zeta_bs= 3.
+         CALL getin_p('zeta_bs',zeta_bs)
+
+         fallv_bs = 0.1
+         CALL getin_p('fallv_bs',fallv_bs)
+
+         coef_eva_bs =  2.e-5 
+         CALL getin_p('coef_eva_bs',coef_eva_bs)
+
+         expo_eva_bs = 0.5
+         CALL getin_p('expo_eva_bs',expo_eva_bs)
+
+         niter_bs = 5
+         CALL getin_p('niter_bs',niter_bs)
+
+
+      end subroutine blowing_snow_ini
+
+end module blowing_snow_ini_mod
Index: /LMDZ6/branches/blowing_snow/libf/phylmd/blowing_snow_sublim_sedim.F90
===================================================================
--- /LMDZ6/branches/blowing_snow/libf/phylmd/blowing_snow_sublim_sedim.F90	(revision 4485)
+++ /LMDZ6/branches/blowing_snow/libf/phylmd/blowing_snow_sublim_sedim.F90	(revision 4485)
@@ -0,0 +1,212 @@
+subroutine blowing_snow_sublim_sedim(ngrid,nlay,dtime,temp,q,qbs,pplay,paprs,dtemp_bs,dq_bs,dqbs_bs,bsfl,precip_bs)
+
+!==============================================================================
+! Routine that calculates the evaporation and sedimentation of blowing snow
+! inspired by what is done in lscp_mod
+! Etienne Vignon, October 2022
+!==============================================================================
+
+
+use blowing_snow_ini_mod, only : coef_eva_bs,RTT,RD,RG,expo_eva_bs, fallv_bs
+use blowing_snow_ini_mod, only : RCPD, RLSTT, RLMLT, RLVTT, RVTMP2, niter_bs
+USE lscp_tools_mod, only : calc_qsat_ecmwf
+
+implicit none
+
+
+!++++++++++++++++++++++++++++++++++++++++++++++++++++
+! Declarations
+!++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+!INPUT
+!=====
+
+integer, intent(in)                     :: ngrid,nlay
+real, intent(in)                        :: dtime
+real, intent(in), dimension(ngrid,nlay) :: temp
+real, intent(in), dimension(ngrid,nlay) :: q
+real, intent(in), dimension(ngrid,nlay) :: qbs
+real, intent(in), dimension(ngrid,nlay) :: pplay
+real, intent(in), dimension(ngrid,nlay+1) :: paprs
+
+
+
+! OUTPUT
+!========
+
+
+real, intent(out), dimension(ngrid,nlay) :: dtemp_bs
+real, intent(out), dimension(ngrid,nlay) :: dq_bs
+real, intent(out), dimension(ngrid,nlay) :: dqbs_bs
+real, intent(out), dimension(ngrid,nlay+1) :: bsfl
+real, intent(out), dimension(ngrid)      :: precip_bs
+
+
+! LOCAL
+!======
+
+
+integer                                  :: k,i,n
+real                                     :: zqev0, zqevi, zqevti, zcpair, zcpeau, dqbsmelt
+real                                     :: dqsedim,precbs
+real, dimension(ngrid)                   :: zt,zq,zqbs,qsi,dqsi,qsl, dqsl,qzero,sedim,sedimn
+real, dimension(ngrid)                   :: zqbsi,zmqc, zmair, zdz
+real, dimension(ngrid,nlay)              :: velo, zrho
+
+!++++++++++++++++++++++++++++++++++++++++++++++++++
+! Initialisation
+!++++++++++++++++++++++++++++++++++++++++++++++++++
+
+qzero(:)=0.
+dtemp_bs(:,:)=0.
+dq_bs(:,:)=0.
+dqbs_bs(:,:)=0.
+velo(:,:)=0.
+zt(:)=0.
+zq(:)=0.
+zqbs(:)=0.
+sedim(:)=0.
+
+! begin of top-down loop
+DO k = nlay, 1, -1
+    
+    DO i=1,ngrid
+        zt(i)=temp(i,k)
+        zq(i)=q(i,k)
+        zqbs(i)=qbs(i,k)
+    ENDDO
+
+
+    IF (k.LE.nlay-1) THEN
+
+        ! thermalization of blowing snow precip coming from above    
+        DO i = 1, ngrid
+            zmair(i)=(paprs(i,k)-paprs(i,k+1))/RG
+            ! RVTMP2=rcpv/rcpd-1
+            zcpair=RCPD*(1.0+RVTMP2*zq(i))
+            zcpeau=RCPD*RVTMP2
+            ! zmqc: precipitation mass that has to be thermalized with 
+            ! layer's air so that precipitation at the ground has the
+            ! same temperature as the lowermost layer
+            zmqc(i) = (sedim(i))*dtime/zmair(i)
+            ! t(i,k+1)+d_t(i,k+1): new temperature of the overlying layer
+            zt(i) = ( (temp(i,k+1)+dtemp_bs(i,k+1))*zmqc(i)*zcpeau + zcpair*zt(i) ) &
+             / (zcpair + zmqc(i)*zcpeau)
+        ENDDO
+    ELSE
+        DO i = 1, ngrid
+            zmair(i)=(paprs(i,k)-paprs(i,k+1))/RG
+            zmqc(i) = 0.
+        ENDDO
+
+    ENDIF
+
+
+
+    ! calulation saturation specific humidity
+    CALL CALC_QSAT_ECMWF(ngrid,zt(:),qzero(:),pplay(:,k),RTT,2,.false.,qsi(:),dqsi(:))
+    CALL CALC_QSAT_ECMWF(ngrid,zt(:),qzero(:),pplay(:,k),RTT,1,.false.,qsl(:),dqsl(:))
+    ! sublimation calculation
+    ! SUndqvist formula dP/dz=beta*(1-q/qsat)*sqrt(P)
+    
+    DO i = 1, ngrid
+       ! if sedimentation:
+       IF (sedim(i) .GT. 0.) THEN
+
+          IF (zt(i) .GT. RTT) THEN
+             ! if positive celcius temperature, we assume
+             ! that all the blowing snow melt and evaporate
+              zqevti=sedim(i)*(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG     
+             ! we ensure that the whole mesh does not exceed saturation wrt liquid
+              zqev0 = MAX(0.0, qsl(i)-zq(i))
+              zqevi = MIN(zqev0,zqevti)
+              ! New solid precipitation fluxes
+              sedimn(i) = Max(0.,sedim(i) - zqevi*(paprs(i,k)-paprs(i,k+1))/RG/dtime)
+              ! vapor, temperature, precip fluxes update
+              zq(i) = zq(i) - (sedimn(i)-sedim(i))*(RG/(paprs(i,k)-paprs(i,k+1)))*dtime
+              zqbs(i) = zqbs(i) + (sedimn(i)-sedim(i)) * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
+              ! melting
+              zt(i) = zt(i)                             &
+                + (sedimn(i)-sedim(i))                      &
+                * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime    &
+                * RLMLT/RCPD/(1.0+RVTMP2*(zq(i)+zqbs(i)))
+              ! evaporation
+              zt(i) = zt(i)                             &
+                + (sedimn(i)-sedim(i))                      &
+                * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime    &
+                * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zqbs(i)))      
+      
+      
+          ELSE     
+              zqevti = coef_eva_bs*(1.0-zq(i)/qsi(i))*(sedim(i)**expo_eva_bs) &
+                    *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
+              zqevti = MAX(0.0,MIN(zqevti,sedim(i)))*RG*dtime/(paprs(i,k)-paprs(i,k+1))
+
+              ! Sublimation limit: we ensure that the whole mesh does not exceed saturation wrt ice
+              zqev0 = MAX(0.0, qsi(i)-zq(i))
+              zqevi = MIN(zqev0,zqevti)
+
+              ! New solid precipitation fluxes
+              sedimn(i) = Max(0.,sedim(i) - zqevi*(paprs(i,k)-paprs(i,k+1))/RG/dtime)
+
+
+              ! vapor, temperature, precip fluxes update
+              zq(i) = zq(i) - (sedimn(i)-sedim(i))*(RG/(paprs(i,k)-paprs(i,k+1)))*dtime                
+              zqbs(i) = zqbs(i) + (sedimn(i)-sedim(i)) * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
+              zt(i) = zt(i)                             &
+                + (sedimn(i)-sedim(i))                      &
+                * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime    &
+                * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zqbs(i)))
+           ENDIF
+
+
+          ! sedim update
+          sedim(i)=sedimn(i)
+
+
+        ELSE
+           sedim(i)=0.
+        ENDIF ! if sedim > 0
+
+        zqbsi(i)=zqbs(i)
+
+    ENDDO ! loop on ngrid
+
+    ! Now sedimention scheme (alike that in lscp)
+    DO n = 1, niter_bs
+       DO i = 1, ngrid
+         zrho(i,k)  = pplay(i,k) / zt(i) / RD
+         zdz(i)   = (paprs(i,k)-paprs(i,k+1)) / (zrho(i,k)*RG)
+         velo(i,k) = fallv_bs 
+         dqsedim = dtime/REAL(niter_bs)/zdz(i)*zqbs(i)*velo(i,k)   ! dqice/dt=1/rho*d(rho*wice*qice)/dz
+         precbs   = MIN(MAX(dqsedim,0.0),zqbs(i))
+         zqbs(i) = MAX(zqbs(i)-1.*precbs  , 0.0)
+       ENDDO !loop on ngrid
+    ENDDO ! loop on niter_bs
+
+    ! add to non sublimated precip
+    DO i=1,ngrid
+    sedim(i) = sedim(i)+max(0.,zqbsi(i)-zqbs(i))*(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
+    ENDDO
+
+    ! Outputs:
+    DO i = 1, ngrid
+        bsfl(i,k)=sedim(i)
+        dqbs_bs(i,k) = zqbs(i)-qbs(i,k)
+        dq_bs(i,k) = zq(i) - q(i,k)
+        dtemp_bs(i,k) = zt(i) - temp(i,k)
+    ENDDO
+
+
+ENDDO ! vertical loop
+
+
+!surface bs flux
+DO i = 1, ngrid
+  precip_bs(i) = sedim(i)
+ENDDO
+
+
+return
+
+end subroutine blowing_snow_sublim_sedim
Index: /LMDZ6/branches/blowing_snow/libf/phylmd/call_blowing_snow_sublim_sedim.F90
===================================================================
--- /LMDZ6/branches/blowing_snow/libf/phylmd/call_blowing_snow_sublim_sedim.F90	(revision 4485)
+++ /LMDZ6/branches/blowing_snow/libf/phylmd/call_blowing_snow_sublim_sedim.F90	(revision 4485)
@@ -0,0 +1,38 @@
+subroutine call_blowing_snow_sublim_sedim(ngrid,nlay,dtime,temp,q,qbs,pplay,paprs, &
+                                  dtemp_bs,dq_bs,dqbs_bs,bsfl,precip_bs)
+
+implicit none
+
+!INPUT
+!=====
+
+integer, intent(in)                     :: ngrid,nlay
+real, intent(in)                        :: dtime
+real, intent(in), dimension(ngrid,nlay) :: temp
+real, intent(in), dimension(ngrid,nlay) :: q
+real, intent(in), dimension(ngrid,nlay) :: qbs
+real, intent(in), dimension(ngrid,nlay) :: pplay
+real, intent(in), dimension(ngrid,nlay+1) :: paprs
+
+
+! OUTPUT
+!========
+
+
+real, intent(out), dimension(ngrid,nlay) :: dtemp_bs
+real, intent(out), dimension(ngrid,nlay) :: dq_bs
+real, intent(out), dimension(ngrid,nlay) :: dqbs_bs
+real, intent(out), dimension(ngrid,nlay+1) :: bsfl
+real, intent(out), dimension(ngrid)      :: precip_bs
+
+      
+
+
+call blowing_snow_sublim_sedim(ngrid,nlay,dtime,temp,q,qbs,pplay,paprs, &
+                                  dtemp_bs,dq_bs,dqbs_bs,bsfl,precip_bs)
+
+
+
+
+
+end subroutine call_blowing_snow_sublim_sedim
Index: /LMDZ6/branches/blowing_snow/libf/phylmd/clesphys.h
===================================================================
--- /LMDZ6/branches/blowing_snow/libf/phylmd/clesphys.h	(revision 4484)
+++ /LMDZ6/branches/blowing_snow/libf/phylmd/clesphys.h	(revision 4485)
@@ -100,4 +100,5 @@
        LOGICAL :: ok_daily_climoz
        LOGICAL :: ok_new_lscp
+       LOGICAL :: ok_bs, ok_rad_bs
 ! flag to bypass or not the phytrac module
        INTEGER :: iflag_phytrac
@@ -150,5 +151,5 @@
      &     , ok_chlorophyll,ok_conserv_q, adjust_tropopause             &
      &     , ok_daily_climoz, ok_all_xml, ok_lwoff                      &
-     &     , iflag_phytrac, ok_new_lscp                                 &
+     &     , iflag_phytrac, ok_new_lscp, ok_bs, ok_rad_bs               &
      &     ,  iflag_thermals,nsplit_thermals, tau_thermals    
        save /clesphys/
Index: /LMDZ6/branches/blowing_snow/libf/phylmd/climb_qbs_mod.F90
===================================================================
--- /LMDZ6/branches/blowing_snow/libf/phylmd/climb_qbs_mod.F90	(revision 4485)
+++ /LMDZ6/branches/blowing_snow/libf/phylmd/climb_qbs_mod.F90	(revision 4485)
@@ -0,0 +1,389 @@
+MODULE climb_qbs_mod
+!
+! Module to solve the verctical diffusion of blowing snow; 
+!
+  USE dimphy
+
+  IMPLICIT NONE
+  SAVE 
+  PRIVATE
+  PUBLIC :: climb_qbs_down, climb_qbs_up
+
+  REAL, DIMENSION(:,:), ALLOCATABLE :: gamaqbs
+  !$OMP THREADPRIVATE(gamaqbs)
+  REAL, DIMENSION(:,:), ALLOCATABLE :: Ccoef_QBS, Dcoef_QBS
+  !$OMP THREADPRIVATE(Ccoef_QBS, Dcoef_QBS)
+  REAL, DIMENSION(:), ALLOCATABLE   :: Acoef_QBS, Bcoef_QBS
+  !$OMP THREADPRIVATE(Acoef_QBS, Bcoef_QBS)
+  REAL, DIMENSION(:,:), ALLOCATABLE :: Kcoefqbs
+  !$OMP THREADPRIVATE(Kcoefqbs)
+
+CONTAINS
+!
+!****************************************************************************************
+!
+  SUBROUTINE climb_qbs_down(knon, coefqbs, paprs, pplay, &
+       delp, temp, qbs, dtime, &
+       Ccoef_QBS_out, Dcoef_QBS_out, &
+       Kcoef_qbs_out, gama_qbs_out, &
+       Acoef_QBS_out, Bcoef_QBS_out)
+
+! This routine calculates recursivly the coefficients C and D
+! for the quantity X=[QBS] in equation X(k) = C(k) + D(k)*X(k-1), where k is
+! the index of the vertical layer.
+!
+! Input arguments
+!****************************************************************************************
+    INTEGER, INTENT(IN)                      :: knon
+    REAL, DIMENSION(klon,klev), INTENT(IN)   :: coefqbs
+    REAL, DIMENSION(klon,klev), INTENT(IN)   :: pplay 
+    REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs 
+    REAL, DIMENSION(klon,klev), INTENT(IN)   :: delp  
+    REAL, DIMENSION(klon,klev), INTENT(IN)   :: temp
+    REAL, DIMENSION(klon,klev), INTENT(IN)   :: qbs
+    REAL, INTENT(IN)                         :: dtime
+
+! Output arguments
+!****************************************************************************************
+    REAL, DIMENSION(klon), INTENT(OUT)       :: Acoef_QBS_out
+    REAL, DIMENSION(klon), INTENT(OUT)       :: Bcoef_QBS_out
+
+    REAL, DIMENSION(klon,klev), INTENT(OUT)  :: Ccoef_QBS_out
+    REAL, DIMENSION(klon,klev), INTENT(OUT)  :: Dcoef_QBS_out
+    REAL, DIMENSION(klon,klev), INTENT(OUT)  :: Kcoef_qbs_out
+    REAL, DIMENSION(klon,klev), INTENT(OUT)  :: gama_qbs_out
+
+! Local variables
+!****************************************************************************************
+    LOGICAL, SAVE                            :: first=.TRUE.
+    !$OMP THREADPRIVATE(first)
+    REAL, DIMENSION(klon)                    :: psref 
+    REAL                                     :: delz, pkh
+    INTEGER                                  :: k, i, ierr
+! Include
+!****************************************************************************************
+    INCLUDE "YOMCST.h"
+    INCLUDE "compbl.h"    
+
+
+!****************************************************************************************
+! 1)
+! Allocation at first time step only
+!   
+!****************************************************************************************
+
+    IF (first) THEN
+       first=.FALSE.
+       ALLOCATE(Ccoef_QBS(klon,klev), STAT=ierr)
+       IF ( ierr /= 0 )  PRINT*,' pb in allloc Ccoef_QBS, ierr=', ierr
+       
+       ALLOCATE(Dcoef_QBS(klon,klev), STAT=ierr)
+       IF ( ierr /= 0 )  PRINT*,' pb in allloc Dcoef_QBS, ierr=', ierr
+       
+       ALLOCATE(Acoef_QBS(klon), Bcoef_QBS(klon), STAT=ierr)
+       IF ( ierr /= 0 )  PRINT*,' pb in allloc Acoef_BS and Bcoef_BS, ierr=', ierr
+       
+       ALLOCATE(Kcoefqbs(klon,klev), STAT=ierr)
+       IF ( ierr /= 0 )  PRINT*,' pb in allloc Kcoefqbs, ierr=', ierr
+       
+       ALLOCATE(gamaqbs(1:klon,2:klev), STAT=ierr)
+       IF ( ierr /= 0 ) PRINT*,' pb in allloc gamaqbs, ierr=', ierr
+       
+    END IF
+
+!****************************************************************************************
+! 2)
+! Definition of the coeficient K 
+!
+!****************************************************************************************
+    Kcoefqbs(:,:) = 0.0
+    DO k = 2, klev
+       DO i = 1, knon
+          Kcoefqbs(i,k) = &
+               coefqbs(i,k)*RG*RG*dtime /(pplay(i,k-1)-pplay(i,k)) &
+               *(paprs(i,k)*2/(temp(i,k)+temp(i,k-1))/RD)**2
+       ENDDO
+    ENDDO
+
+!****************************************************************************************
+! 3)
+! Calculation of gama for "Q" and "H"
+!
+!****************************************************************************************
+!   surface pressure is used as reference
+    psref(:) = paprs(:,1) 
+
+!   definition of gama
+    IF (iflag_pbl == 1) THEN
+       gamaqbs(:,:) = 0.0
+ 
+! conversion de gama
+       DO k = 2, klev
+          DO i = 1, knon
+             delz = RD * (temp(i,k-1)+temp(i,k)) / & 
+                    2.0 / RG / paprs(i,k) * (pplay(i,k-1)-pplay(i,k))
+             pkh  = (psref(i)/paprs(i,k))**RKAPPA
+          
+! convertie gradient verticale de contenu en neige soufflee en difference de neige soufflee entre centre de couches
+             gamaqbs(i,k) = gamaqbs(i,k) * delz    
+          ENDDO
+       ENDDO
+
+    ELSE
+       gamaqbs(:,:) = 0.0
+    ENDIF
+    
+
+!****************************************************************************************    
+! 4)
+! Calculte the coefficients C and D for specific content of blowing snow, qbs
+!
+!****************************************************************************************
+    
+    CALL calc_coef_qbs(knon, Kcoefqbs(:,:), gamaqbs(:,:), delp(:,:), qbs(:,:), &
+         Ccoef_QBS(:,:), Dcoef_QBS(:,:), Acoef_QBS, Bcoef_QBS)
+
+ 
+!****************************************************************************************
+! 5)
+! Return the first layer in output variables
+!
+!****************************************************************************************
+    Acoef_QBS_out = Acoef_QBS
+    Bcoef_QBS_out = Bcoef_QBS
+
+!****************************************************************************************
+! 6)
+! If Pbl is split, return also the other layers in output variables
+!
+!****************************************************************************************
+    IF (mod(iflag_pbl_split,10) .ge.1) THEN
+    DO k= 1, klev
+      DO i= 1, klon
+        Ccoef_QBS_out(i,k) = Ccoef_QBS(i,k)
+        Dcoef_QBS_out(i,k) = Dcoef_QBS(i,k)
+        Kcoef_qbs_out(i,k) = Kcoefqbs(i,k)
+          IF (k.eq.1) THEN
+            gama_qbs_out(i,k)  = 0.
+          ELSE
+            gama_qbs_out(i,k)  = gamaqbs(i,k)
+          ENDIF
+      ENDDO
+    ENDDO
+!!!      
+       ENDIF  ! (mod(iflag_pbl_split,2) .ge.1)
+!!!
+
+  END SUBROUTINE climb_qbs_down
+!
+!****************************************************************************************
+!
+  SUBROUTINE calc_coef_qbs(knon, Kcoef, gama, delp, X, Ccoef, Dcoef, Acoef, Bcoef)
+!
+! Calculate the coefficients C and D in : X(k) = C(k) + D(k)*X(k-1)
+! where X is QQBS, and k the vertical level k=1,klev
+!
+    INCLUDE "YOMCST.h"
+! Input arguments
+!****************************************************************************************
+    INTEGER, INTENT(IN)                      :: knon
+    REAL, DIMENSION(klon,klev), INTENT(IN)   :: Kcoef, delp
+    REAL, DIMENSION(klon,klev), INTENT(IN)   :: X
+    REAL, DIMENSION(klon,2:klev), INTENT(IN) :: gama
+
+! Output arguments
+!****************************************************************************************
+    REAL, DIMENSION(klon), INTENT(OUT)       :: Acoef, Bcoef
+    REAL, DIMENSION(klon,klev), INTENT(OUT)  :: Ccoef, Dcoef
+
+! Local variables
+!****************************************************************************************
+    INTEGER                                  :: k, i
+    REAL                                     :: buf
+
+!****************************************************************************************
+! Niveau au sommet, k=klev
+!
+!****************************************************************************************
+    Ccoef(:,:) = 0.0
+    Dcoef(:,:) = 0.0
+
+    DO i = 1, knon
+       buf = delp(i,klev) + Kcoef(i,klev)
+       
+       Ccoef(i,klev) = (X(i,klev)*delp(i,klev) - Kcoef(i,klev)*gama(i,klev))/buf
+       Dcoef(i,klev) = Kcoef(i,klev)/buf
+    END DO
+
+
+!****************************************************************************************
+! Niveau  (klev-1) <= k <= 2
+!
+!****************************************************************************************
+
+    DO k=(klev-1),2,-1
+       DO i = 1, knon
+          buf = delp(i,k) + Kcoef(i,k) + Kcoef(i,k+1)*(1.-Dcoef(i,k+1))
+          Ccoef(i,k) = (X(i,k)*delp(i,k) + Kcoef(i,k+1)*Ccoef(i,k+1) + &
+               Kcoef(i,k+1)*gama(i,k+1) - Kcoef(i,k)*gama(i,k))/buf
+          Dcoef(i,k) = Kcoef(i,k)/buf
+       END DO
+    END DO
+
+!****************************************************************************************
+! Niveau k=1
+!
+!****************************************************************************************
+
+    DO i = 1, knon
+       buf = delp(i,1) + Kcoef(i,2)*(1.-Dcoef(i,2))
+       Acoef(i) = (X(i,1)*delp(i,1) + Kcoef(i,2)*(gama(i,2)+Ccoef(i,2)))/buf
+       Bcoef(i) = -1. * RG / buf
+    END DO
+
+  END SUBROUTINE calc_coef_qbs
+!
+!****************************************************************************************
+!
+  SUBROUTINE climb_qbs_up(knon, dtime, qbs_old, &
+       flx_qbs1, paprs, pplay, &
+       Acoef_QBS_in, Bcoef_QBS_in, &
+       Ccoef_QBS_in, Dcoef_QBS_in, &
+       Kcoef_qbs_in, gama_qbs_in, &
+       flux_qbs, d_qbs)
+! 
+! This routine calculates the flux and tendency of the specific content of blowing snow qbs 
+! The quantity qbs is calculated according to 
+! X(k) = C(k) + D(k)*X(k-1) for X=[qbs], where the coefficients 
+! C and D are known from before and k is index of the vertical layer.
+!   
+
+! Input arguments
+!****************************************************************************************
+    INTEGER, INTENT(IN)                      :: knon
+    REAL, INTENT(IN)                         :: dtime
+    REAL, DIMENSION(klon,klev), INTENT(IN)   :: qbs_old
+    REAL, DIMENSION(klon), INTENT(IN)        :: flx_qbs1
+    REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs
+    REAL, DIMENSION(klon,klev), INTENT(IN)   :: pplay
+
+!!! nrlmd le 02/05/2011
+    REAL, DIMENSION(klon), INTENT(IN)        :: Acoef_QBS_in, Bcoef_QBS_in
+    REAL, DIMENSION(klon,klev), INTENT(IN)   :: Ccoef_QBS_in, Dcoef_QBS_in
+    REAL, DIMENSION(klon,klev), INTENT(IN)   :: Kcoef_qbs_in, gama_qbs_in
+!!!
+
+! Output arguments
+!****************************************************************************************
+    REAL, DIMENSION(klon,klev), INTENT(OUT)  :: flux_qbs, d_qbs
+
+! Local variables
+!****************************************************************************************
+    LOGICAL, SAVE                            :: last=.FALSE.
+    REAL, DIMENSION(klon,klev)               :: qbs_new
+    REAL, DIMENSION(klon)                    :: psref         
+    INTEGER                                  :: k, i, ierr
+ 
+! Include
+!****************************************************************************************
+    INCLUDE "YOMCST.h"
+    INCLUDE "compbl.h"    
+
+!****************************************************************************************
+! 1) 
+! Definition of some variables
+    REAL, DIMENSION(klon,klev)               :: zairm
+!
+!****************************************************************************************
+    flux_qbs(:,:) = 0.0
+    d_qbs(:,:)    = 0.0
+
+    psref(1:knon) = paprs(1:knon,1)  
+
+       IF (mod(iflag_pbl_split,10) .ge.1) THEN
+    DO i = 1, knon
+      Acoef_QBS(i)=Acoef_QBS_in(i)
+      Bcoef_QBS(i)=Bcoef_QBS_in(i)
+    ENDDO
+    DO k = 1, klev
+      DO i = 1, knon
+        Ccoef_QBS(i,k)=Ccoef_QBS_in(i,k)
+        Dcoef_QBS(i,k)=Dcoef_QBS_in(i,k)
+        Kcoefqbs(i,k)=Kcoef_qbs_in(i,k)
+          IF (k.gt.1) THEN
+            gamaqbs(i,k)=gama_qbs_in(i,k)
+          ENDIF
+      ENDDO
+    ENDDO
+!!!      
+       ENDIF  ! (mod(iflag_pbl_split,2) .ge.1)
+!!!
+
+!****************************************************************************************
+! 2)
+! Calculation of QBS
+!
+!****************************************************************************************
+
+!- First layer
+    qbs_new(1:knon,1) = Acoef_QBS(1:knon) + Bcoef_QBS(1:knon)*flx_qbs1(1:knon)*dtime
+!- All the other layers 
+    DO k = 2, klev
+       DO i = 1, knon
+          qbs_new(i,k) = Ccoef_QBS(i,k) + Dcoef_QBS(i,k)*qbs_new(i,k-1)
+       END DO
+    END DO
+!****************************************************************************************
+! 3)
+! Calculation of the flux for QBS
+!
+!****************************************************************************************
+
+!- The flux at first layer, k=1
+    flux_qbs(1:knon,1)=flx_qbs1(1:knon)
+
+!- The flux at all layers above surface
+    DO k = 2, klev
+       DO i = 1, knon
+          flux_qbs(i,k) = (Kcoefqbs(i,k)/RG/dtime) * &
+               (qbs_new(i,k)-qbs_new(i,k-1)+gamaqbs(i,k))
+       END DO
+    END DO
+
+!****************************************************************************************
+! 4)
+! Calculation of tendency for QBS
+!
+!****************************************************************************************
+    DO k = 1, klev
+       DO i = 1, knon
+          d_qbs(i,k) = qbs_new(i,k) - qbs_old(i,k)
+          zairm(i, k) = (paprs(i,k)-paprs(i,k+1))/rg
+        END DO
+    END DO
+
+!****************************************************************************************
+! Some deallocations
+!
+!****************************************************************************************
+    IF (last) THEN
+       DEALLOCATE(Ccoef_QBS, Dcoef_QBS,stat=ierr)    
+       IF ( ierr /= 0 )  PRINT*,' pb in dealllocate Ccoef_QBS, Dcoef_QBS, ierr=', ierr
+       DEALLOCATE(Acoef_QBS, Bcoef_QBS,stat=ierr)    
+       IF ( ierr /= 0 )  PRINT*,' pb in dealllocate Acoef_QBS, Bcoef_QBS, ierr=', ierr
+       DEALLOCATE(gamaqbs,stat=ierr)
+       IF ( ierr /= 0 )  PRINT*,' pb in dealllocate gamaqbs, ierr=', ierr
+       DEALLOCATE(Kcoefqbs,stat=ierr)
+       IF ( ierr /= 0 )  PRINT*,' pb in dealllocate Kcoefqbs, ierr=', ierr
+    END IF
+  END SUBROUTINE climb_qbs_up
+!
+!****************************************************************************************
+!
+END MODULE climb_qbs_mod
+
+ 
+
+
+
+
Index: /LMDZ6/branches/blowing_snow/libf/phylmd/conf_phys_m.F90
===================================================================
--- /LMDZ6/branches/blowing_snow/libf/phylmd/conf_phys_m.F90	(revision 4484)
+++ /LMDZ6/branches/blowing_snow/libf/phylmd/conf_phys_m.F90	(revision 4485)
@@ -255,4 +255,5 @@
     LOGICAL, SAVE :: ok_new_lscp_omp
     LOGICAL, SAVE :: ok_icefra_lscp_omp
+    LOGICAL, SAVE :: ok_bs_omp, ok_rad_bs_omp
 
 
@@ -2359,4 +2360,22 @@
     !Config Help = ...
 
+
+
+    ok_bs_omp = .FALSE.
+    CALL getin('ok_bs', ok_bs_omp)
+    !
+    !Config Key  = ok_bs_omp
+    !Config Desc = include blowing snow
+    !Config Def  = .FALSE.
+    !Config Help = ...
+
+
+    ok_rad_bs_omp = .FALSE.
+    CALL getin('ok_rad_bs', ok_rad_bs_omp)
+    !
+    !Config Key  = ok_rad_bs_omp
+    !Config Desc = include blowing snow radiative effect
+    !Config Def  = .FALSE.
+    !Config Help = ...
 
     ecrit_LES_omp = 1./8.
@@ -2473,4 +2492,6 @@
     inertie_lic = inertie_lic_omp
     inertie_sno = inertie_sno_omp
+    ok_bs = ok_bs_omp
+    ok_rad_bs=ok_rad_bs_omp
     rad_froid = rad_froid_omp
     rad_chau1 = rad_chau1_omp
@@ -2982,4 +3003,6 @@
     WRITE(lunout,*) ' inertie_lic = ', inertie_lic
     WRITE(lunout,*) ' inertie_sno = ', inertie_sno
+    WRITE(lunout,*) ' ok_bs = ', ok_bs
+    WRITE(lunout,*) ' ok_rad_bs = ', ok_rad_bs
     WRITE(lunout,*) ' f_cdrag_ter = ',f_cdrag_ter
     WRITE(lunout,*) ' f_cdrag_oce = ',f_cdrag_oce
Index: /LMDZ6/branches/blowing_snow/libf/phylmd/fonte_neige_mod.F90
===================================================================
--- /LMDZ6/branches/blowing_snow/libf/phylmd/fonte_neige_mod.F90	(revision 4484)
+++ /LMDZ6/branches/blowing_snow/libf/phylmd/fonte_neige_mod.F90	(revision 4485)
@@ -164,6 +164,6 @@
     REAL, DIMENSION(klon), INTENT(IN)    :: precip_rain
     REAL, DIMENSION(klon), INTENT(IN)    :: precip_snow
-    
-! Input/Output variables
+
+    ! Input/Output variables
 !****************************************************************************************
 
Index: /LMDZ6/branches/blowing_snow/libf/phylmd/infotrac_phy.F90
===================================================================
--- /LMDZ6/branches/blowing_snow/libf/phylmd/infotrac_phy.F90	(revision 4484)
+++ /LMDZ6/branches/blowing_snow/libf/phylmd/infotrac_phy.F90	(revision 4485)
@@ -68,5 +68,5 @@
 !  | longName    | Long name (with adv. scheme suffix) for outputs      | ttext       |                        |
 !  | type        | Type (so far: tracer or tag)                         | /           | tracer,tag             |
-!  | phase       | Phases list ("g"as / "l"iquid / "s"olid)             | /           | [g][l][s]              |
+!  | phase       | Phases list ("g"as / "l"iquid / "s"olid / "b"lowing) | /           | [g][l][s][b]           |
 !  | component   | Name(s) of the merged/cumulated section(s)           | /           | coma-separated names   |
 !  | iGeneration | Generation (>=1)                                     | /           |                        |
@@ -92,5 +92,5 @@
 !  | trac   | ntiso  | Isotopes + tagging tracers list + number         | / | ntraciso       |                 |
 !  | zone   | nzone  | Geographic tagging zones   list + number         | / | ntraceurs_zone |                 |
-!  | phase  | nphas  | Phases                     list + number         |                    | [g][l][s], 1:3  |
+!  | phase  | nphas  | Phases                     list + number         |                    |[g][l][s][b] 1:4 |
 !  | iqIsoPha        | Index in "qx"           = f(name(1:ntiso)),phas) | iqiso              | 1:nqtot         |
 !  | itZonIso        | Index in "trac(1:ntiso)"= f(zone, name(1:niso))  | index_trac         | 1:ntiso         |
Index: /LMDZ6/branches/blowing_snow/libf/phylmd/ocean_forced_mod.F90
===================================================================
--- /LMDZ6/branches/blowing_snow/libf/phylmd/ocean_forced_mod.F90	(revision 4484)
+++ /LMDZ6/branches/blowing_snow/libf/phylmd/ocean_forced_mod.F90	(revision 4485)
@@ -167,5 +167,5 @@
   SUBROUTINE ocean_forced_ice( &
        itime, dtime, jour, knon, knindex, &
-       tsurf_in, p1lay, cdragh, cdragm, precip_rain, precip_snow, temp_air, spechum, &
+       tsurf_in, p1lay, cdragh, cdragm, precip_rain, precip_snow, temp_air,spechum, &
        AcoefH, AcoefQ, BcoefH, BcoefQ, &
        AcoefU, AcoefV, BcoefU, BcoefV, &
Index: /LMDZ6/branches/blowing_snow/libf/phylmd/pbl_surface_mod.F90
===================================================================
--- /LMDZ6/branches/blowing_snow/libf/phylmd/pbl_surface_mod.F90	(revision 4484)
+++ /LMDZ6/branches/blowing_snow/libf/phylmd/pbl_surface_mod.F90	(revision 4485)
@@ -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
@@ -184,7 +185,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,           &
@@ -198,8 +199,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,                              &
@@ -224,5 +225,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,                    &
@@ -230,5 +231,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
@@ -317,4 +318,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
@@ -351,4 +353,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
@@ -356,4 +359,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
@@ -420,4 +424,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 ???
@@ -436,4 +441,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)
@@ -495,5 +502,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
@@ -517,6 +524,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
+
 
 ! Martin
@@ -563,16 +573,18 @@
     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
     REAL, DIMENSION(klon)              :: ysolsw, ysollw
     REAL, DIMENSION(klon)              :: yfder
     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
     REAL, DIMENSION(klon)              :: y_dflux_t, y_dflux_q
     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
@@ -601,4 +613,5 @@
     REAL, DIMENSION(klon)              :: AcoefH, AcoefQ, BcoefH, BcoefQ
     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
@@ -609,12 +622,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
     REAL, DIMENSION(klon,klev)         :: ypplay, ydelp
     REAL, DIMENSION(klon,klev)         :: delp
@@ -678,4 +691,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
@@ -683,4 +697,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
@@ -862,4 +877,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)
@@ -980,10 +996,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.
@@ -1005,11 +1021,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.
@@ -1048,12 +1066,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 
@@ -1070,4 +1089,5 @@
     ycldt = 0.0      ; yrmu0 = 0.0
     ! Martin
+    y_d_qbs(:,:)=0.0
 
 !!! nrlmd+jyg le 02/05/2011 et le 20/02/2012
@@ -1395,4 +1415,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)
@@ -1406,4 +1427,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<
@@ -1419,5 +1441,6 @@
 !!! nrlmd le 13/06/2011
           y_delta_tsurf(j)=delta_tsurf(i,nsrf)
-!!!
+          yfluxbs(j)=0.0
+          y_flux_bs(j) = 0.0
        ENDDO
 ! >> PC
@@ -1458,4 +1481,5 @@
              yt(j,k) = t(i,k)
              yq(j,k) = q(i,k)
+             yqbs(j,k)=qbs(i,k)
           ENDDO
         ENDDO
@@ -1880,4 +1904,18 @@
        ENDIF  ! (iflag_split .eq.0)
 !!!
+
+ ! For blowing snow:
+    IF (ok_bs) THEN
+     ! following Bintanja et al 2000, part II
+     ! we assume that the eddy diffuvisity coefficient for
+     ! suspended particles is larger than Km by a factor zeta_bs
+     ! which is equal to 3 by default
+     ycoefqbs=ycoefm*zeta_bs 
+     CALL climb_qbs_down(knon, ycoefqbs, ypaprs, ypplay, &
+     ydelp, yt, yqbs, dtime, & 
+     CcoefQBS, DcoefQBS, &
+     Kcoef_qbs, gama_qbs, &
+     AcoefQBS, BcoefQBS)
+    ENDIF
 
 !****************************************************************************************
@@ -2057,5 +2095,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, & 
@@ -2063,5 +2101,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, &
@@ -2106,12 +2144,11 @@
                   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, &
                   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, &
@@ -2147,5 +2184,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, &
@@ -2277,4 +2314,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
 !
 ! ------------------------------------------------------------------------------
@@ -2583,4 +2627,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
+
 !!!
 !!
@@ -2696,12 +2750,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)
@@ -2873,4 +2935,14 @@
           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
 
 !      print*,'Dans pbl OK4'
@@ -3327,4 +3399,20 @@
        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
+
 !!!
 
Index: /LMDZ6/branches/blowing_snow/libf/phylmd/phyaqua_mod.F90
===================================================================
--- /LMDZ6/branches/blowing_snow/libf/phylmd/phyaqua_mod.F90	(revision 4484)
+++ /LMDZ6/branches/blowing_snow/libf/phylmd/phyaqua_mod.F90	(revision 4485)
@@ -360,4 +360,5 @@
     ql_ancien   = 0.
     qs_ancien   = 0.
+    qbs_ancien  = 0.
     u_ancien    = 0.
     v_ancien    = 0.
@@ -365,4 +366,5 @@
     prlw_ancien = 0.
     prsw_ancien = 0.  
+    prbsw_ancien= 0.
 
     ale_wake    = 0.
Index: /LMDZ6/branches/blowing_snow/libf/phylmd/phyredem.F90
===================================================================
--- /LMDZ6/branches/blowing_snow/libf/phylmd/phyredem.F90	(revision 4484)
+++ /LMDZ6/branches/blowing_snow/libf/phylmd/phyredem.F90	(revision 4485)
@@ -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, rneb_ancien, u_ancien, &
+                                prw_ancien, prlw_ancien, prsw_ancien, prbsw_ancien,      &
+                                ql_ancien, qs_ancien, qbs_ancien, rneb_ancien, u_ancien, &
                                 v_ancien, clwcon, rnebcon, ratqs, pbl_tke,   &
                                 wake_delta_pbl_tke, zmax0, f0, sig1, w01,    &
@@ -243,4 +243,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,"RNEBANCIEN", "RNEBANCIEN", rneb_ancien)
Index: /LMDZ6/branches/blowing_snow/libf/phylmd/phys_local_var_mod.F90
===================================================================
--- /LMDZ6/branches/blowing_snow/libf/phylmd/phys_local_var_mod.F90	(revision 4484)
+++ /LMDZ6/branches/blowing_snow/libf/phylmd/phys_local_var_mod.F90	(revision 4485)
@@ -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(:,:)
@@ -307,6 +311,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)
@@ -325,4 +329,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)
@@ -446,6 +452,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)
@@ -603,5 +609,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))
@@ -610,6 +616,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
@@ -634,5 +640,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))
@@ -760,6 +767,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))
@@ -776,4 +783,6 @@
       ALLOCATE(zxtsol(klon), snow_lsc(klon), zxfqfonte(klon), zxqsurf(klon))
       ALLOCATE(zxrunofflic(klon))
+      ALLOCATE(zxustartlic(klon), zxrhoslic(klon))
+      zxustartlic(:)=0. ; zxrhoslic(:)=0.
       ALLOCATE(rain_lsc(klon))
       ALLOCATE(rain_num(klon))
@@ -852,4 +861,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))
 
@@ -933,5 +943,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)
@@ -939,6 +949,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
@@ -963,4 +973,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)
       DEALLOCATE(d_u_vdf,d_v_vdf)
       DEALLOCATE(d_t_oli,d_t_oro)
@@ -1072,6 +1084,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)
 !
@@ -1086,4 +1098,5 @@
       DEALLOCATE(zxfqcalving, zxfluxlat)
       DEALLOCATE(zxrunofflic)
+      DEALLOCATE(zxustartlic, zxrhoslic)
       DEALLOCATE(zxtsol, snow_lsc, zxfqfonte, zxqsurf)
       DEALLOCATE(rain_lsc)
@@ -1154,5 +1167,5 @@
 
 
-      DEALLOCATE(prfl, psfl, fraca, Vprecip)
+      DEALLOCATE(prfl, psfl, bsfl, fraca, Vprecip)
       DEALLOCATE(zw2)
 
Index: /LMDZ6/branches/blowing_snow/libf/phylmd/phys_output_ctrlout_mod.F90
===================================================================
--- /LMDZ6/branches/blowing_snow/libf/phylmd/phys_output_ctrlout_mod.F90	(revision 4484)
+++ /LMDZ6/branches/blowing_snow/libf/phylmd/phys_output_ctrlout_mod.F90	(revision 4485)
@@ -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) /))
@@ -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) /))
@@ -1494,4 +1505,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 +1517,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) /))
@@ -1599,4 +1616,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 +1694,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
Index: /LMDZ6/branches/blowing_snow/libf/phylmd/phys_output_var_mod.F90
===================================================================
--- /LMDZ6/branches/blowing_snow/libf/phylmd/phys_output_var_mod.F90	(revision 4484)
+++ /LMDZ6/branches/blowing_snow/libf/phylmd/phys_output_var_mod.F90	(revision 4485)
@@ -30,5 +30,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)
@@ -36,8 +37,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/blowing_snow/libf/phylmd/phys_output_write_mod.F90
===================================================================
--- /LMDZ6/branches/blowing_snow/libf/phylmd/phys_output_write_mod.F90	(revision 4484)
+++ /LMDZ6/branches/blowing_snow/libf/phylmd/phys_output_write_mod.F90	(revision 4485)
@@ -17,5 +17,5 @@
   SUBROUTINE 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, flag_aerosol, flag_aerosol_strat, ok_cdnc)
@@ -47,5 +47,6 @@
          o_psol, o_mass, o_qsurf, o_qsol, &
          o_precip, o_rain_fall, o_rain_con, o_ndayrain, o_plul, o_pluc, o_plun, &
-         o_snow, o_msnow, o_fsnow, o_evap, o_ep,o_epmax_diag, & ! epmax_cape
+         o_snow, o_msnow, o_fsnow, o_evap, o_snowerosion, o_ustart_lic, o_rhosnow_lic, o_bsfall, & 
+         o_ep,o_epmax_diag, & ! epmax_cape
          o_tops, o_tops0, o_topl, o_topl0, &
          o_SWupTOA, o_SWupTOAclr, o_SWupTOAcleanclr, o_SWdnTOA, o_fdiffSWdnSFC, &
@@ -81,5 +82,5 @@
          o_Mipsh, o_Ma, &
          o_dnwd, o_dnwd0, o_ftime_deepcv, o_ftime_con, o_mc, &
-         o_prw, o_prlw, o_prsw, o_s_pblh, o_s_pblt, o_s_lcl, &
+         o_prw, o_prlw, o_prsw, o_prbsw, o_s_pblh, o_s_pblt, o_s_lcl, &
          o_s_therm, o_uSTDlevs, o_vSTDlevs, &
          o_wSTDlevs, o_zSTDlevs, o_qSTDlevs, &
@@ -106,5 +107,5 @@
          o_weakinv, o_dthmin, o_cldtau, &
          o_cldemi, o_pr_con_l, o_pr_con_i, &
-         o_pr_lsc_l, o_pr_lsc_i, o_re, o_fl, &
+         o_pr_lsc_l, o_pr_lsc_i, o_pr_bs, o_re, o_fl, &
          o_rh2m, &
          o_qsat2m, o_tpot, o_tpote, o_SWnetOR, &
@@ -134,14 +135,14 @@
          o_cldicemxrat, o_cldwatmxrat, o_reffclwtop, o_ec550aer, &
          o_lwcon, o_iwcon, o_temp, o_theta, &
-         o_ovapinit, o_ovap, o_oliq, o_ocond, o_geop, &
+         o_ovapinit, o_ovap, o_oliq, o_ocond, o_geop,o_qbs, &
          o_vitu, o_vitv, o_vitw, o_pres, o_paprs, &
          o_zfull, o_zhalf, o_rneb, o_rnebjn, o_rnebcon, &
          o_rnebls, o_rneblsvol, o_rhum, o_rhl, o_rhi, o_ozone, o_ozone_light, &
          o_duphy, o_dtphy, o_dqphy, o_dqphy2d, o_dqlphy, o_dqlphy2d, &
-         o_dqsphy, o_dqsphy2d, o_albe_srf, o_z0m_srf, o_z0h_srf, &
+         o_dqsphy, o_dqsphy2d, o_dqbsphy, o_dqbsphy2d, o_albe_srf, o_z0m_srf, o_z0h_srf, &
          o_ages_srf, o_snow_srf, o_alb1, o_alb2, o_tke, o_tke_dissip, &
          o_tke_max, o_kz, o_kz_max, o_clwcon, &
          o_dtdyn, o_dqdyn, o_dqdyn2d, o_dqldyn, o_dqldyn2d, & 
-         o_dqsdyn, o_dqsdyn2d, o_dudyn, o_dvdyn, &
+         o_dqsdyn, o_dqsdyn2d, o_dqbsdyn, o_dqbsdyn2d, o_dudyn, o_dvdyn, &
          o_dtcon, o_tntc, o_ducon, o_dvcon, &
          o_dqcon, o_dqcon2d, o_tnhusc, o_tnhusc, o_dtlsc, &
@@ -151,4 +152,5 @@
          o_ptconvth, o_lmaxth, o_dtvdf, &
          o_dtdis, o_dqvdf, o_dqvdf2d, o_dteva, o_dqeva, o_dqeva2d, &
+         o_dqbsvdf, o_dtbs, o_dqbs, o_dqbsbs, &
          o_ptconv, o_ratqs, o_dtthe, & 
          o_duthe, o_dvthe, o_ftime_th, &
@@ -241,5 +243,5 @@
     USE phys_state_var_mod, ONLY: heat_volc, cool_volc !NL
 
-    USE phys_state_var_mod, ONLY: pctsrf, rain_fall, snow_fall, &
+    USE phys_state_var_mod, ONLY: pctsrf, rain_fall, snow_fall, bs_fall,&
          qsol, z0m, z0h, fevap, agesno, &
          nday_rain, rain_con, snow_con, &
@@ -275,4 +277,5 @@
     USE phys_local_var_mod, ONLY: zxfluxlat, slp, ptstar, pt0, zxtsol, zt2m, &
          zn2mout, t2m_min_mon, t2m_max_mon, evap, &
+         snowerosion, zxustartlic, zxrhoslic, &
          l_mixmin,l_mix, tke_dissip, &
          zu10m, zv10m, zq2m, zustar, zxqsurf, &
@@ -288,5 +291,5 @@
          uwat, vwat, &
          rneb_seri, d_rneb_dyn, &
-         plcl, plfc, wbeff, convoccur, upwd, dnwd, dnwd0, prw, prlw, prsw, &
+         plcl, plfc, wbeff, convoccur, upwd, dnwd, dnwd0, prw, prlw, prsw,prbsw, &
          s_pblh, s_pblt, s_lcl, s_therm, uwriteSTD, &
          vwriteSTD, wwriteSTD, phiwriteSTD, qwriteSTD, &
@@ -308,5 +311,5 @@
          alp_bl_stat, alp_bl_fluct_tke, slab_wfbils, &
          weak_inversion, dthmin, cldtau, cldemi, &
-         pmflxr, pmflxs, prfl, psfl, re, fl, rh2m, &
+         pmflxr, pmflxs, prfl, psfl,bsfl, re, fl, rh2m, &
          qsat2m, tpote, tpot, d_ts, od443aer, od550aer, dryod550aer, &
          od865aer, abs550aer, od550lt1aer, sconcso4, sconcno3, &
@@ -324,9 +327,9 @@
          icc3dcon, icc3dstra, zfice, reffclwtop, &
          ec550aer, flwc, fiwc, t_seri, theta, q_seri, &
-         ql_seri, qs_seri, tr_seri, &
+         ql_seri, qs_seri, qbs_seri, tr_seri, qbs_seri,&
          zphi, u_seri, v_seri, omega, cldfra, &
          rneb, rnebjn, rneblsvol, zx_rh, zx_rhl, zx_rhi, d_t_dyn,  & 
-         d_q_dyn,  d_ql_dyn, d_qs_dyn, &
-         d_q_dyn2d,  d_ql_dyn2d, d_qs_dyn2d, &
+         d_q_dyn,  d_ql_dyn, d_qs_dyn, d_qbs_dyn,  &
+         d_q_dyn2d,  d_ql_dyn2d, d_qs_dyn2d, d_qbs_dyn2d, &
          d_u_dyn, d_v_dyn, d_t_con, d_t_ajsb, d_t_ajs, &
          d_u_ajs, d_v_ajs, &
@@ -336,4 +339,5 @@
          d_t_lscst, d_q_lscth, d_q_lscst, plul_th, &
          plul_st, d_t_vdf, d_t_diss, d_q_vdf, d_q_eva, &
+         d_t_bs, d_q_bs, d_qbs_bs, d_qbs_vdf, &
          zw2, fraca, zmax_th, d_q_ajsb, d_t_ec, d_u_vdf, &
          d_v_vdf, d_u_oro, d_v_oro, d_t_oro, d_u_lif, &
@@ -431,5 +435,5 @@
 
     ! Input
-    INTEGER :: itap, ivap, iliq, isol, read_climoz
+    INTEGER :: itap, ivap, iliq, isol, ibs, read_climoz
     INTEGER, DIMENSION(klon) :: lmax_th
     LOGICAL :: aerosol_couple, ok_sync
@@ -887,7 +891,13 @@
 
        IF (vars_defined) THEN
-          DO i = 1, klon
+          IF (ok_bs) THEN
+             DO i = 1, klon
+             zx_tmp_fi2d(i) = rain_fall(i) + snow_fall(i) + bs_fall(i)
+             ENDDO
+          ELSE
+             DO i = 1, klon
              zx_tmp_fi2d(i) = rain_fall(i) + snow_fall(i)
-          ENDDO
+             ENDDO
+          ENDIF
        ENDIF
 
@@ -919,4 +929,11 @@
        CALL histwrite_phy(o_fsnow, zfra_o)
        CALL histwrite_phy(o_evap, evap)
+
+       IF (ok_bs) THEN
+           CALL histwrite_phy(o_bsfall, bs_fall)     
+           CALL histwrite_phy(o_snowerosion, snowerosion)
+           CALL histwrite_phy(o_ustart_lic, zxustartlic)
+           CALL histwrite_phy(o_rhosnow_lic, zxrhoslic)
+       ENDIF
 
        IF (vars_defined) THEN
@@ -1105,15 +1122,4 @@
        ENDIF
        CALL histwrite_phy(o_tauy, zx_tmp_fi2d)
-
-       ! Etienne: test sorties pour compil sur JZ
-!       IF (landice_opt .GE. 1) THEN
-!          CALL histwrite_phy(o_snowsrf, snow_o)
-!          CALL histwrite_phy(o_qsnow, qsnow)
-!          CALL histwrite_phy(o_snowhgt,snowhgt)
-!          CALL histwrite_phy(o_toice,to_ice)
-!          CALL histwrite_phy(o_sissnow,sissnow)
-!          CALL histwrite_phy(o_runoff,runoff)
-!          CALL histwrite_phy(o_albslw3,albsol3_lic)
-!       ENDIF
 
        DO nsrf = 1, nbsrf
@@ -1281,4 +1287,7 @@
        CALL histwrite_phy(o_prlw, prlw)
        CALL histwrite_phy(o_prsw, prsw)
+       IF (ok_bs) THEN
+       CALL histwrite_phy(o_prbsw, prbsw)
+       ENDIF
        CALL histwrite_phy(o_s_pblh, s_pblh)
        CALL histwrite_phy(o_s_pblt, s_pblt)
@@ -1526,4 +1535,9 @@
        CALL histwrite_phy(o_re, re)
        CALL histwrite_phy(o_fl, fl)
+
+       IF (ok_bs) THEN
+         CALL histwrite_phy(o_pr_bs, bsfl(:,1:klev))
+       ENDIF
+
        IF (vars_defined) THEN
           DO i=1, klon
@@ -1795,4 +1809,8 @@
        CALL histwrite_phy(o_zfull,zx_tmp_fi3d)
 
+       IF (ok_bs) THEN
+          CALL histwrite_phy(o_qbs, qbs_seri)
+       ENDIF
+
 #ifdef CPP_XIOS
 !solbnd begin
@@ -1905,4 +1923,16 @@
        zx_tmp_fi2d=0.0
        CALL histwrite_phy(o_dqsphy2d,  zx_tmp_fi2d)
+       ENDIF
+
+
+       IF (ok_bs) THEN
+       CALL histwrite_phy(o_dqbsphy,  d_qx(:,:,ibs))
+       IF (vars_defined) CALL water_int(klon,klev,d_qx(:,:,ibs),zmasse,zx_tmp_fi2d)
+       CALL histwrite_phy(o_dqbsphy2d,  zx_tmp_fi2d)
+       ELSE
+       zx_tmp_fi3d=0.0
+       CALL histwrite_phy(o_dqbsphy,  zx_tmp_fi3d)
+       zx_tmp_fi2d=0.0
+       CALL histwrite_phy(o_dqbsphy2d,  zx_tmp_fi2d)
        ENDIF
 
@@ -1955,4 +1985,9 @@
 
        CALL histwrite_phy(o_dqsdyn2d, d_qs_dyn2d)
+
+       IF (ok_bs) THEN
+       CALL histwrite_phy(o_dqbsdyn, d_qbs_dyn)
+       CALL histwrite_phy(o_dqbsdyn2d, d_qbs_dyn2d)
+       ENDIF
 
        CALL histwrite_phy(o_dudyn, d_u_dyn)
@@ -2098,4 +2133,15 @@
        ENDIF
        CALL histwrite_phy(o_dvthe, zx_tmp_fi3d)
+
+       IF (ok_bs) THEN
+          IF (vars_defined) zx_tmp_fi3d(1:klon,1:klev)=d_qbs_vdf(1:klon,1:klev)/pdtphys
+          CALL histwrite_phy(o_dqbsvdf, zx_tmp_fi3d)
+          IF (vars_defined) zx_tmp_fi3d(1:klon,1:klev)=d_qbs_bs(1:klon,1:klev)/pdtphys
+          CALL histwrite_phy(o_dqbsbs, zx_tmp_fi3d)
+          IF (vars_defined) zx_tmp_fi3d(1:klon,1:klev)=d_q_bs(1:klon,1:klev)/pdtphys
+          CALL histwrite_phy(o_dqbs, zx_tmp_fi3d)
+          IF (vars_defined) zx_tmp_fi3d(1:klon,1:klev)=d_t_bs(1:klon,1:klev)/pdtphys
+          CALL histwrite_phy(o_dtbs, zx_tmp_fi3d)
+       ENDIF
 
        IF (iflag_thermals>=1) THEN
Index: /LMDZ6/branches/blowing_snow/libf/phylmd/phys_state_var_mod.F90
===================================================================
--- /LMDZ6/branches/blowing_snow/libf/phylmd/phys_state_var_mod.F90	(revision 4484)
+++ /LMDZ6/branches/blowing_snow/libf/phylmd/phys_state_var_mod.F90	(revision 4485)
@@ -54,6 +54,6 @@
 
 
-      REAL, ALLOCATABLE, SAVE :: rain_fall(:), snow_fall(:)
-!$OMP THREADPRIVATE( rain_fall, snow_fall)
+      REAL, ALLOCATABLE, SAVE :: rain_fall(:), snow_fall(:), bs_fall(:)
+!$OMP THREADPRIVATE( rain_fall, snow_fall, bs_fall)
       REAL, ALLOCATABLE, SAVE :: solsw(:), solswfdiff(:), sollw(:)
 !$OMP THREADPRIVATE(solsw, solswfdiff, sollw)
@@ -82,8 +82,8 @@
       REAL, ALLOCATABLE, SAVE :: t_ancien(:,:), q_ancien(:,:)
 !$OMP THREADPRIVATE(t_ancien, q_ancien)
-      REAL, ALLOCATABLE, SAVE :: ql_ancien(:,:), qs_ancien(:,:)
-!$OMP THREADPRIVATE(ql_ancien, qs_ancien)
-      REAL, ALLOCATABLE, SAVE :: prw_ancien(:), prlw_ancien(:), prsw_ancien(:)
-!$OMP THREADPRIVATE(prw_ancien, prlw_ancien, prsw_ancien)
+      REAL, ALLOCATABLE, SAVE :: ql_ancien(:,:), qs_ancien(:,:), qbs_ancien(:,:)
+!$OMP THREADPRIVATE(ql_ancien, qs_ancien, qbs_ancien)
+      REAL, ALLOCATABLE, SAVE :: prw_ancien(:), prlw_ancien(:), prsw_ancien(:), prbsw_ancien(:)
+!$OMP THREADPRIVATE(prw_ancien, prlw_ancien, prsw_ancien, prbsw_ancien)
 #ifdef ISO
       REAL, ALLOCATABLE, SAVE :: xt_ancien(:,:,:),xtl_ancien(:,:,:),xts_ancien(:,:,:)
@@ -547,4 +547,5 @@
       ALLOCATE(rain_fall(klon))
       ALLOCATE(snow_fall(klon))
+      ALLOCATE(bs_fall(klon))
       ALLOCATE(solsw(klon), solswfdiff(klon), sollw(klon))
       sollw=0.0
@@ -556,6 +557,6 @@
       ALLOCATE(rugoro(klon))
       ALLOCATE(t_ancien(klon,klev), q_ancien(klon,klev))
-      ALLOCATE(ql_ancien(klon,klev), qs_ancien(klon,klev))
-      ALLOCATE(prw_ancien(klon), prlw_ancien(klon), prsw_ancien(klon))
+      ALLOCATE(ql_ancien(klon,klev), qs_ancien(klon,klev), qbs_ancien(klon,klev))
+      ALLOCATE(prw_ancien(klon), prlw_ancien(klon), prsw_ancien(klon), prbsw_ancien(klon))
       ALLOCATE(u_ancien(klon,klev), v_ancien(klon,klev))
 !!! Rom P >>>
@@ -763,10 +764,10 @@
 !FC
       DEALLOCATE(treedrg)
-      DEALLOCATE(rain_fall, snow_fall, solsw, solswfdiff, sollw, radsol, swradcorr)
+      DEALLOCATE(rain_fall, snow_fall, bs_fall,solsw, solswfdiff, sollw, radsol, swradcorr)
       DEALLOCATE(zmea, zstd, zsig, zgam)
       DEALLOCATE(zthe, zpic, zval)
       DEALLOCATE(rugoro, t_ancien, q_ancien, clwcon, rnebcon)
-      DEALLOCATE(qs_ancien, ql_ancien, rneb_ancien)
-      DEALLOCATE(prw_ancien, prlw_ancien, prsw_ancien)
+      DEALLOCATE(qs_ancien, ql_ancien, qbs_ancien, rneb_ancien)
+      DEALLOCATE(prw_ancien, prlw_ancien, prsw_ancien, prbsw_ancien)
       DEALLOCATE(qtc_cv,sigt_cv)
       DEALLOCATE(u_ancien, v_ancien)
Index: /LMDZ6/branches/blowing_snow/libf/phylmd/physiq_mod.F90
===================================================================
--- /LMDZ6/branches/blowing_snow/libf/phylmd/physiq_mod.F90	(revision 4484)
+++ /LMDZ6/branches/blowing_snow/libf/phylmd/physiq_mod.F90	(revision 4485)
@@ -84,4 +84,5 @@
     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
 
@@ -144,8 +145,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, &
@@ -164,8 +165,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, &
@@ -215,6 +218,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,                   &
@@ -304,5 +307,5 @@
        fsolsw, wfbils, wfbilo,  &
        wfevap, wfrain, wfsnow,  &  
-       prfl, psfl, fraca, Vprecip,  &
+       prfl, psfl,bsfl, fraca, Vprecip,  &
        zw2,  &
        !
@@ -460,7 +463,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)
     !
     !
@@ -833,7 +836,9 @@
     !XXX PB 
     REAL fluxq(klon,klev, nbsrf)   ! flux turbulent d'humidite
+    REAL fluxqbs(klon,klev, nbsrf)   ! flux turbulent de neige soufflee
     !
     REAL zxfluxt(klon, klev)
     REAL zxfluxq(klon, klev)
+    REAL zxfluxqbs(klon,klev)
     REAL zxfluxu(klon, klev)
     REAL zxfluxv(klon, klev)
@@ -919,5 +924,5 @@
     !
     ! tendance nulles
-    REAL, dimension(klon,klev):: du0, dv0, dt0, dq0, dql0, dqi0
+    REAL, dimension(klon,klev):: du0, dv0, dt0, dq0, dql0, dqi0, dqbs0
     REAL, dimension(klon)     :: dsig0, ddens0
     INTEGER, dimension(klon)  :: wkoccur1
@@ -1060,5 +1065,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
     !IM: t2m, q2m, ustar, u10m, v10m et t2mincels, t2maxcels
     CHARACTER*40 t2mincels, t2maxcels       !t2m min., t2m max
@@ -1270,4 +1276,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)
@@ -1319,5 +1326,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.'
@@ -1337,4 +1344,13 @@
           CALL abort_physic(modname,abort_message,1)
        ENDIF
+
+        IF (ok_bs) THEN
+         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
@@ -1754,4 +1770,8 @@
            CALL lscp_ini(pdtphys,ok_ice_sursat)
        endif
+       CALL blowing_snow_ini(prt_level,lunout, &
+                             RCPD, RLSTT, RLVTT, RLMLT, &
+                             RVTMP2, RTT,RD,RG)
+
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
@@ -1779,5 +1799,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,           &
@@ -2250,4 +2270,5 @@
     dql0(:,:)=0.
     dqi0(:,:)=0.
+    dqbs0(:,:)=0.
     dsig0(:) = 0.
     ddens0(:) = 0.
@@ -2302,4 +2323,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
@@ -2309,8 +2331,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
@@ -2321,5 +2351,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
     ENDIF
@@ -2377,4 +2407,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
@@ -2383,4 +2414,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
@@ -2398,8 +2431,10 @@
        d_ql_dyn2d(:) = 0.0
        d_qs_dyn2d(:) = 0.0
+       d_qbs_dyn2d(:)= 0.0
        ! !! RomP >>>   td dyn traceur
        IF (nqtot > nqo) d_tr_dyn(:,:,:)= 0.0
        ! !! RomP <<<
        d_rneb_dyn(:,:)=0.0
+       d_qbs_dyn(:,:)=0.0
        ancien_ok = .TRUE.
     ENDIF
@@ -2517,5 +2552,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)
     CALL prt_enerbil('eva',itap)
@@ -2682,7 +2717,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, &
@@ -2695,9 +2730,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<
@@ -2725,5 +2760,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, &
@@ -2750,12 +2785,13 @@
        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)
        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)
        ENDIF
        CALL prt_enerbil('vdf',itap)
+
        !--------------------------------------------------------------------
 
@@ -2808,4 +2844,24 @@
 
     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)
+
+    ENDIF
+
     ! =================================================================== c
     !   Calcul de Qsat
@@ -3232,5 +3288,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)
     CALL prt_enerbil('convection',itap)
@@ -3363,5 +3419,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)
        CALL prt_enerbil('wake',itap)
@@ -3521,5 +3577,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)
           CALL prt_enerbil('thermals',itap)
           !
@@ -3583,5 +3639,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)
           CALL prt_enerbil('ajsb',itap)
@@ -3672,5 +3728,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)
     CALL prt_enerbil('lsc',itap)
@@ -3695,12 +3751,40 @@
     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?
+          !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
+
     IF (check) THEN
        za = qcheck(klon,klev,paprs,q_seri,ql_seri,cell_area)
@@ -4489,7 +4573,7 @@
     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)
     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)
     CALL prt_enerbil('LW',itap)
 
@@ -4566,5 +4650,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)
        CALL prt_enerbil('oro',itap)
@@ -4617,5 +4701,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)
        CALL prt_enerbil('lif',itap)
@@ -4642,5 +4726,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)
        CALL prt_enerbil('hin',itap)
     ENDIF
@@ -4660,5 +4744,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)
        CALL prt_enerbil('front_gwd_rando',itap)
@@ -4669,5 +4753,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)
        CALL prt_enerbil('flott_gwd_rando',itap)
@@ -4721,5 +4805,5 @@
        ! ajout de la tendance d'humidite due au methane
        d_q_ch4_dtime(:,:) = d_q_ch4(:,:)*phys_tstep
-       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)
        d_q_ch4(:,:) = d_q_ch4_dtime(:,:)/phys_tstep
@@ -5026,5 +5110,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'
@@ -5103,11 +5187,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
     !
@@ -5168,7 +5255,12 @@
           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
@@ -5217,8 +5309,10 @@
     ql_ancien(:,:) = ql_seri(:,:)
     qs_ancien(:,:) = qs_seri(:,:)
+    qbs_ancien(:,:) = qbs_seri(:,:)
     rneb_ancien(:,:) = rneb_seri(:,:)
     CALL water_int(klon,klev,q_ancien,zmasse,prw_ancien)
     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(:,:,:)
@@ -5345,5 +5439,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/blowing_snow/libf/phylmd/surf_land_mod.F90
===================================================================
--- /LMDZ6/branches/blowing_snow/libf/phylmd/surf_land_mod.F90	(revision 4484)
+++ /LMDZ6/branches/blowing_snow/libf/phylmd/surf_land_mod.F90	(revision 4485)
@@ -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 , & 
@@ -81,5 +81,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
@@ -109,5 +109,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
@@ -125,4 +125,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
 
@@ -130,4 +131,16 @@
     REAL, DIMENSION(klon)      :: alb1_new,alb2_new
 !albedo SB <<<
+
+
+!****************************************************************************************
+!Total solid precip
+
+IF (ok_bs) THEN
+precip_totsnow=precip_snow+precip_bs
+ELSE
+precip_totsnow=precip_snow
+ENDIF
+!****************************************************************************************
+
 
 !**************************************************************************************** 
@@ -165,7 +178,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,fluxbs, &              
             tsol_rad, tsurf_new, alb1_new, alb2_new, &
             emis_new, z0m, z0h, qsurf, &
@@ -184,5 +197,5 @@
 !****************************************************************************************
        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, &
@@ -192,5 +205,10 @@
         z0h(1:knon)=z0m(1:knon) ! En attendant mieux
 
+
     ENDIF ! ok_veget
+
+        ! blowing snow not treated yet over land
+        fluxbs(:)=0.
+
 
 !****************************************************************************************
Index: /LMDZ6/branches/blowing_snow/libf/phylmd/surf_land_orchidee_mod.F90
===================================================================
--- /LMDZ6/branches/blowing_snow/libf/phylmd/surf_land_orchidee_mod.F90	(revision 4484)
+++ /LMDZ6/branches/blowing_snow/libf/phylmd/surf_land_orchidee_mod.F90	(revision 4485)
@@ -47,5 +47,5 @@
        precip_rain, precip_snow, lwdown, swnet, swdown, &
        ps, q2m, t2m, &
-       evap, fluxsens, fluxlat, &              
+       evap, fluxsens, fluxlat, fluxbs, &              
        tsol_rad, tsurf_new, alb1_new, alb2_new, &
        emis_new, z0m_new, z0h_new, qsurf, &
@@ -140,5 +140,5 @@
 ! Parametres de sortie
 !****************************************************************************************
-    REAL, DIMENSION(klon), INTENT(OUT)        :: evap, fluxsens, fluxlat, qsurf
+    REAL, DIMENSION(klon), INTENT(OUT)        :: evap, fluxsens, fluxlat, fluxbs, qsurf
     REAL, DIMENSION(klon), INTENT(OUT)        :: tsol_rad, tsurf_new
     REAL, DIMENSION(klon), INTENT(OUT)        :: alb1_new, alb2_new
@@ -632,5 +632,8 @@
     fluxsens(1:knon) = -1. * fluxsens(1:knon)
     fluxlat(1:knon)  = -1. * fluxlat(1:knon)
-    
+   
+! temporary until blowing snow flux is calculated in orchidee
+    fluxbs(:)=0.
+
 !  evap     = -1. * evap
 
Index: /LMDZ6/branches/blowing_snow/libf/phylmd/surf_landice_mod.F90
===================================================================
--- /LMDZ6/branches/blowing_snow/libf/phylmd/surf_landice_mod.F90	(revision 4484)
+++ /LMDZ6/branches/blowing_snow/libf/phylmd/surf_landice_mod.F90	(revision 4485)
@@ -12,10 +12,10 @@
        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, &
        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, &
@@ -30,8 +30,9 @@
     USE cpl_mod,          ONLY : cpl_send_landice_fields
     USE calcul_fluxs_mod
+    USE phys_local_var_mod, ONLY : zxrhoslic, zxustartlic
     USE phys_output_var_mod, ONLY : snow_o,zfra_o
 !FC
     USE ioipsl_getin_p_mod, ONLY : getin_p
-
+    USE blowing_snow_ini_mod, ONLY : zeta_bs, pbst_bs, prt_bs
 
 #ifdef CPP_INLANDSIS
@@ -56,5 +57,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
@@ -62,5 +63,5 @@
     REAL, DIMENSION(klon), INTENT(IN)             :: AcoefU, AcoefV, BcoefU, BcoefV
     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
@@ -94,4 +95,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      
@@ -134,5 +136,9 @@
 
     REAL,DIMENSION(klon) :: alb1,alb2
+    REAL,DIMENSION(klon) :: precip_totsnow, evap_totsnow
     REAL, DIMENSION (klon,6) :: alb6
+    REAL                   :: rho0, rhoice, ustart0, hsalt, esalt, qsalt
+    REAL                   :: tau_dens, tau_dens0, tau_densmin, rhomax, rhohard
+    REAL, DIMENSION(klon)  :: ws1, rhos, ustart
 ! End definition
 !****************************************************************************************
@@ -165,5 +171,5 @@
     alb2(:) = 999999.
     alb1(:) = 999999.
-    
+    fluxbs(:)=0.  
     runoff(:) = 0.
 !****************************************************************************************
@@ -173,4 +179,6 @@
     radsol(:) = 0.0
     radsol(1:knon) = swnet(1:knon) + lwnet(1:knon)
+
+!****************************************************************************************
 
 !****************************************************************************************
@@ -264,5 +272,4 @@
 
 
-
     ELSE 
 
@@ -314,12 +321,4 @@
          flux_u1, flux_v1)
 
-!****************************************************************************************
-! Calculate snow height, age, run-off,..
-!    
-!****************************************************************************************
-    CALL fonte_neige(knon, is_lic, knindex, dtime, &
-         tsurf, precip_rain, precip_snow, &
-         snow, qsol, tsurf_new, evap)
-
 
 !****************************************************************************************
@@ -327,10 +326,12 @@
 !
 !****************************************************************************************
+ 
     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))
+
+! 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"
@@ -354,4 +355,63 @@
     !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
+       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)
+           esalt=1./(3.25*max(ustar(i),0.001))
+           hsalt=0.08436*ustar(i)**1.27
+           qsalt=(max(ustar(i)**2-ustart(i)**2,0.))/(RG*hsalt)*esalt
+           !ep=qsalt*cdragm(i)*sqrt(u1(i)**2+v1(i)**2)
+           fluxbs(i)= zeta_bs*p1lay(i)/RD/temp_air(i)*ws1(i)*cdragm(i)*(qbs1(i)-qsalt)
+       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_totsnow,  &
+         snow, qsol, tsurf_new, evap_totsnow)
+
+    WHERE (snow(1 : knon) .LT. 0.0001) agesno(1 : knon) = 0.                                         
+    zfra(1:knon) = MAX(0.0,MIN(1.0,snow(1:knon)/(snow(1:knon)+10.0)))  
 
 
@@ -378,6 +438,4 @@
     runoff(1:knon)=run_off_lic(1:knon)/dtime
 
-  
-!****************************************************************************************
        snow_o=0.
        zfra_o = 0.
@@ -420,6 +478,4 @@
 !albedo SB <<<
 
-  
-  
 
   END SUBROUTINE surf_landice
Index: /LMDZ6/branches/blowing_snow/libf/phylmd/surf_ocean_mod.F90
===================================================================
--- /LMDZ6/branches/blowing_snow/libf/phylmd/surf_ocean_mod.F90	(revision 4484)
+++ /LMDZ6/branches/blowing_snow/libf/phylmd/surf_ocean_mod.F90	(revision 4485)
@@ -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, &
@@ -60,5 +60,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
@@ -144,4 +144,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)
@@ -160,4 +161,15 @@
     radsol(1:klon) = 0.0 ! initialisation a priori inutile
     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
+
 
     !******************************************************************************
@@ -188,5 +200,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, &
@@ -200,5 +212,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, &
@@ -211,5 +223,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, &
@@ -326,5 +338,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), &
Index: /LMDZ6/branches/blowing_snow/libf/phylmdiso/blowing_snow_ini_mod.F90
===================================================================
--- /LMDZ6/branches/blowing_snow/libf/phylmdiso/blowing_snow_ini_mod.F90	(revision 4485)
+++ /LMDZ6/branches/blowing_snow/libf/phylmdiso/blowing_snow_ini_mod.F90	(revision 4485)
@@ -0,0 +1,1 @@
+link ../phylmd/blowing_snow_ini_mod.F90
Index: /LMDZ6/branches/blowing_snow/libf/phylmdiso/blowing_snow_sublim_sedim.F90
===================================================================
--- /LMDZ6/branches/blowing_snow/libf/phylmdiso/blowing_snow_sublim_sedim.F90	(revision 4485)
+++ /LMDZ6/branches/blowing_snow/libf/phylmdiso/blowing_snow_sublim_sedim.F90	(revision 4485)
@@ -0,0 +1,1 @@
+link ../phylmd/blowing_snow_sublim_sedim.F90
Index: /LMDZ6/branches/blowing_snow/libf/phylmdiso/call_blowing_snow_sublim_sedim.F90
===================================================================
--- /LMDZ6/branches/blowing_snow/libf/phylmdiso/call_blowing_snow_sublim_sedim.F90	(revision 4485)
+++ /LMDZ6/branches/blowing_snow/libf/phylmdiso/call_blowing_snow_sublim_sedim.F90	(revision 4485)
@@ -0,0 +1,1 @@
+link ../phylmd/call_blowing_snow_sublim_sedim.F90
Index: /LMDZ6/branches/blowing_snow/libf/phylmdiso/climb_qbs_mod.F90
===================================================================
--- /LMDZ6/branches/blowing_snow/libf/phylmdiso/climb_qbs_mod.F90	(revision 4485)
+++ /LMDZ6/branches/blowing_snow/libf/phylmdiso/climb_qbs_mod.F90	(revision 4485)
@@ -0,0 +1,1 @@
+link ../phylmd/climb_qbs_mod.F90
