Changeset 4669 for LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmd/surf_landice_mod.F90

Timestamp:

Sep 4, 2023, 10:17:16 AM (15 months ago)

Author:

Laurent Fairhead

Message:

Merged with trunk revision 4586 corresponding to june 2023 testing

Location:

LMDZ6/branches/LMDZ_cdrag_LSCE

Files:

• . (modified) (1 prop)
• libf/phylmd/surf_landice_mod.F90 (modified) (13 diffs)

LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmd/surf_landice_mod.F90

@@ -12 +12 @@
        rmu0, lwdownm, albedo, pphi1, &
        swnet, lwnet, tsurf, p1lay, &
-       cdragh, cdragm, precip_rain, precip_snow, temp_air, spechum, &
+       cdragh, cdragm, precip_rain, precip_snow, precip_bs, temp_air, spechum, &
        AcoefH, AcoefQ, BcoefH, BcoefQ, &
        AcoefU, AcoefV, BcoefU, BcoefV, &
+       AcoefQBS, BcoefQBS, &
        ps, u1, v1, gustiness, rugoro, pctsrf, &
-       snow, qsurf, qsol, agesno, &
-       tsoil, z0m, z0h, SFRWL, alb_dir, alb_dif, evap, fluxsens, fluxlat, &
+       snow, qsurf, qsol, qbs1, agesno, &
+       tsoil, z0m, z0h, SFRWL, alb_dir, alb_dif, evap, fluxsens, fluxlat, fluxbs, &
        tsurf_new, dflux_s, dflux_l, &
        alt, slope, cloudf, &
@@ -30 +31 @@
     USE cpl_mod,          ONLY : cpl_send_landice_fields
     USE calcul_fluxs_mod
+    USE phys_local_var_mod, ONLY : zxrhoslic, zxustartlic
     USE phys_output_var_mod, ONLY : snow_o,zfra_o
 !FC
     USE ioipsl_getin_p_mod, ONLY : getin_p
-
+    USE blowing_snow_ini_mod, ONLY : zeta_bs, pbst_bs, prt_bs, iflag_saltation_bs
 
 #ifdef CPP_INLANDSIS
@@ -56 +58 @@
     REAL, DIMENSION(klon), INTENT(IN)             :: p1lay
     REAL, DIMENSION(klon), INTENT(IN)             :: cdragh, cdragm
-    REAL, DIMENSION(klon), INTENT(IN)             :: precip_rain, precip_snow
+    REAL, DIMENSION(klon), INTENT(IN)             :: precip_rain, precip_snow, precip_bs
     REAL, DIMENSION(klon), INTENT(IN)             :: temp_air, spechum
     REAL, DIMENSION(klon), INTENT(IN)             :: AcoefH, AcoefQ
     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
+    REAL, DIMENSION(klon), INTENT(IN)             :: u1, v1, gustiness, qbs1
     REAL, DIMENSION(klon), INTENT(IN)             :: rugoro
     REAL, DIMENSION(klon,nbsrf), INTENT(IN)       :: pctsrf
@@ -94 +97 @@
 !albedo SB <<<
     REAL, DIMENSION(klon), INTENT(OUT)            :: evap, fluxsens, fluxlat
+    REAL, DIMENSION(klon), INTENT(OUT)            :: fluxbs
     REAL, DIMENSION(klon), INTENT(OUT)            :: tsurf_new
     REAL, DIMENSION(klon), INTENT(OUT)            :: dflux_s, dflux_l      
@@ -134 +138 @@
 
     REAL,DIMENSION(klon) :: alb1,alb2
+    REAL,DIMENSION(klon) :: precip_totsnow, evap_totsnow
     REAL, DIMENSION (klon,6) :: alb6
+    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, qsalt
 ! End definition
 !****************************************************************************************
@@ -165 +174 @@
     alb2(:) = 999999.
     alb1(:) = 999999.
-    
+    fluxbs(:)=0.  
     runoff(:) = 0.
 !****************************************************************************************
@@ -173 +182 @@
     radsol(:) = 0.0
     radsol(1:knon) = swnet(1:knon) + lwnet(1:knon)
+
+!****************************************************************************************
 
 !****************************************************************************************
@@ -264 +275 @@
 
 
-
     ELSE 
 
@@ -314 +324 @@
          flux_u1, flux_v1)
 
-!****************************************************************************************
-! Calculate snow height, age, run-off,..
-!    
-!****************************************************************************************
-    CALL fonte_neige(knon, is_lic, knindex, dtime, &
-         tsurf, precip_rain, precip_snow, &
-         snow, qsol, tsurf_new, evap)
-
 
 !****************************************************************************************
@@ -327 +329 @@
 !
 !****************************************************************************************
+ 
     CALL albsno(klon,knon,dtime,agesno(:),alb_neig(:), precip_snow(:))  
 
-    WHERE (snow(1 : knon) .LT. 0.0001) agesno(1 : knon) = 0.
-    zfra(1:knon) = MAX(0.0,MIN(1.0,snow(1:knon)/(snow(1:knon)+10.0)))
-    alb1(1:knon) = alb_neig(1:knon)*zfra(1:knon) + &
-         0.6 * (1.0-zfra(1:knon))
+
+! 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"
@@ -354 +358 @@
     !z0m = SQRT(z0m**2+rugoro**2)
 
+
+
+  ! Simple blowing snow param
+  if (ok_bs) then
+       ustart0 = 0.211
+       rhoice = 920.0
+       rho0 = 200.0
+       rhomax=450.0
+       rhohard=400.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
+
+       ! computation of threshold friction velocity
+       ! which depends on surface snow density
+       do i = 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))
+           ! 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)) 
+           ! we have multiplied by exp to prevent erosion when rhos>rhomax (usefull till
+           ! rhohard<450)
+       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*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)
+       enddo
+
+  endif
+
+
+
+!****************************************************************************************
+! Calculate surface snow amount
+!    
+!****************************************************************************************
+    IF (ok_bs) THEN
+      precip_totsnow(:)=precip_snow(:)+precip_bs(:)
+      evap_totsnow(:)=evap(:)-fluxbs(:) ! flux bs is positive towards the surface (snow erosion)
+    ELSE
+      precip_totsnow(:)=precip_snow(:)
+      evap_totsnow(:)=evap(:)
+    ENDIF
+
+    CALL fonte_neige(knon, is_lic, knindex, dtime, &
+         tsurf, precip_rain, precip_totsnow,  &
+         snow, qsol, tsurf_new, evap_totsnow)
+
+    WHERE (snow(1 : knon) .LT. 0.0001) agesno(1 : knon) = 0.                                         
+    zfra(1:knon) = MAX(0.0,MIN(1.0,snow(1:knon)/(snow(1:knon)+10.0)))  
 
 
@@ -378 +482 @@
     runoff(1:knon)=run_off_lic(1:knon)/dtime
 
-  
-!****************************************************************************************
        snow_o=0.
        zfra_o = 0.
@@ -420 +522 @@
 !albedo SB <<<
 
-  
-  
 
   END SUBROUTINE surf_landice