Index: /LMDZ6/trunk/libf/phylmd/lmdz_blowing_snow_ini.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/lmdz_blowing_snow_ini.F90	(revision 4820)
+++ /LMDZ6/trunk/libf/phylmd/lmdz_blowing_snow_ini.F90	(revision 4821)
@@ -3,12 +3,10 @@
 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 :: iflag_saltation_bs
+   real, save, protected :: RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT,RD,RG
+   real, save, protected :: coef_eva_bs, expo_eva_bs, fallv_bs, zeta_bs
+   real, save, protected :: prt_bs, pbst_bs, qbst_bs
+   
+   integer, save, protected :: iflag_saltation_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)
@@ -16,28 +14,22 @@
    !$OMP THREADPRIVATE(iflag_saltation_bs)
 
-    real tbsmelt                  ! parameter to calculate melting fraction of BS sedimentation
-    parameter (tbsmelt=278.15)
+   real, save, protected :: tbsmelt=278.15     ! parameter to calculate melting fraction of BS sedimentation
+   real, save, protected :: taumeltbs0=1800.0  ! Melting time scale of blowing snow at 273.15K
+   real, save, protected :: qbsmin=1.E-10      ! Minimum blowing snow specific content
 
-    real taumeltbs0               ! Melting time scale of blowing snow at 273.15K
-    parameter (taumeltbs0=1800.0)
-
-    real qbsmin                   ! Minimum blowing snow specific content
-    parameter (qbsmin=1.E-10)
+   !$OMP THREADPRIVATE(tbsmelt, taumeltbs0, qbsmin)
 
 
-      contains
+    contains
 
-      subroutine blowing_snow_ini(prt_level_in,lunout_in, &
-                                  RCPD_in, RLSTT_in, RLVTT_in, RLMLT_in,&
+   subroutine blowing_snow_ini(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
@@ -49,5 +41,4 @@
          RG=RG_in
          RVTMP2=RVTMP2_in
-         lunout=lunout_in
 
 
@@ -60,5 +51,4 @@
          prt_bs= 0.0003
          CALL getin_p('prt_bs',prt_bs)
-
 
          zeta_bs= 3.
Index: /LMDZ6/trunk/libf/phylmd/lmdz_blowing_snow_sublim_sedim.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/lmdz_blowing_snow_sublim_sedim.F90	(revision 4820)
+++ /LMDZ6/trunk/libf/phylmd/lmdz_blowing_snow_sublim_sedim.F90	(revision 4821)
@@ -51,10 +51,10 @@
 
 integer                                  :: k,i,n
-real                                     :: zqev0, zqevi, zqevti, zcpair, zcpeau, dqbsmelt
-real                                     :: dqsedim,precbs, deltaqchaud, zmelt, taumeltbs
-real                                     :: maxdeltaqchaud
+real                                     :: cpd, cpw, dqsub, maxdqsub, dqbsmelt
+real                                     :: dqsedim,precbs, dqmelt, zmelt, taumeltbs
+real                                     :: maxdqmelt, rhoair, dz
 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
+real, dimension(ngrid)                   :: zqbsi, zqbs_up, zmair
+real, dimension(ngrid)                   :: zvelo 
 
 !++++++++++++++++++++++++++++++++++++++++++++++++++
@@ -66,9 +66,10 @@
 dq_bs(:,:)=0.
 dqbs_bs(:,:)=0.
-velo(:,:)=0.
+zvelo(:)=0.
 zt(:)=0.
 zq(:)=0.
 zqbs(:)=0.
 sedim(:)=0.
+sedimn(:)=0.
 
 ! begin of top-down loop
@@ -81,27 +82,20 @@
     ENDDO
 
-
+     ! thermalization of blowing snow precip coming from above     
     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 
+            cpd=RCPD*(1.0+RVTMP2*zq(i))
+            cpw=RCPD*RVTMP2
+            ! zqbs_up: blowing snow 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)
+            zqbs_up(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)
+            zt(i) = ( (temp(i,k+1)+dtemp_bs(i,k+1))*zqbs_up(i)*cpw + cpd*zt(i) ) &
+             / (cpd + zqbs_up(i)*cpw)
         ENDDO
-    ELSE
-        DO i = 1, ngrid
-            zmair(i)=(paprs(i,k)-paprs(i,k+1))/RG
-            zmqc(i) = 0.
-        ENDDO
-
     ENDIF
 
@@ -111,107 +105,68 @@
     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)
     
+    
+    ! 3 processes: melting, sublimation and precipitation of blowing snow
     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)
-       ! BS fall velocity
-       velo(i,k) = fallv_bs
-
 
           IF (zt(i) .GT. RTT) THEN
              
