Index: LMDZ6/trunk/libf/phylmd/physiq_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/physiq_mod.F90	(revision 4519)
+++ LMDZ6/trunk/libf/phylmd/physiq_mod.F90	(revision 4523)
@@ -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,9 +145,9 @@
        ! [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, &
        rhcl, &        
        ! 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, &
@@ -165,8 +166,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, &
@@ -216,6 +219,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,                   &
@@ -305,5 +308,5 @@
        fsolsw, wfbils, wfbilo,  &
        wfevap, wfrain, wfsnow,  &  
-       prfl, psfl, fraca, Vprecip,  &
+       prfl, psfl,bsfl, fraca, Vprecip,  &
        zw2,  &
        !
@@ -461,7 +464,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)
     !
     !
@@ -834,7 +837,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)
@@ -920,5 +925,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
@@ -1066,5 +1071,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
@@ -1276,4 +1282,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)
@@ -1325,5 +1332,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.'
@@ -1343,4 +1350,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
@@ -1761,4 +1777,8 @@
            CALL lscp_ini(pdtphys,ok_ice_sursat)
        endif
+       CALL blowing_snow_ini(prt_level,lunout, &
+                             RCPD, RLSTT, RLVTT, RLMLT, &
+                             RVTMP2, RTT,RD,RG)
+
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
@@ -1786,5 +1806,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,           &
@@ -2257,4 +2277,5 @@
     dql0(:,:)=0.
     dqi0(:,:)=0.
+    dqbs0(:,:)=0.
     dsig0(:) = 0.
     ddens0(:) = 0.
@@ -2309,4 +2330,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
@@ -2316,8 +2338,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
@@ -2328,5 +2358,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
@@ -2384,4 +2414,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
@@ -2390,4 +2421,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
@@ -2405,8 +2438,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
@@ -2524,5 +2559,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)
@@ -2689,7 +2724,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, &
@@ -2702,9 +2737,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<
@@ -2732,5 +2767,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, &
@@ -2757,12 +2792,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)
+
        !--------------------------------------------------------------------
 
@@ -2815,4 +2851,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
@@ -3239,5 +3295,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)
@@ -3370,5 +3426,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)
@@ -3528,5 +3584,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)
           !
@@ -3590,5 +3646,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)
@@ -3679,5 +3735,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)
@@ -3702,12 +3758,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)
@@ -4510,7 +4594,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)
 
@@ -4587,5 +4671,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)
@@ -4638,5 +4722,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)
@@ -4663,5 +4747,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
@@ -4681,5 +4765,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)
@@ -4690,5 +4774,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)
@@ -4742,5 +4826,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
@@ -5047,5 +5131,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'
@@ -5124,11 +5208,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
     !
@@ -5189,7 +5276,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
@@ -5238,8 +5330,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(:,:,:)
@@ -5366,5 +5460,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,           &
