Changeset 5912

Timestamp:

Dec 4, 2025, 3:16:33 PM (111 minutes ago)

Author:

evignon

Message:

changes in blowing snow routines to comply with the revised version of the paper

Location:

LMDZ6/trunk/libf

Files:

• phylmd/lmdz_blowing_snow_ini.f90 (modified) (2 diffs)
• phylmd/lmdz_blowing_snow_sublim_sedim.f90 (modified) (11 diffs)
• phylmd/lmdz_call_blowing_snow.f90 (modified) (3 diffs)
• phylmd/physiq_mod.F90 (modified) (1 diff)
• phylmdiso/physiq_mod.F90 (modified) (1 diff)

LMDZ6/trunk/libf/phylmd/lmdz_blowing_snow_ini.f90

@@ -70 +70 @@
          CALL getin_p('fallv_bs',fallv_bs)
 
-         coef_sub_bs =  0.1 
+         coef_sub_bs =  0.01 
          CALL getin_p('coef_sub_bs',coef_sub_bs)
 
@@ -79 +79 @@
          CALL getin_p('iflag_sedim_bs',iflag_sedim_bs)
 
-         r_bs=150.0e-6
+         r_bs=50.0e-6
          CALL getin_p('r_bs',r_bs)
 

LMDZ6/trunk/libf/phylmd/lmdz_blowing_snow_sublim_sedim.f90

@@ -2 +2 @@
 
 contains
-subroutine blowing_snow_sublim_sedim(ngrid,nlay,dtime,temp,qv,qb,pplay,paprs,dtemp_bs,dqv_bs,dqb_bs,bsfl,precip_bs)
+
+!==============================================================================        
+subroutine blowing_snow_sublim_sedim(ngrid,nlay,dtime,ustar,temp,qv,qb,pplay,paprs,dtemp_bs,dqv_bs,dqb_bs,bsfl,precip_bs)
 
 !==============================================================================
-! Routine that calculates the evaporation and sedimentation of blowing snow
-! inspired by what is done in lscp_mod
+! Routine that calculates the sublimation, melting and sedimentation 
+! of blowing snow
+! Reference: Vignon et al 2025, GMD, https://doi.org/10.5194/egusphere-2025-2871
+! 
 ! Etienne Vignon, October 2022
 !==============================================================================
@@ -12 +16 @@
 
 use lmdz_blowing_snow_ini, only : iflag_sublim_bs, iflag_sedim_bs, coef_sub_bs,RTT,RD,RG,fallv_bs
-use lmdz_blowing_snow_ini, only : qbmin, RCPD, RLSTT, RLMLT, RLVTT, RVTMP2, RV, RPI, tbsmelt, taumeltbs0, rhobs, r_bs
+use lmdz_blowing_snow_ini, only : qbmin, RCPD, RLSTT, RLMLT, RLVTT, RVTMP2, RV, RPI
+use lmdz_blowing_snow_ini, only : tbsmelt, taumeltbs0, rhobs, r_bs
 USE lmdz_lscp_tools, only : calc_qsat_ecmwf
 
@@ -25 +30 @@
 !=====
 
-integer, intent(in)                     :: ngrid,nlay
-real, intent(in)                        :: dtime
-real, intent(in), dimension(ngrid,nlay) :: temp
-real, intent(in), dimension(ngrid,nlay) :: qv
-real, intent(in), dimension(ngrid,nlay) :: qb
-real, intent(in), dimension(ngrid,nlay) :: pplay
-real, intent(in), dimension(ngrid,nlay+1) :: paprs
+integer, intent(in)                     :: ngrid ! number of horizontal grid points
+integer, intent(in)                     :: nlay  ! number of vertical layers
+real, intent(in)                        :: dtime ! physics time step [s]
+real, intent(in), dimension(ngrid)      :: ustar ! surface friction velocity [m/s]
+real, intent(in), dimension(ngrid,nlay) :: temp  ! temperature of the air [K]
+real, intent(in), dimension(ngrid,nlay) :: qv    ! specific content of water [kg/kg]
+real, intent(in), dimension(ngrid,nlay) :: qb    ! specific content of blowing snow [kg/kg]
+real, intent(in), dimension(ngrid,nlay) :: pplay ! pressure at the middle of layers [Pa]
+real, intent(in), dimension(ngrid,nlay+1) :: paprs ! pressure at the layer bottom interface [Pa]
 
 
@@ -39 +46 @@
 
 
-real, intent(out), dimension(ngrid,nlay) :: dtemp_bs
-real, intent(out), dimension(ngrid,nlay) :: dqv_bs
-real, intent(out), dimension(ngrid,nlay) :: dqb_bs
-real, intent(out), dimension(ngrid,nlay+1) :: bsfl
-real, intent(out), dimension(ngrid)      :: precip_bs
-
+real, intent(out), dimension(ngrid,nlay) :: dtemp_bs ! temperature tendency [K/s] 
+real, intent(out), dimension(ngrid,nlay) :: dqv_bs   ! water vapor tendency [kg/kg/s]
+real, intent(out), dimension(ngrid,nlay) :: dqb_bs   ! blowing snow mass tendancy [kg/kg/s]
+real, intent(out), dimension(ngrid,nlay+1) :: bsfl   ! vertical profile of blowing snow vertical flux [kg/m2/s]
+real, intent(out), dimension(ngrid)      :: precip_bs ! surface sedimentation flux of blowing snow [kg/s/m2]
 
 ! LOCAL
