Changeset 4916

Timestamp:

Apr 22, 2024, 9:06:52 PM (7 months ago)

Author:

evignon

Message:

correction formulation de l'erosion de la neige soufflee suite aux investigations de Nicolas Chiabrando,
prise en compte d'un temps d'erosion de la neige fraiche nouvellement accumulee

Location:

LMDZ6/trunk/libf/phylmd

Files:

• lmdz_blowing_snow_sublim_sedim.F90 (modified) (3 diffs)
• pbl_surface_mod.F90 (modified) (1 diff)
• surf_landice_mod.F90 (modified) (8 diffs)

LMDZ6/trunk/libf/phylmd/lmdz_blowing_snow_sublim_sedim.F90

@@ -74 +74 @@
 
 
-Ka=2.4e-2      ! thermal conductivity of the air, SI
 p0=101325.0    ! ref pressure
 
@@ -229 +228 @@
               IF (zqv(i) .LT. qsi(i)) THEN
                  rhoair=zpres(i)/ztemp(i)/RD
-                 !diffusivity of water vapor
                  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
@@ -241 +239 @@
                  dqbsub=(-b_p+sqrt(delta_p))/(2.*a_p) - zqb(i)       
                  dqbsub = MIN(0.0,MAX(dqbsub,-zqb(i)))
-
                  ! Sublimation limit: we ensure that the whole mesh does not exceed saturation wrt ice
                  maxdqbsub = MAX(0.0, qsi(i)-zqv(i))

LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90

@@ -1915 +1915 @@
 ! For blowing snow:
     IF (ok_bs) THEN
-     ! following Bintanja et al 2000, part II
+     ! following Bintanja et al 2000, part II and Vionnet V PhD thesis
      ! we assume that the eddy diffsivity coefficient for
-     ! suspended particles is larger than Km by a factor zeta_bs
-     ! which is equal to 3 by default
+     ! suspended particles is a fraction of Kh 
      do k=1,klev
         do j=1,knon
-           ycoefqbs(j,k)=ycoefm(j,k)*zeta_bs 
+           ycoefqbs(j,k)=ycoefh(j,k)*zeta_bs 
         enddo
      enddo

LMDZ6/trunk/libf/phylmd/surf_landice_mod.F90

@@ -142 +142 @@
     REAL                   :: esalt
     REAL                   :: lambdasalt,fluxsalt, csalt, nunu, aa, bb, cc
-    REAL                   :: tau_dens, maxerosion, fluxbs_2
+    REAL                   :: tau_dens, maxerosion
     REAL, DIMENSION(klon)  :: ws1, rhod, rhos, ustart0, ustart, qsalt, hsalt
+    REAL, DIMENSION(klon)  :: fluxbs_1, fluxbs_2, bsweight_fresh
+    LOGICAL, DIMENSION(klon) :: ok_remaining_freshsnow
 
 ! End definition
@@ -373 +375 @@
        ! 1st step: erosion of fresh snow accumulated during the time step
        do j=1, knon
-
+       if (precip_snow(j) .GT. 0.) then
            rhos(j)=rhofresh_bs
            ! blowing snow flux formula used in MAR
@@ -392 +394 @@
            maxerosion=min(precip_snow(j),hsalt(j)*qsalt(j)*rhod(j)/tau_eqsalt_bs)
 
