use StieglitzSnow_calc_asnow() instead of writing out equation each t…

…ime (StieglitzSnow.F90)
GEOS-ESM · Oct 15, 2023 · 7856ac7 · 7856ac7
1 parent b036f48
commit 7856ac7
Show file tree

Hide file tree

Showing 2 changed files with 13 additions and 15 deletions.
diff --git a/...dComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/Shared/lsm_routines.F90 b/...dComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/Shared/lsm_routines.F90
@@ -1298,9 +1298,6 @@ SUBROUTINE SIBALB (NCH, ITYP, VLAI, VGRN, ZTH,                 &
       REAL, PARAMETER :: ALVDRI = 0.700
       REAL, PARAMETER :: ALIDRI = 0.700
 
-
-!      REAL, PARAMETER :: WEMIN  = 13.0   ! [KG/M2]
-
 ! ALVDRS:  Albedo of soil for visible   direct  solar radiation.
 ! ALIDRS:  Albedo of soil for infra-red direct  solar radiation.
 ! ALVDFS:  Albedo of soil for visible   diffuse solar radiation.

diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Shared/StieglitzSnow.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Shared/StieglitzSnow.F90
@@ -195,8 +195,8 @@ subroutine StieglitzSnow_snowrt(N_zones, N_snow, tileType,                     &
        hlwout, fices, tpsn, rmelt,                                               &  ! out
        areasc, areasc0, pre, fhgnd, evap, shflux, lhflux, hcorr, ghfluxsno,      &  ! out
        sndzsc, wesnprec, sndzprec, sndz1perc,                                    &  ! out
-       wesnperc, wesndens, wesnrepar, mltwtr,                                    &  ! out
-       excs, drho0, wesnbot, tksno, dtss        )
+       wesnperc, wesndens, wesnrepar, mltwtr,                                    &  ! out  
+       excs, drho0, wesnbot, tksno, dtss        )                                   ! out
 
     !*********************************************************************
     ! AUTHORS:  M. Stieglitz, M. Suarez, R. Koster & S. Dery.
@@ -404,7 +404,7 @@ subroutine StieglitzSnow_snowrt(N_zones, N_snow, tileType,                     &
        ! Melt off initial (very small) snowpack; new snow pack is based
        !   on new snowfall only (if any)
 
-       areasc = min(snowd/wemin,1.)
+       call StieglitzSnow_calc_asnow( snowd, areasc )
        areasc0 = 0.
        pre = snowd/dts + areasc*rainf       ! pre = melted snowpack plus rainfall
        wesn  = 0.
@@ -422,7 +422,7 @@ subroutine StieglitzSnow_snowrt(N_zones, N_snow, tileType,                     &
 
           wesn    = snowf*dts/float(N_snow)  
           htsnn   = (tsx-alhm)*wesn
-          areasc0 = min((snowf*dts)/wemin,1.)
+          call StieglitzSnow_calc_asnow( snowf*dts, areasc0 )
 
           !*** should have fractional snow cover taken into account
           sndz = wesn/(max(areasc0,small)*rhofs)
@@ -453,7 +453,7 @@ subroutine StieglitzSnow_snowrt(N_zones, N_snow, tileType,                     &
 
     !**** Determine the fractional snow coverage
 
-    areasc = min(snowd/wemin,1.)
+    call StieglitzSnow_calc_asnow( snowd, areasc )
 
     !**** Set the mean density & diffusivity of the layers
 
@@ -748,9 +748,9 @@ subroutine StieglitzSnow_snowrt(N_zones, N_snow, tileType,                     &
 
        if(snowd > wemin) then
 
-          icedens=wesn(i)*fices(i)/(sndz(i)+1.e-20)
-          densfac=amax1(0., amin1(1., icedens/rhofs))
-          term=densfac*snfr*(sndz(i)*rhow-wesn(i)*fices(i))
+          icedens = wesn(i)*fices(i)/(sndz(i)+1.e-20)
+          densfac = amax1(0., amin1(1., icedens/rhofs))
+          term    = densfac*snfr*(sndz(i)*rhow-wesn(i)*fices(i))
 
           if(pre > term) then
              pre = min(pre - term, wesn(i))        ! when asnow=1, retain some liquid water in snow pack
@@ -881,7 +881,8 @@ subroutine StieglitzSnow_snowrt(N_zones, N_snow, tileType,                     &
     excs = excs * fices / dts
 
     snowd=sum(wesn)
-    areasc0 = max(small, min(snowd/wemin,1.) )
+    call StieglitzSnow_calc_asnow( snowd, areasc0 )
+    areasc0 = max(small, areasc0 )
     sndz = (wesn/areasc0)/dens
 
     sndzsum = sum(sndz)
@@ -914,7 +915,7 @@ subroutine StieglitzSnow_snowrt(N_zones, N_snow, tileType,                     &
 
     !**** Reset fractional area coverage.
 
-    areasc0 = min(sum(wesn)/wemin,1.)
+    call StieglitzSnow_calc_asnow( sum(wesn), areasc0 )
 
     !**** Final check for water balance.
 
@@ -1653,7 +1654,7 @@ SUBROUTINE StieglitzSnow_snow_albedo(                                  &
        SWE=SUM(WESN(:,I))
 
        TOTDEP=SNDZ(1,I)
-       AREASC = MIN(SWE/WEMIN,1.)
+       call StieglitzSnow_calc_asnow( SWE, AREASC )
        !DENSITY=(SWE/(AREASC+1.e-20)) / (TOTDEP+1.e-20)
        !*** only use top layer density to dentermine albedo 
        DENSITY=(WESN(1,I)/(AREASC+1.e-20)) / (TOTDEP+1.e-20)
@@ -1803,7 +1804,7 @@ SUBROUTINE ALB_WITH_IMPURITY (N_snow, ZTH,               &
 
     SWE=SUM(WESN(:))
     TOTDEP=SUM(SNDZ(:))
-    AREASC = MIN(SWE/WEMIN,1.)
+    call StieglitzSnow_calc_asnow( SWE, AREASC )
     DENSITY=(SWE/(AREASC+1.e-20)) / (TOTDEP+1.e-20)
 
     WSS=UM