@@ -55 +61 @@
 real                                     :: maxdqvmelt, rhoair, dz, dqbsedim
 real                                     :: delta_p, b_p, a_p, c_p, c_sub, qvsub
-real                                     :: p0, T0, Dv, Aprime, Bprime, Ka
+real                                     :: p0, T0, Dv, Aprime, Bprime, Ka, radius0, radius
 real, dimension(ngrid)                   :: ztemp,zqv,zqb,zpres,qsi,dqsi,qsl,dqsl,qzero,sedim
 real, dimension(ngrid)                   :: zpaprsup, zpaprsdn, ztemp_up, zqb_up, zvelo
-real, dimension(ngrid,nlay)              :: temp_seri, qb_seri, qv_seri
+real, dimension(ngrid,nlay)              :: temp_seri, qb_seri, qv_seri, zz
 
 !++++++++++++++++++++++++++++++++++++++++++++++++++
@@ -72 +78 @@
 precip_bs(:)=0.
 bsfl(:,:)=0.
+zz(:,:)=0.
 
 
@@ -82 +89 @@
       qv_seri(i,k)=qv(i,k)
       qb_seri(i,k)=qb(i,k)
+      zz(i,k)=zz(i,k)+(paprs(i,k)-paprs(i,k+1)) / (pplay(i,k)/RD/temp(i,k)*RG)
    ENDDO
 ENDDO
@@ -175 +183 @@
         zpaprsup(i)=paprs(i,k+1)
         zpaprsdn(i)=paprs(i,k)
+
+    ! computation of blowing snow mean radius. Either constant value = r_bs if >0
+    ! or use of the height-dependent formulation from Sauter et al. 2013 and Saigger et al. 2024
+        IF (r_bs .GT. 0.) THEN
+            radius=r_bs
+        ELSE
+            radius0=0.5*(ustar(i)*(7.8e-6)/0.036+31.e-6)
+            radius=radius0*zz(i,k)**(-0.258)
+        ENDIF
+
     ENDDO
 
@@ -213 +231 @@
               ! Sublimation formulation for ice crystals from Pruppacher & Klett, Rutledge & Hobbs 1983
               ! assuming monodispered crystal distrib
-              ! dqb/dt_sub=-coef_sub_bs*(1-qv/qsi)*nc*8*r_bs/(Aprime+Bprime)
-              ! assuming Mi=rhobs*4/3*pi*r_bs**3 
-              ! rhoair qb=nc*Mi -> nc=rhoair qb/Mi
-              ! dqb/dt_sub=-coef_sub_bs*(1-qv/qsi)*6*rhoair*qb/(rhobs*pi*r_bs**2)/(Aprime+Bprime)
+              ! dqb/dt_sub=-coef_sub_bs*(1-qv/qsi)*nc*8*radius/(Aprime+Bprime)
+              ! assuming Mi=rhobs*4/3*pi*radius**3 
+              ! qb=nc*Mi -> nc=qb/Mi
+              ! dqb/dt_sub=-coef_sub_bs*(1-qv/qsi)*6*qb/(rhobs*pi*radius**2)/(Aprime+Bprime)
               ! dqb/dt_sub=-c_sub(1-qv/qsi)*qb
