Changeset 4835 for LMDZ6/trunk/libf/phylmdiso/surf_landice_mod.F90

Timestamp:

Feb 29, 2024, 7:42:12 PM (3 months ago)

Author:

evignon

Message:

commission pour la suite du travail sur la mise a jour
de la param de neige soufflee

File:

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

LMDZ6/trunk/libf/phylmdiso/surf_landice_mod.F90

@@ -35 +35 @@
     USE cpl_mod,          ONLY : cpl_send_landice_fields
     USE calcul_fluxs_mod
-    USE phys_local_var_mod, ONLY : zxrhoslic, zxustartlic
+    USE phys_local_var_mod, ONLY : zxrhoslic, zxustartlic, zxqsaltlic
     USE phys_output_var_mod, ONLY : snow_o,zfra_o
 #ifdef ISO   
@@ -50 +50 @@
 !FC
     USE ioipsl_getin_p_mod, ONLY : getin_p
-    USE lmdz_blowing_snow_ini, ONLY : zeta_bs, pbst_bs, prt_bs, iflag_saltation_bs
-
+    USE lmdz_blowing_snow_ini, ONLY : c_esalt_bs, zeta_bs, pbst_bs, prt_bs, rhoice_bs, rhohard_bs
+    USE lmdz_blowing_snow_ini, ONLY : rhofresh_bs, tau_eqsalt_bs, tau_dens0_bs, tau_densmin_bs
 #ifdef CPP_INLANDSIS
     USE surf_inlandsis_mod,  ONLY : surf_inlandsis
@@ -57 +57 @@
 
     USE indice_sol_mod
-
 
 !    INCLUDE "indicesol.h"
@@ -185 +184 @@
     REAL,DIMENSION(klon) :: precip_totsnow, evap_totsnow
     REAL, DIMENSION (klon,6) :: alb6
-    REAL                   :: rho0, rhoice, ustart0,esalt, rhod
+    REAL                   :: esalt
     REAL                   :: lambdasalt,fluxsalt, csalt, nunu, aa, bb, cc
-    REAL                   :: tau_dens, tau_dens0, tau_densmin, rhomax, rhohard
-    REAL, DIMENSION(klon)  :: ws1, rhos, ustart, qsalt, hsalt, snowini
-    REAL                   :: maxerosion ! in kg/s
+    REAL                   :: tau_dens, maxerosion, fluxbs_2
+    REAL, DIMENSION(klon)  :: ws1, rhod, rhos, ustart0, ustart, qsalt, hsalt
 
 ! End definition
@@ -420 +418 @@
 !
 !****************************************************************************************
- 
-    CALL albsno(klon,knon,dtime,agesno(:),alb_neig(:), precip_snow(:))  
-
-
-! EV: following lines are obsolete since we set alb1 and alb2 to constant values
-! I therefore comment them
-!    alb1(1:knon) = alb_neig(1:knon)*zfra(1:knon) + &
-!         0.6 * (1.0-zfra(1:knon))
+
 !
 !IM: plusieurs choix/tests sur l'albedo des "glaciers continentaux"
@@ -450 +441 @@
 
 
-
-  ! Simple blowing snow param
-  if (ok_bs) then
-       ustart0 = 0.211
-       rhoice = 920.0
-       rho0 = 300.0
-       rhomax=450.0
-       rhohard=450.0
-       tau_dens0=86400.0*10.  ! 10 days by default, in s
-       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)
-
+!****************************************************************************************
+! Simple blowing snow param
+!****************************************************************************************
+! we proceed in 2 steps:
+! first we erode - if possible -the accumulated snow during the time step
+! then we update the density of the underlying layer and see if we can also erode
+! this layer
+
+
+   if (ok_bs) then
+       fluxbs(:)=0.
+       do j=1,klon
+          ws1(j)=(u1(j)**2+v1(j)**2)**0.5
+          ustar(j)=(cdragm(j)*(u1(j)**2+v1(j)**2))**0.5
+          rhod(j)=p1lay(j)/RD/temp_air(j)
+          ustart0(j) =(log(2.868)-log(1.625))/0.085*sqrt(cdragm(j))
+       enddo
+
+       ! 1st step: erosion of fresh snow accumulated during the time step
+       do j=1, knon
+
+           rhos(j)=rhofresh_bs
+           ! blowing snow flux formula used in MAR
+           ustart(j)=ustart0(j)*exp(max(rhoice_bs/rhofresh_bs-rhoice_bs/rhos(j),0.))*exp(max(0.,rhos(j)-rhohard_bs))
+           ! we have multiplied by exp to prevent erosion when rhos>rhohard_bs 
+           ! computation of qbs at the top of the saltation layer
+           ! default formulation from MAR model (Amory et al. 2021, Gallee et al. 2001)
+           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
+           ! calculation of erosion (flux positive towards the surface here) 
+           ! consistent with implicit resolution of turbulent mixing equation
+           ! Nemoto and Nishimura 2004 show that steady-state saltation is achieved within a time tau_eqsalt_bs of about 10s
+           ! we thus prevent snowerosion (snow particle transfer from the saltation layer to the first model level)
+           ! integrated over tau_eqsalt_bs to exceed the total mass of snow particle in the saltation layer
+           ! (rho*qsalt*hsalt)
+           ! during this first step we also lower bound the erosion to the amount of fresh snow accumulated during the time step
+           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))
+       enddo
+
+
+       ! 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(:))
+
+       ! 2nd step:
        ! computation of threshold friction velocity
        ! which depends on surface snow density
@@ -468 +498 @@
            ! snow density increases with snow age and
            ! 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))
+           tau_dens=max(tau_densmin_bs, tau_dens0_bs*exp(-abs(precip_bs(j))/pbst_bs - & 
+                    abs(precip_rain(j))/prt_bs)*exp(-max(tsurf(j)-RTT,0.)))
+           rhos(j)=rhofresh_bs+(rhohard_bs-rhofresh_bs)*(1.-exp(-agesno(j)*86400.0/tau_dens))
            ! blowing snow flux formula used in MAR
-           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)
+           ustart(j)=ustart0(j)*exp(max(rhoice_bs/rhofresh_bs-rhoice_bs/rhos(j),0.))*exp(max(0.,rhos(j)-rhohard_bs)) 
+           ! we have multiplied by exp to prevent erosion when rhos>rhohard_bs 
+           ! 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)
+           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
+           ! calculation of erosion (flux positive towards the surface here) 
+           ! consistent with implicit resolution of turbulent mixing equation
+           ! Nemoto and Nishimura 2004 show that steady-state saltation is achieved within a time tau_eqsalt_bs of about 10s
+           ! we thus prevent snowerosion (snow particle transfer from the saltation layer to the first model level)
+           ! integrated over tau_eqsalt_bs to exceed the total mass of snow particle in the saltation layer
+           ! (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
        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 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(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
-
-
-       else
-       ! default formulation from MAR model (Amory et al. 2021, Gallee et al. 2001)        
-          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 (flux positive towards the surface here) 
-        ! consistent with implicit resolution of turbulent mixing equation
-       do j=1, knon
+
+       ! for outputs       
+        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
-
-  endif
-
-
-
-!****************************************************************************************
-! Calculate surface snow amount
+              zxqsaltlic(i)=qsalt(j)
+        enddo
+
+
+  else
+  ! those lines are useful to calculate the snow age
+       CALL albsno(klon,knon,dtime,agesno(:),alb_neig(:), precip_snow(:))
+
+  endif ! if ok_bs
+
+
+
+
+
+
+!****************************************************************************************
+! Calculate snow amount
 !    
 !****************************************************************************************