Index: LMDZ6/trunk/libf/phylmd/lmdz_blowing_snow_ini.f90 =================================================================== --- LMDZ6/trunk/libf/phylmd/lmdz_blowing_snow_ini.f90 (revision 5911) +++ LMDZ6/trunk/libf/phylmd/lmdz_blowing_snow_ini.f90 (revision 5912) @@ -70,5 +70,5 @@ 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,5 +79,5 @@ 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) Index: LMDZ6/trunk/libf/phylmd/lmdz_blowing_snow_sublim_sedim.f90 =================================================================== --- LMDZ6/trunk/libf/phylmd/lmdz_blowing_snow_sublim_sedim.f90 (revision 5911) +++ LMDZ6/trunk/libf/phylmd/lmdz_blowing_snow_sublim_sedim.f90 (revision 5912) @@ -2,9 +2,13 @@ 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,5 +16,6 @@ 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,11 +30,13 @@ !===== -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,10 +46,9 @@ -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,8 +61,8 @@ 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,4 +78,5 @@ precip_bs(:)=0. bsfl(:,:)=0. +zz(:,:)=0. @@ -82,4 +89,5 @@ 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,4 +183,14 @@ 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,10 +231,10 @@ ! 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,9 +246,9 @@ 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,3 +311,7 @@ end subroutine blowing_snow_sublim_sedim +!============================================================================== + + + end module lmdz_blowing_snow_sublim_sedim Index: LMDZ6/trunk/libf/phylmd/lmdz_call_blowing_snow.f90 =================================================================== --- LMDZ6/trunk/libf/phylmd/lmdz_call_blowing_snow.f90 (revision 5911) +++ LMDZ6/trunk/libf/phylmd/lmdz_call_blowing_snow.f90 (revision 5912) @@ -4,5 +4,5 @@ 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,4 +15,5 @@ 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,5 +36,5 @@ -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) Index: LMDZ6/trunk/libf/phylmd/physiq_mod.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/physiq_mod.F90 (revision 5911) +++ LMDZ6/trunk/libf/phylmd/physiq_mod.F90 (revision 5912) @@ -3106,5 +3106,5 @@ 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) Index: LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90 =================================================================== --- LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90 (revision 5911) +++ LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90 (revision 5912) @@ -3619,6 +3619,6 @@ 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 &