-             ! if positive celcius temperature, we assume
-             ! that part of the the blowing snow flux melts and evaporates
-             
-             ! vapor, bs, temperature, precip fluxes update
-             zmelt = ((zt(i)-RTT)/(tbsmelt-RTT))
-             zmelt = MIN(MAX(zmelt,0.),1.)
-             sedimn(i)=sedim(i)*zmelt
-             deltaqchaud=-(sedimn(i)-sedim(i))*(RG/(paprs(i,k)-paprs(i,k+1)))*dtime
-             ! max evap such as celcius temperature cannot become negative
-             maxdeltaqchaud= max(RCPD*(1.0+RVTMP2*(zq(i)+zqbs(i)))*(zt(i)-RTT)/(RLMLT+RLVTT),0.)
-
-             deltaqchaud=min(max(deltaqchaud,0.),maxdeltaqchaud)
-             zq(i) = zq(i) + deltaqchaud  
-              
-             ! melting + evaporation
-             zt(i) = zt(i) - deltaqchaud * (RLMLT+RLVTT)/RCPD/(1.0+RVTMP2*(zq(i)+zqbs(i)))
-
-             sedim(i)=sedimn(i)
-
-             ! if temperature still positive, we assume that part of the blowing snow 
-             ! already present in the mesh melts and evaporates with a typical time 
-             ! constant between taumeltbs0 and 0
- 
-             IF ( zt(i) .GT. RTT ) THEN
-                taumeltbs=taumeltbs0*exp(-max(0.,(zt(i)-RTT)/(tbsmelt-RTT)))
-                deltaqchaud=min(zqbs(i),zqbs(i)/taumeltbs*dtime)
-                maxdeltaqchaud= max(RCPD*(1.0+RVTMP2*(zq(i)+zqbs(i)))*(zt(i)-RTT)/(RLMLT+RLVTT),0.)
-                deltaqchaud=min(max(deltaqchaud,0.),maxdeltaqchaud)
-                zq(i) = zq(i) + deltaqchaud
-
-                ! melting + evaporation
-                zt(i) = zt(i) - deltaqchaud * (RLMLT+RLVTT)/RCPD/(1.0+RVTMP2*(zq(i)+zqbs(i)))
-                ! qbs update
-                zqbs(i)=max(zqbs(i)-deltaqchaud,0.)
-             ENDIF
-
-             ! now sedimentation scheme with an exact numerical resolution 
-             ! (correct if fall velocity is constant)
-             zqbsi(i)=zqbs(i)
-             zqbs(i) = zqbsi(i)*exp(-velo(i,k)/zdz(i)*dtime)+sedim(i)/zrho(i,k)/velo(i,k)
-             zqbs(i) = max(zqbs(i),0.)
-
-             ! flux update
-             sedim(i) = sedim(i) + zrho(i,k)*zdz(i)/dtime*(zqbsi(i)-zqbs(i))
-             sedim(i) = max(0.,sedim(i))
-      
+             ! if temperature is positive, we assume that part of the blowing snow 
+             ! already present  melts and evaporates with a typical time 
+             ! constant taumeltbs
+           
+             taumeltbs=taumeltbs0*exp(-max(0.,(zt(i)-RTT)/(tbsmelt-RTT)))
+             dqmelt=min(zqbs(i),zqbs(i)/taumeltbs*dtime)
+             maxdqmelt= max(RCPD*(1.0+RVTMP2*(zq(i)+zqbs(i)))*(zt(i)-RTT)/(RLMLT+RLVTT),0.)
+             dqmelt=min(max(dqmelt,0.),maxdqmelt)
+             ! q update, melting + evaporation
+             zq(i) = zq(i) + dqmelt
+             ! temp update melting + evaporation
+             zt(i) = zt(i) - dqmelt * (RLMLT+RLVTT)/RCPD/(1.0+RVTMP2*(zq(i)+zqbs(i)))
+             ! qbs update melting + evaporation
+             zqbs(i)=zqbs(i)-dqmelt
+
           ELSE 
               ! negative celcius temperature     
               ! Sublimation scheme   
 
-              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))
+              rhoair=pplay(i,k)/zt(i)/RD 
+              dqsub = 1./rhoair*coef_eva_bs*(1.0-zq(i)/qsi(i))*((rhoair*zvelo(i)*zqbs(i))**expo_eva_bs)*dtime
+              dqsub = MAX(0.0,MIN(dqsub,zqbs(i)))
 
               ! 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)
+              maxdqsub = MAX(0.0, qsi(i)-zq(i))
+              dqsub = MIN(dqsub,maxdqsub)
 
 
               ! vapor, temperature, precip fluxes update following sublimation