-           fluxbs(j)=rhod(j)*ws1(j)*cdragm(j)*(AcoefQBS(j)-qsalt(j)) &
-                   / (1.-rhod(j)*ws1(j)*cdragm(j)*BcoefQBS(j)*dtime)
-           !fluxbs(j)= rhod*ws1(j)*cdragm(j)*(qbs1(j)-qsalt(j))
-           fluxbs(j)=max(-maxerosion,fluxbs(j))
+           fluxbs_1(j)=rhod(j)*ws1(j)*cdragh(j)*zeta_bs*(AcoefQBS(j)-qsalt(j)) &
+                   / (1.-rhod(j)*ws1(j)*cdragh(j)*zeta_bs*BcoefQBS(j)*dtime)
+           fluxbs_1(j)=max(-maxerosion,fluxbs_1(j))
+
+           if (precip_snow(j) .gt. abs(fluxbs_1(j))) then
+               ok_remaining_freshsnow(j)=.true.
+               bsweight_fresh(j)=1.
+           else
+               ok_remaining_freshsnow(j)=.false.
+               bsweight_fresh(j)=exp(-(abs(fluxbs_1(j))-precip_snow(j))/precip_snow(j))
+           endif
+       else
+           ok_remaining_freshsnow(j)=.false.
+           fluxbs_1(j)=0.
+           bsweight_fresh(j)=0.
+       endif
        enddo
 
@@ -401 +415 @@
        ! we now compute the snow age of the overlying layer (snow surface after erosion of the fresh snow accumulated during the time step)
        ! this is done through the routine albsno
-       CALL albsno(klon,knon,dtime,agesno(:),alb_neig(:), precip_snow(:)+fluxbs(:))
+       CALL albsno(klon,knon,dtime,agesno(:),alb_neig(:), precip_snow(:)+fluxbs_1(:))
 
        ! 2nd step:
@@ -407 +421 @@
        ! which depends on surface snow density
        do j = 1, knon
+        if (ok_remaining_freshsnow(j)) then
+           fluxbs_2(j)=0.
+        else
+           ! we start eroding the underlying layer
            ! estimation of snow density
            ! snow density increases with snow age and
@@ -418 +436 @@
            ! computation of qbs at the top of the saltation layer
            ! default formulation from MAR model (Amory et al. 2021, Gallee et al. 2001)
-           esalt=c_esalt_bs*ustar(j)
+           esalt=1./(c_esalt_bs*max(1.e-6,ustar(j)))
            hsalt(j)=0.08436*(max(1.e-6,ustar(j))**1.27)
            qsalt(j)=(max(ustar(j)**2-ustart(j)**2,0.))/(RG*hsalt(j))*esalt
@@ -428 +446 @@
            ! (rho*qsalt*hsalt)
            maxerosion=hsalt(j)*qsalt(j)*rhod(j)/tau_eqsalt_bs
-           fluxbs_2=rhod(j)*ws1(j)*cdragm(j)*(AcoefQBS(j)-qsalt(j)) &
-                   / (1.-rhod(j)*ws1(j)*cdragm(j)*BcoefQBS(j)*dtime)
-           !fluxbs_2= rhod*ws1(j)*cdragm(j)*(qbs1(j)-qsalt(j))
-           fluxbs_2=max(-maxerosion,fluxbs_2)
-           fluxbs(j)=fluxbs(j)+fluxbs_2
+           fluxbs_2(j)=rhod(j)*ws1(j)*cdragh(j)*zeta_bs*(AcoefQBS(j)-qsalt(j)) &
+                   / (1.-rhod(j)*ws1(j)*cdragh(j)*zeta_bs*BcoefQBS(j)*dtime)
+           fluxbs_2(j)=max(-maxerosion,fluxbs_2(j))
+         endif
        enddo
 
-       ! for outputs       
+
+
+
+       ! final flux and outputs       
         do j=1, knon
+              ! total flux is the erosion of fresh snow + 
+              ! a fraction of the underlying snow (if all the fresh snow has been eroded)
+              ! the calculation of the fraction is quite delicate since we do not know
+              ! how much time was needed to erode the fresh snow. We assume that this time
+              ! is dt*exp(-(abs(fluxbs1)-precipsnow)/precipsnow)=dt*bsweight_fresh
+
+              fluxbs(j)=fluxbs_1(j)+fluxbs_2(j)*(1.-bsweight_fresh(j))
               i = knindex(j)
               zxustartlic(i) = ustart(j)
@@ -444 +471 @@
 
 
-  else
+  else ! not ok_bs
   ! those lines are useful to calculate the snow age
        CALL albsno(klon,knon,dtime,agesno(:),alb_neig(:), precip_snow(:))