diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/Process_Library.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/Process_Library.F90
index 9ec85d331..220b3009c 100644
--- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/Process_Library.F90
+++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/Process_Library.F90
@@ -32,6 +32,13 @@ module GEOSmoist_Process_Library
     module procedure ICE_FRACTION_1D
     module procedure ICE_FRACTION_SC
   end interface ICE_FRACTION
+
+  ! SRF_TYPE constants
+  integer, parameter :: SRF_TYPE_LAND  = 1
+  integer, parameter :: SRF_TYPE_SNOW  = 2
+  integer, parameter :: SRF_TYPE_ICE   = 3
+  integer, parameter :: SRF_TYPE_OCEAN = 0
+
   ! ICE_FRACTION constants
    ! In anvil/convective clouds
    real, parameter :: aT_ICE_ALL = 252.16
@@ -403,7 +410,8 @@ function ICE_FRACTION_SC (TEMP,CNV_FRACTION,SRF_TYPE) RESULT(ICEFRCT)
       ICEFRCT_C = MAX(ICEFRCT_C,0.00)
       ICEFRCT_C = ICEFRCT_C**aICEFRPWR
      ! Sigmoidal functions like figure 6b/6c of Hu et al 2010, doi:10.1029/2009JD012384
-      if (SRF_TYPE >= 2.0) then
+      select case (nint(SRF_TYPE))
+      case (SRF_TYPE_SNOW, SRF_TYPE_ICE)
         ! Over snow (SRF_TYPE == 2.0) and ice (SRF_TYPE == 3.0)
         if (ICE_RADII_PARAM == 1) then
           ! Jason formula
@@ -424,8 +432,8 @@ function ICE_FRACTION_SC (TEMP,CNV_FRACTION,SRF_TYPE) RESULT(ICEFRCT)
         ICEFRCT_M = MIN(ICEFRCT_M,1.00)
         ICEFRCT_M = MAX(ICEFRCT_M,0.00)
         ICEFRCT_M = ICEFRCT_M**iICEFRPWR
-      else if (SRF_TYPE > 1.0) then
-        ! Over Land
+      case (SRF_TYPE_LAND)
+        ! Over Land (SRF_TYPE == 1)
         ICEFRCT_M  = 0.00
         if ( TEMP <= lT_ICE_ALL ) then
            ICEFRCT_M = 1.000
@@ -435,8 +443,8 @@ function ICE_FRACTION_SC (TEMP,CNV_FRACTION,SRF_TYPE) RESULT(ICEFRCT)
         ICEFRCT_M = MIN(ICEFRCT_M,1.00)
         ICEFRCT_M = MAX(ICEFRCT_M,0.00)
         ICEFRCT_M = ICEFRCT_M**lICEFRPWR
-      else
-        ! Over Oceans
+      case (SRF_TYPE_OCEAN)
+        ! Over Oceans (SRF_TYPE == 0)
         ICEFRCT_M  = 0.00
         if ( TEMP <= oT_ICE_ALL ) then
            ICEFRCT_M = 1.000
@@ -446,7 +454,11 @@ function ICE_FRACTION_SC (TEMP,CNV_FRACTION,SRF_TYPE) RESULT(ICEFRCT)
         ICEFRCT_M = MIN(ICEFRCT_M,1.00)
         ICEFRCT_M = MAX(ICEFRCT_M,0.00)
         ICEFRCT_M = ICEFRCT_M**oICEFRPWR
-      endif
+      case default
+        ! You should not be here
+        print *, 'ICE_FRACTION_SC: Unknown SRF_TYPE = ',SRF_TYPE
+        error stop
+      end select
       ! Combine the Convective and MODIS functions
         ICEFRCT  = ICEFRCT_M*(1.0-CNV_FRACTION) + ICEFRCT_C*(CNV_FRACTION)
 #endif