Changeset 4529 for LMDZ6/trunk/libf/phylmd

Timestamp:

May 11, 2023, 10:53:52 AM (3 years ago)

Author:

evignon

Message:

nouvelle formulation du flux d'emission de neige soufflee

Location:

LMDZ6/trunk/libf/phylmd

Files:

• blowing_snow_ini_mod.F90 (modified) (3 diffs)
• pbl_surface_mod.F90 (modified) (1 diff)
• surf_landice_mod.F90 (modified) (6 diffs)

LMDZ6/trunk/libf/phylmd/blowing_snow_ini_mod.F90

@@ -8 +8 @@
    real coef_eva_bs, expo_eva_bs, fallv_bs, zeta_bs
    real prt_bs, pbst_bs, qbst_bs
-   integer :: niter_bs
+   integer :: niter_bs, iflag_saltation_bs
 
    !$OMP THREADPRIVATE(prt_level, lunout)
@@ -14 +14 @@
    !$OMP THREADPRIVATE(coef_eva_bs, expo_eva_bs, fallv_bs, zeta_bs)
    !$OMP THREADPRIVATE(pbst_bs, prt_bs, qbst_bs)
-   !$OMP THREADPRIVATE(niter_bs)
+   !$OMP THREADPRIVATE(niter_bs, iflag_saltation_bs)
 
       contains
@@ -67 +67 @@
          CALL getin_p('niter_bs',niter_bs)
 
+         iflag_saltation_bs=0
+         CALL getin_p('iflag_saltation_bs',iflag_saltation_bs)
+
 
       end subroutine blowing_snow_ini

LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90

@@ -2147 +2147 @@
                   AcoefH, AcoefQ, BcoefH, BcoefQ, &
                   AcoefU, AcoefV, BcoefU, BcoefV, &
+                  AcoefQBS, BcoefQBS, &
                   ypsref, yu1, yv1, ygustiness, yrugoro, pctsrf, &
                   ysnow, yqsurf, yqsol,yqbs1, yagesno, &

LMDZ6/trunk/libf/phylmd/surf_landice_mod.F90

@@ -15 +15 @@
        AcoefH, AcoefQ, BcoefH, BcoefQ, &
        AcoefU, AcoefV, BcoefU, BcoefV, &
+       AcoefQBS, BcoefQBS, &
        ps, u1, v1, gustiness, rugoro, pctsrf, &
        snow, qsurf, qsol, qbs1, agesno, &
@@ -34 +35 @@
 !FC
     USE ioipsl_getin_p_mod, ONLY : getin_p
-    USE blowing_snow_ini_mod, ONLY : zeta_bs, pbst_bs, prt_bs
+    USE blowing_snow_ini_mod, ONLY : zeta_bs, pbst_bs, prt_bs, iflag_saltation_bs
 
 #ifdef CPP_INLANDSIS
@@ -62 +63 @@
     REAL, DIMENSION(klon), INTENT(IN)             :: BcoefH, BcoefQ
     REAL, DIMENSION(klon), INTENT(IN)             :: AcoefU, AcoefV, BcoefU, BcoefV
+    REAL, DIMENSION(klon), INTENT(IN)             :: AcoefQBS, BcoefQBS
     REAL, DIMENSION(klon), INTENT(IN)             :: ps
     REAL, DIMENSION(klon), INTENT(IN)             :: u1, v1, gustiness, qbs1
@@ -138 +140 @@
     REAL,DIMENSION(klon) :: precip_totsnow, evap_totsnow
     REAL, DIMENSION (klon,6) :: alb6
-    REAL                   :: rho0, rhoice, ustart0, hsalt, esalt, qsalt
+    REAL                   :: rho0, rhoice, ustart0, hsalt, esalt, rhod
+    REAL                   :: lambdasalt,fluxsalt, csalt, nunu, aa, bb, cc
     REAL                   :: tau_dens, tau_dens0, tau_densmin, rhomax, rhohard
-    REAL, DIMENSION(klon)  :: ws1, rhos, ustart
+    REAL, DIMENSION(klon)  :: ws1, rhos, ustart, qsalt
 ! End definition
 !****************************************************************************************
@@ -366 +369 @@
        tau_dens0=86400.0*10.  ! 10 days by default, in s
        tau_densmin=86400.0 ! 1 days according to in situ obs by C. Amory
+
+       ! computation of threshold friction velocity
+       ! which depends on surface snow density
        do i = 1, knon
            ! estimation of snow density
@@ -378 +384 @@
            ! we have multiplied by exp to prevent erosion when rhos>rhomax (usefull till
            ! rhohard<450)
-           esalt=1./(3.25*max(ustar(i),0.001))
-           hsalt=0.08436*ustar(i)**1.27
-           qsalt=(max(ustar(i)**2-ustart(i)**2,0.))/(RG*hsalt)*esalt
-           !ep=qsalt*cdragm(i)*sqrt(u1(i)**2+v1(i)**2)
-           !Pay attention that this is an explicit formulation for the surface flux
-           !Formulation fully consistent with the implicit resolution of the turbulent
-           !mixing scheme remains to be coded
-           fluxbs(i)= zeta_bs*p1lay(i)/RD/temp_air(i)*ws1(i)*cdragm(i)*(qbs1(i)-qsalt)
+       enddo
+       
+       ! computation of qbs at the top of the saltation layer
+       ! two formulations possible
+       ! pay attention that qbs is a mixing ratio and has to be converted
+       ! to specific content
+       
+       if (iflag_saltation_bs .eq. 1) then
+       ! expression from CRYOWRF (Sharma et al. 2022)
+          aa=2.6
+          bb=2.5
+          cc=2.0
+          lambdasalt=0.45
+          do i =1, knon
+               rhod=p1lay(i)/RD/temp_air(i)
+               nunu=max(ustar(i)/ustart(i),1.e-3)
+               fluxsalt=rhod/RG*(ustar(i)**3)*(1.-nunu**(-2)) * &
+                        (aa+bb*nunu**(-2)+cc*nunu**(-1))  
+               csalt=fluxsalt/(2.8*ustart(i))
+               hsalt=0.08436*ustar(i)**1.27
+               qsalt(i)=1./rhod*csalt*lambdasalt*RG/(max(ustar(i)**2,1E-6)) &
+                       * exp(-lambdasalt*RG*hsalt/max(ustar(i)**2,1E-6))
+               qsalt(i)=max(qsalt(i),0.)
+          enddo
+
+
+       else
+       ! default formulation from MAR model (Amory et al. 2021, Gallee et al. 2001)        
+          do i=1, knon
+              esalt=1./(3.25*max(ustar(i),0.001))
+              hsalt=0.08436*ustar(i)**1.27
+              qsalt(i)=(max(ustar(i)**2-ustart(i)**2,0.))/(RG*hsalt)*esalt
+              !ep=qsalt*cdragm(i)*sqrt(u1(i)**2+v1(i)**2)
+          enddo
+       endif
+
+        ! calculation of erosion (emission flux towards the first atmospheric level)
+        ! consistent with implicit resolution of turbulent mixing equation
+       do i=1, knon
+              rhod=p1lay(i)/RD/temp_air(i)
+              fluxbs(i)=rhod*ws1(i)*cdragm(i)*zeta_bs*(AcoefQBS(i)-qsalt(i)) &
+                       / (1.-rhod*ws1(i)*zeta_bs*BcoefQBS(i)*dtime)
+              !fluxbs(i)= zeta_bs*rhod*ws1(i)*cdragm(i)*(qbs1(i)-qsalt(i))
        enddo