Changeset 5871

Timestamp:

Nov 18, 2025, 11:04:21 AM (3 weeks ago)

Author:

evignon

Message:

slight changes on blowing snow

Location:

LMDZ6/branches/lmdz-snow/libf/phylmd

Files:

• lmdz_blowing_snow_ini.f90 (modified) (2 diffs)
• lmdz_blowing_snow_sublim_sedim.f90 (modified) (8 diffs)
• lmdz_call_blowing_snow.f90 (modified) (3 diffs)
• physiq_mod.F90 (modified) (1 diff)

LMDZ6/branches/lmdz-snow/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/branches/lmdz-snow/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)
 
 !==============================================================================
@@ -27 +27 @@
 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) :: qv
@@ -55 +56 @@
 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 +73 @@
 precip_bs(:)=0.
 bsfl(:,:)=0.
+zz(:,:)=0.
 
 
@@ -82 +84 @@
       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 +178 @@
         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 +226 @@
               ! 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 +241 @@
               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)

LMDZ6/branches/lmdz-snow/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/branches/lmdz-snow/libf/phylmd/physiq_mod.F90

@@ -3098 +3098 @@
     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)