Changeset 4672 for LMDZ6/trunk/libf/phylmdiso

Timestamp:

Sep 4, 2023, 4:35:35 PM (17 months ago)

Author:

evignon

Message:

inclusion des recents developpements sur la neige soufflee

File:

• LMDZ6/trunk/libf/phylmdiso/surf_landice_mod.F90 (modified) (4 diffs)

LMDZ6/trunk/libf/phylmdiso/surf_landice_mod.F90

@@ -185 +185 @@
     REAL,DIMENSION(klon) :: precip_totsnow, evap_totsnow
     REAL, DIMENSION (klon,6) :: alb6
-    REAL                   :: rho0, rhoice, ustart0, hsalt, esalt, rhod
+    REAL                   :: rho0, rhoice, ustart0,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, qsalt
+    REAL, DIMENSION(klon)  :: ws1, rhos, ustart, qsalt, hsalt, snowini
+    REAL                   :: maxerosion ! in kg/s
+
 ! End definition
 !****************************************************************************************
@@ -453 +455 @@
        ustart0 = 0.211
        rhoice = 920.0
-       rho0 = 200.0
+       rho0 = 300.0
        rhomax=450.0
-       rhohard=400.0
+       rhohard=450.0
        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
+       tau_densmin=21600.0    ! 1/4 days according to in situ obs by C. Amory during blowing snow +
+                              ! Marshall et al. 1999 (snow densification during rain)
 
        ! computation of threshold friction velocity
        ! which depends on surface snow density
-       do i = 1, knon
+       do j = 1, knon
            ! estimation of snow density
            ! snow density increases with snow age and
-           ! increases even faster in case of sedimentation of blowing snow
-           tau_dens=max(tau_densmin, tau_dens0*exp(-abs(precip_bs(i))/pbst_bs-abs(precip_rain(i))/prt_bs))
-           rhos(i)=rho0+(rhohard-rho0)*(1.-exp(-agesno(i)*86400.0/tau_dens))
+           ! increases even faster in case of sedimentation of blowing snow or rain
+           tau_dens=max(tau_densmin, & tau_dens0*exp(-abs(precip_bs(j))/pbst_bs-abs(precip_rain(j))/prt_bs) &
+                    *exp(-max(tsurf(j)-272.0,0.)))
+           rhos(j)=rho0+(rhohard-rho0)*(1.-exp(-agesno(j)*86400.0/tau_dens))
            ! blowing snow flux formula used in MAR
-           ws1(i)=(u1(i)**2+v1(i)**2)**0.5
-           ustar(i)=(cdragm(i)*(u1(i)**2+v1(i)**2))**0.5
-           ustart(i)=ustart0*exp(max(rhoice/rho0-rhoice/rhos(i),0.))*exp(max(0.,rhos(i)-rhomax)) 
+           ws1(j)=(u1(j)**2+v1(j)**2)**0.5
+           ustar(j)=(cdragm(j)*(u1(j)**2+v1(j)**2))**0.5
+           ustart(j)=ustart0*exp(max(rhoice/rho0-rhoice/rhos(j),0.))*exp(max(0.,rhos(j)-rhomax)) 
            ! we have multiplied by exp to prevent erosion when rhos>rhomax (usefull till
            ! rhohard<450)
@@ -486 +490 @@
           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)) * &
+          do j =1, knon
+               rhod=p1lay(j)/RD/temp_air(j)
+               nunu=max(ustar(j)/ustart(j),1.e-3)
+               fluxsalt=rhod/RG*(ustar(j)**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.)
+               csalt=fluxsalt/(2.8*ustart(j))
+               hsalt(j)=0.08436*ustar(j)**1.27
+               qsalt(j)=1./rhod*csalt*lambdasalt*RG/(max(ustar(j)**2,1E-6)) &
+                       * exp(-lambdasalt*RG*hsalt(j)/max(ustar(j)**2,1E-6))
+               qsalt(j)=max(qsalt(j),0.)
           enddo
 
@@ -501 +505 @@
        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)
+          do j=1, knon
+              esalt=1./(3.25*max(ustar(j),0.001))
+              hsalt(j)=0.08436*ustar(j)**1.27
+              qsalt(j)=(max(ustar(j)**2-ustart(j)**2,0.))/(RG*hsalt(j))*esalt
+              !ep=qsalt*cdragm(j)*sqrt(u1(j)**2+v1(j)**2)
           enddo
        endif
 
-        ! calculation of erosion (emission flux towards the first atmospheric level)
+        ! calculation of erosion (flux positive towards the surface here) 
         ! 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*cdragm(i)*BcoefQBS(i)*dtime)
-              !fluxbs(i)= zeta_bs*rhod*ws1(i)*cdragm(i)*(qbs1(i)-qsalt(i))
-       enddo
-
-       ! for outputs
-       do j = 1, knon
-          i = knindex(j)
-          zxustartlic(i) = ustart(j)
-          zxrhoslic(i) = rhos(j)
+       do j=1, knon
+              i = knindex(j)
+              rhod=p1lay(j)/RD/temp_air(j)
+              ! Nemoto and Nishimura 2004 show that steady-state saltation is achieved within an interval of about 10s
+              ! we thus prevent snowerosion (snow particle transfer from the saltation layer to the first model level)
+              ! to exceed (in intensity) integrated over 10s to exceed the total mass of snow particle in the saltation layer
+              ! (rho*qsalt*hsalt)
+              maxerosion=hsalt(j)*qsalt(j)*rhod/10.
+              fluxbs(j)=rhod*ws1(j)*cdragm(j)*(AcoefQBS(j)-qsalt(j)) &
+                       / (1.-rhod*ws1(j)*cdragm(j)*BcoefQBS(j)*dtime)
+              !fluxbs(j)= rhod*ws1(j)*cdragm(j)*(qbs1(j)-qsalt(j))
+              fluxbs(j)=max(-maxerosion,fluxbs(j))
+
+              ! for outputs
+
+              zxustartlic(i) = ustart(j)
+              zxrhoslic(i) = rhos(j)
        enddo