-              zq(i) = zq(i) - (sedimn(i)-sedim(i))*(RG/(paprs(i,k)-paprs(i,k+1)))*dtime                
-              zq(i) = max(0., zq(i))
-              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)))
-
-               sedim(i)=sedimn(i)
-               zqbsi(i)=zqbs(i)
-
-               ! now sedimentation scheme with an exact numerical resolution 
-               ! (correct if fall velocity is constant)
-
-               zqbs(i) = zqbsi(i)*exp(-velo(i,k)/zdz(i)*dtime)+sedim(i)/zrho(i,k)/velo(i,k)
-               zqbs(i) = max(zqbs(i),0.)
-
-               ! flux update
-               sedim(i) = sedim(i) + zrho(i,k)*zdz(i)/dtime*(zqbsi(i)-zqbs(i))
-               sedim(i) = max(0.,sedim(i))
-
-
-           ENDIF
-
-
-           ! if qbs<qbsmin, sublimate or melt and evaporate qbs
-           ! see Gerber et al. 2023, JGR Atmos for the choice of qbsmin
-           IF (zqbs(i) .LT. qbsmin) THEN                       
+              zq(i) = zq(i) + dqsub     
+              zqbs(i)=zqbs(i)-dqsub           
+              zt(i) = zt(i) - dqsub*RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zqbs(i)))
+
+          ENDIF
+
+        ! Sedimentation scheme
+           
+        rhoair  = pplay(i,k) / zt(i) / RD
+        dz      = (paprs(i,k)-paprs(i,k+1)) / (rhoair*RG)
+        ! BS fall velocity assumed to be constant for now
+        zvelo(i) = fallv_bs
+
+        sedimn(i) = rhoair*zqbs(i)*zvelo(i)
+
+        ! exact numerical resolution of sedimentation
+        ! assuming fall velocity is constant 
+
+        zqbs(i) = zqbs(i)*exp(-zvelo(i)/dz*dtime)+sedim(i)/rhoair/zvelo(i)
+
+        ! flux update
+        sedim(i) = sedimn(i)
+
+
+
+        ! if qbs<qbsmin, sublimate or melt and evaporate qbs
+        ! see Gerber et al. 2023, JGR Atmos for the choice of qbsmin
+        IF (zqbs(i) .LT. qbsmin) THEN                       
               zq(i) = zq(i)+zqbs(i)
               IF (zt(i) .LT. RTT) THEN
@@ -221,19 +176,14 @@
               ENDIF
               zqbs(i)=0.
-           ENDIF
-
-
-    ENDDO ! loop on ngrid
-
-
+        ENDIF
 
 
     ! 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) - qv(i,k)
         dtemp_bs(i,k) = zt(i) - temp(i,k)
-    ENDDO
+
+    ENDDO ! Loop on ngrid
 
 
Index: /LMDZ6/trunk/libf/phylmd/physiq_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/physiq_mod.F90	(revision 4820)
+++ /LMDZ6/trunk/libf/phylmd/physiq_mod.F90	(revision 4821)
@@ -1859,6 +1859,5 @@
        CALL ratqs_ini(klon,klev,iflag_thermals,lunout,nbsrf,is_lic,is_ter,RG,RV,RD,RCPD,RLSTT,RLVTT,RTT)
        CALL lscp_ini(pdtphys,lunout,prt_level,ok_ice_sursat,iflag_ratqs,fl_cor_ebil,RCPD,RLSTT,RLVTT,RLMLT,RVTMP2,RTT,RD,RG,RPI)
-       CALL blowing_snow_ini(prt_level,lunout, &
-                             RCPD, RLSTT, RLVTT, RLMLT, &
+       CALL blowing_snow_ini(RCPD, RLSTT, RLVTT, RLMLT, &
                              RVTMP2, RTT,RD,RG)
        ! Test de coherence sur oc_cdnc utilisé uniquement par cloud_optics_prop
Index: /LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90	(revision 4820)
+++ /LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90	(revision 4821)
@@ -1951,6 +1951,5 @@
        CALL ratqs_ini(klon,klev,iflag_thermals,lunout,nbsrf,is_lic,is_ter,RG,RV,RD,RCPD,RLSTT,RLVTT,RTT)
        CALL lscp_ini(pdtphys,lunout,prt_level,ok_ice_sursat,iflag_ratqs,fl_cor_ebil,RCPD,RLSTT,RLVTT,RLMLT,RVTMP2,RTT,RD,RG,RPI)
-       CALL blowing_snow_ini(prt_level,lunout, &
-                             RCPD, RLSTT, RLVTT, RLMLT, &
+       CALL blowing_snow_ini(RCPD, RLSTT, RLVTT, RLMLT, &
                              RVTMP2, RTT,RD,RG)
        ! Test de coherence sur oc_cdnc utilisé uniquement par cloud_optics_prop