-              ! c_sub=coef_sub_bs*6*rhoair/(rhobs*pi*r_bs**2)/(Aprime+Bprime)
+              ! c_sub=coef_sub_bs*6/(rhobs*pi*radius**2)/(Aprime+Bprime)
               ! 
               ! Note the strong coupling between specific contents of water vapor and blowing snow during sublimation
@@ -228 +246 @@
               IF (zqv(i) .LT. qsi(i)) THEN
                  rhoair=zpres(i)/ztemp(i)/RD
-                 Dv=0.0001*0.211*(p0/zpres(i))*((ztemp(i)/RTT)**1.94) ! water vapor diffusivity in air, SI
+                 Dv=0.0001*0.211*(p0/zpres(i))*((ztemp(i)/RTT)**1.94)           ! water vapor diffusivity in air, SI
                  Ka=(5.69+0.017*(ztemp(i)-RTT))*1.e-5*100.*4.184                ! thermal conductivity of the air, SI
                  Aprime=RLSTT/Ka/ztemp(i)*(RLSTT/RV/ztemp(i) -1.)
                  Bprime=1./(rhoair*Dv*qsi(i)) 
-                 c_sub=coef_sub_bs*6.*rhoair/(rhobs*RPI*r_bs**2)/(Aprime+Bprime)
+                 c_sub=coef_sub_bs*3./(rhobs*radius**2)/(Aprime+Bprime)
                  c_p=-zqb(i)
                  b_p=1.+c_sub*dtime-c_sub*dtime/qsi(i)*zqb(i)-c_sub*dtime/qsi(i)*zqv(i)
@@ -293 +311 @@
 
 end subroutine blowing_snow_sublim_sedim
+!==============================================================================
+
+
+
 end module lmdz_blowing_snow_sublim_sedim

LMDZ6/trunk/libf/phylmd/lmdz_call_blowing_snow.f90

@@ -4 +4 @@
 contains
 
-subroutine call_blowing_snow_sublim_sedim(ngrid,nlay,dtime,temp,q,qbs,pplay,paprs, &
+subroutine call_blowing_snow_sublim_sedim(ngrid,nlay,dtime,ustar,temp,q,qbs,pplay,paprs, &
                                   dtemp_bs,dq_bs,dqbs_bs,bsfl,precip_bs)
 
@@ -15 +15 @@
 integer, intent(in)                     :: ngrid,nlay
 real, intent(in)                        :: dtime
+real, intent(in), dimension(ngrid)      :: ustar
 real, intent(in), dimension(ngrid,nlay) :: temp
 real, intent(in), dimension(ngrid,nlay) :: q
@@ -35 +36 @@
 
 
-call blowing_snow_sublim_sedim(ngrid,nlay,dtime,temp,q,qbs,pplay,paprs, &
+call blowing_snow_sublim_sedim(ngrid,nlay,dtime,ustar,temp,q,qbs,pplay,paprs, &
                                   dtemp_bs,dq_bs,dqbs_bs,bsfl,precip_bs)
 

LMDZ6/trunk/libf/phylmd/physiq_mod.F90

@@ -3106 +3106 @@
     IF (ok_bs) THEN
 
-     CALL call_blowing_snow_sublim_sedim(klon,klev,phys_tstep,t_seri,q_seri,qbs_seri,pplay,paprs, &
+     CALL call_blowing_snow_sublim_sedim(klon,klev,phys_tstep,zustar,t_seri,q_seri,qbs_seri,pplay,paprs, &
                                         d_t_bsss,d_q_bsss,d_qbs_bsss,bsfl,bs_fall)
 

LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90

@@ -3619 +3619 @@
     IF (ok_bs) THEN
 
-     CALL call_blowing_snow_sublim_sedim(klon,klev,phys_tstep,t_seri,q_seri,qbs_seri,pplay,paprs, &
-                                        d_t_bsss,d_q_bsss,d_qbs_bsss,bsfl,bs_fall)
+     CALL call_blowing_snow_sublim_sedim(klon,klev,phys_tstep,zustar,t_seri,q_seri,qbs_seri,pplay,paprs, &
+                                         d_t_bsss,d_q_bsss,d_qbs_bsss,bsfl,bs_fall)
 
      CALL add_phys_tend &