diff --git a/modulefiles/build.hera.gnu.lua b/modulefiles/build.hera.gnu.lua
index 3756fd3e0..daeeb9cbc 100644
--- a/modulefiles/build.hera.gnu.lua
+++ b/modulefiles/build.hera.gnu.lua
@@ -31,7 +31,7 @@ load(pathJoin("bacio", bacio_ver))
 g2_ver=os.getenv("g2_ver") or "3.4.3"
 load(pathJoin("g2", g2_ver))
 
-ip_ver=os.getenv("ip_ver") or "3.3.3"
+ip_ver=os.getenv("ip_ver") or "4.0.0"
 load(pathJoin("ip", ip_ver))
 
 nemsio_ver=os.getenv("nemsio_ver") or "2.5.2"
diff --git a/modulefiles/build.hera.intel.lua b/modulefiles/build.hera.intel.lua
index 6637da08a..a3c823ad6 100644
--- a/modulefiles/build.hera.intel.lua
+++ b/modulefiles/build.hera.intel.lua
@@ -25,7 +25,7 @@ load(pathJoin("bacio", bacio_ver))
 g2_ver=os.getenv("g2_ver") or "3.4.5"
 load(pathJoin("g2", g2_ver))
 
-ip_ver=os.getenv("ip_ver") or "3.3.3"
+ip_ver=os.getenv("ip_ver") or "4.0.0"
 load(pathJoin("ip", ip_ver))
 
 nemsio_ver=os.getenv("nemsio_ver") or "2.5.4"
diff --git a/modulefiles/build.jet.intel.lua b/modulefiles/build.jet.intel.lua
index c85f8e498..d6a0a877b 100644
--- a/modulefiles/build.jet.intel.lua
+++ b/modulefiles/build.jet.intel.lua
@@ -37,7 +37,7 @@ load(pathJoin("w3nco", w3nco_ver))
 sp_ver=os.getenv("sp_ver") or "2.3.3"
 load(pathJoin("sp", sp_ver))
 
-ip_ver=os.getenv("ip_ver") or "3.3.3"
+ip_ver=os.getenv("ip_ver") or "4.0.0"
 load(pathJoin("ip", ip_ver))
 
 bacio_ver=os.getenv("bacio_ver") or "2.4.1"
diff --git a/modulefiles/build.orion.intel.lua b/modulefiles/build.orion.intel.lua
index b2efc0045..553069024 100644
--- a/modulefiles/build.orion.intel.lua
+++ b/modulefiles/build.orion.intel.lua
@@ -5,7 +5,7 @@ Load environment to compile UFS_UTILS on Orion
 cmake_ver=os.getenv("cmake_ver") or "3.17.3"
 load(pathJoin("cmake", cmake_ver))
 
-prepend_path("MODULEPATH", "/apps/contrib/NCEP/libs/hpc-stack/modulefiles/stack")
+prepend_path("MODULEPATH", "/apps/contrib/NCEP/hpc-stack/libs/hpc-stack/modulefiles/stack")
 
 hpc_ver=os.getenv("hpc_ver") or "1.2.0"
 load(pathJoin("hpc", hpc_ver))
@@ -22,7 +22,7 @@ load(pathJoin("bacio", bacio_ver))
 g2_ver=os.getenv("g2_ver") or "3.4.5"
 load(pathJoin("g2", g2_ver))
 
-ip_ver=os.getenv("ip_ver") or "3.3.3"
+ip_ver=os.getenv("ip_ver") or "4.0.0"
 load(pathJoin("ip", ip_ver))
 
 nemsio_ver=os.getenv("nemsio_ver") or "2.5.4"
diff --git a/modulefiles/build.wcoss2.intel.lua b/modulefiles/build.wcoss2.intel.lua
index 7fa476848..1c8953c56 100644
--- a/modulefiles/build.wcoss2.intel.lua
+++ b/modulefiles/build.wcoss2.intel.lua
@@ -17,7 +17,6 @@ load(pathJoin("intel", intel_ver))
 cray_mpich_ver=os.getenv("cray_mpich_ver") or "8.1.7"
 load(pathJoin("cray-mpich", cray_mpich_ver))
 
-
 libjpeg_ver=os.getenv("libjpeg_ver") or "9c"
 load(pathJoin("libjpeg", libjpeg_ver))
 
@@ -51,7 +50,7 @@ load(pathJoin("sigio", sigio_ver))
 sp_ver=os.getenv("sp_ver") or "2.3.3"
 load(pathJoin("sp", sp_ver))
 
-ip_ver=os.getenv("ip_ver") or "3.3.3"
+ip_ver=os.getenv("ip_ver") or "4.0.0"
 load(pathJoin("ip", ip_ver))
 
 g2_ver=os.getenv("g2_ver") or "3.4.5"
diff --git a/sorc/emcsfc_snow2mdl.fd/CMakeLists.txt b/sorc/emcsfc_snow2mdl.fd/CMakeLists.txt
index b20a722ee..057422d83 100644
--- a/sorc/emcsfc_snow2mdl.fd/CMakeLists.txt
+++ b/sorc/emcsfc_snow2mdl.fd/CMakeLists.txt
@@ -17,6 +17,9 @@ if(CMAKE_Fortran_COMPILER_ID MATCHES "^(Intel)$")
 elseif(CMAKE_Fortran_COMPILER_ID MATCHES "^(GNU)$")
   set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -fdefault-real-8 -fconvert=big-endian")
 endif()
+if(ip_VERSION GREATER_EQUAL 4.0.0)
+  set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -DIP_V4")
+endif()
 
 set(exe_name emcsfc_snow2mdl)
 
diff --git a/sorc/emcsfc_snow2mdl.fd/snow2mdl.F90 b/sorc/emcsfc_snow2mdl.fd/snow2mdl.F90
index 967117b75..802ca27b5 100755
--- a/sorc/emcsfc_snow2mdl.fd/snow2mdl.F90
+++ b/sorc/emcsfc_snow2mdl.fd/snow2mdl.F90
@@ -157,22 +157,29 @@ module snow2mdl
  subroutine interp
  use gdswzd_mod
 
+! Required when using iplib v4.0 or higher.
+#ifdef IP_V4
+ use ipolates_mod
+#endif
+
  implicit none
 
+ integer, parameter        :: km=1
+
  integer                   :: i, j, ii, jj, ij
  integer                   :: ijmdl2, istart, iend, imid, iii
  integer, allocatable      :: idum(:,:)
- integer                   :: int_opt, ipopt(20)
+ integer                   :: int_opt, ipopt(20), ibi(km)
  integer                   :: kgds_mdl_tmp(200)
- integer                   :: no, ibo, iret, nret
+ integer                   :: no, ibo(km), iret, nret
 
- logical*1, allocatable    :: bitmap_mdl(:)
+ logical*1, allocatable    :: bitmap_mdl(:,:)
 
- real                      :: gridi(1)
- real                      :: gridj(1)
+ real                      :: gridi(1), gridj(1)
+ real                      :: lats(1), lons(1)
  real, allocatable         :: lsmask_1d(:)
- real, allocatable         :: snow_cvr_mdl_1d(:)
- real, allocatable         :: snow_dep_mdl_tmp(:) 
+ real, allocatable         :: snow_cvr_mdl_1d(:,:)
+ real, allocatable         :: snow_dep_mdl_tmp(:,:) 
  real                      :: sumc, sumd, x1, r, fraction, gridis, gridie
  real, parameter           :: undefined_value = -999.
 
@@ -229,6 +236,7 @@ subroutine interp
  NESDIS_IMS : if (use_nesdis) then
 
    ipopt = 0
+   ibi = 1
    if (nesdis_res < (0.5*resol_mdl)) then
      print*,"- INTERPOLATE NH NESDIS/IMS DATA TO MODEL GRID USING BUDGET METHOD."
      ipopt(1)=2  ! break model grid cell into 25 points.
@@ -248,16 +256,17 @@ subroutine interp
      no = ijmdl  ! an input when kgds(1) < 0 (subset of grid)
    end if
 
-   allocate (snow_cvr_mdl_1d(ijmdl))
+   allocate (snow_cvr_mdl_1d(ijmdl,km))
    snow_cvr_mdl_1d = 0.0
 
-   allocate (bitmap_mdl(ijmdl))
+   allocate (bitmap_mdl(ijmdl,km))
    bitmap_mdl=.false.  ! if interpolation routine can't find data
                        ! at a point, this flag is false.
 
+ 
    call ipolates(int_opt, ipopt, kgds_nesdis, kgds_mdl_tmp,   &
                 (inesdis*jnesdis), ijmdl,               &
-                 1, 1, bitmap_nesdis, snow_cvr_nesdis,  &
+                 km, ibi, bitmap_nesdis, snow_cvr_nesdis,  &
                  no, lats_mdl, lons_mdl, ibo, bitmap_mdl,     &
                  snow_cvr_mdl_1d, iret)
 
@@ -279,22 +288,24 @@ subroutine interp
 
    do ij = 1, ijmdl
      if (lats_mdl(ij) < 0.0) cycle  ! only consider nh model points
-     if (.not. bitmap_mdl(ij)) then 
+     if (.not. bitmap_mdl(ij,km)) then 
        if (lats_mdl(ij) <= lat_threshold) then
-         snow_cvr_mdl_1d(ij) = 0.0
+         snow_cvr_mdl_1d(ij,km) = 0.0
        else 
+         lats(1)=lats_mdl(ij)
+         lons(1)=lons_mdl(ij)
          call gdswzd(kgds_nesdis,-1,1,undefined_value,gridi,gridj, &
-                     lons_mdl(ij),lats_mdl(ij),nret) 
+                     lons,lats,nret) 
          if (nret /= 1) then
            print*,"- WARNING: MODEL POINT OUTSIDE NESDIS/IMS GRID: ", ipts_mdl(ij), jpts_mdl(ij)
-           snow_cvr_mdl_1d(ij) = 0.0
+           snow_cvr_mdl_1d(ij,km) = 0.0
          else
            ii = nint(gridi(1))
            jj = nint(gridj(1))
            if (sea_ice_nesdis(ii,jj) == 1) then
-             snow_cvr_mdl_1d(ij) = 100.0
+             snow_cvr_mdl_1d(ij,km) = 100.0
            else
-             snow_cvr_mdl_1d(ij) = 0.0
+             snow_cvr_mdl_1d(ij,km) = 0.0
            end if
          end if
        end if
@@ -318,6 +329,7 @@ subroutine interp
 !----------------------------------------------------------------------
 
    ipopt = 0
+   ibi = 1
    if (afwa_res < (0.5*resol_mdl)) then
      print*,"- INTERPOLATE GLOBAL AFWA DATA TO MODEL GRID USING BUDGET METHOD."
      ipopt(1)=-1  ! break model grid cell into 25 points.
@@ -336,15 +348,15 @@ subroutine interp
      no = ijmdl  ! an input when kgds(1) < 0 (subset of grid)
    end if
 
-   allocate (snow_dep_mdl_tmp(ijmdl))
+   allocate (snow_dep_mdl_tmp(ijmdl,km))
    snow_dep_mdl_tmp = 0.0
 
-   allocate (bitmap_mdl(ijmdl))
+   allocate (bitmap_mdl(ijmdl,km))
    bitmap_mdl = .false.
 
    call ipolates(int_opt, ipopt, kgds_afwa_global, kgds_mdl_tmp,    &
                 (iafwa*jafwa), ijmdl,  &
-                 1, 1, bitmap_afwa_global, snow_dep_afwa_global, &
+                 km, ibi, bitmap_afwa_global, snow_dep_afwa_global, &
                  no, lats_mdl, lons_mdl, ibo, bitmap_mdl,     &
                  snow_dep_mdl_tmp, iret)
 
@@ -364,11 +376,11 @@ subroutine interp
 !----------------------------------------------------------------------
 
    do ij = 1, ijmdl
-     if (.not. bitmap_mdl(ij)) then
+     if (.not. bitmap_mdl(ij,km)) then
        if (abs(lats_mdl(ij)) >= lat_threshold) then
-         snow_dep_mdl_tmp(ij) = min_snow_depth
+         snow_dep_mdl_tmp(ij,km) = min_snow_depth
        else
-         snow_dep_mdl_tmp(ij) = 0.0
+         snow_dep_mdl_tmp(ij,km) = 0.0
        endif
      endif
    enddo
@@ -389,6 +401,7 @@ subroutine interp
 !----------------------------------------------------------------------
 
    ipopt = 0
+   ibi = 1
    if (afwa_res < (0.5*resol_mdl)) then
      print*,"- INTERPOLATE NH AFWA DATA TO MODEL GRID USING BUDGET METHOD."
      ipopt(1)=-1  ! break model grid cell into 25 points.
@@ -407,15 +420,15 @@ subroutine interp
      no = ijmdl  ! an input when kgds(1) < 0 (subset of grid)
    end if
 
-   allocate (snow_dep_mdl_tmp(ijmdl))
+   allocate (snow_dep_mdl_tmp(ijmdl,km))
    snow_dep_mdl_tmp = 0.0
 
-   allocate (bitmap_mdl(ijmdl))
+   allocate (bitmap_mdl(ijmdl,km))
    bitmap_mdl = .false.
 
    call ipolates(int_opt, ipopt, kgds_afwa_nh, kgds_mdl_tmp,    &
                 (iafwa*jafwa), ijmdl,  &
-                 1, 1, bitmap_afwa_nh, snow_dep_afwa_nh, &
+                 1, ibi, bitmap_afwa_nh, snow_dep_afwa_nh, &
                  no, lats_mdl, lons_mdl, ibo, bitmap_mdl,     &
                  snow_dep_mdl_tmp, iret)
 
@@ -436,11 +449,11 @@ subroutine interp
 
    do ij = 1, ijmdl
      if (lats_mdl(ij) >= 0.) then  ! only consider model pts in n hemi.
-       if (.not. bitmap_mdl(ij)) then
+       if (.not. bitmap_mdl(ij,km)) then
          if (abs(lats_mdl(ij)) >= lat_threshold) then
-           snow_dep_mdl_tmp(ij) = min_snow_depth
+           snow_dep_mdl_tmp(ij,km) = min_snow_depth
          else
-           snow_dep_mdl_tmp(ij) = 0.0
+           snow_dep_mdl_tmp(ij,km) = 0.0
          endif
        endif
      endif
@@ -464,11 +477,11 @@ subroutine interp
    print*,"- BLEND NESDIS/IMS AND AFWA DATA IN NH."
    do ij = 1, ijmdl
      if (lats_mdl(ij) >= 0.0) then
-       if (snow_cvr_mdl_1d(ij) >= snow_cvr_threshold) then
+       if (snow_cvr_mdl_1d(ij,km) >= snow_cvr_threshold) then
          snow_dep_mdl(ipts_mdl(ij),jpts_mdl(ij)) =   &
-                      max(snow_dep_mdl_tmp(ij), min_snow_depth)
+                      max(snow_dep_mdl_tmp(ij,km), min_snow_depth)
        endif   
-       snow_cvr_mdl(ipts_mdl(ij),jpts_mdl(ij)) = snow_cvr_mdl_1d(ij)
+       snow_cvr_mdl(ipts_mdl(ij),jpts_mdl(ij)) = snow_cvr_mdl_1d(ij,km)
      endif
    enddo
    deallocate (snow_cvr_mdl_1d)
@@ -476,11 +489,11 @@ subroutine interp
    print*,"- BLEND NESDIS/IMS AND AFWA DATA IN NH."
    do ij = 1, ijmdl
      if (lats_mdl(ij) >= 0.0) then
-       if (snow_cvr_mdl_1d(ij) >= snow_cvr_threshold) then
+       if (snow_cvr_mdl_1d(ij,km) >= snow_cvr_threshold) then
          snow_dep_mdl(ipts_mdl(ij),jpts_mdl(ij)) =   &
-                      max(snow_dep_mdl_tmp(ij), min_snow_depth)
+                      max(snow_dep_mdl_tmp(ij,km), min_snow_depth)
        endif   
-       snow_cvr_mdl(ipts_mdl(ij),jpts_mdl(ij)) = snow_cvr_mdl_1d(ij)
+       snow_cvr_mdl(ipts_mdl(ij),jpts_mdl(ij)) = snow_cvr_mdl_1d(ij,km)
      endif
    enddo
    deallocate (snow_cvr_mdl_1d)
@@ -489,9 +502,9 @@ subroutine interp
    print*,"- SET DEPTH/COVER FROM AFWA DATA IN NH."
    do ij = 1, ijmdl
      if (lats_mdl(ij) >= 0.0) then
-       if (snow_dep_mdl_tmp(ij) > 0.0) then
+       if (snow_dep_mdl_tmp(ij,km) > 0.0) then
          snow_cvr_mdl(ipts_mdl(ij),jpts_mdl(ij)) = 100.0
-         snow_dep_mdl(ipts_mdl(ij),jpts_mdl(ij)) = snow_dep_mdl_tmp(ij)
+         snow_dep_mdl(ipts_mdl(ij),jpts_mdl(ij)) = snow_dep_mdl_tmp(ij,km)
        endif
      endif
    enddo
@@ -499,9 +512,9 @@ subroutine interp
    print*,"- SET DEPTH/COVER FROM AFWA DATA IN NH."
    do ij = 1, ijmdl
      if (lats_mdl(ij) >= 0.0) then
-       if (snow_dep_mdl_tmp(ij) > 0.0) then
+       if (snow_dep_mdl_tmp(ij,km) > 0.0) then
          snow_cvr_mdl(ipts_mdl(ij),jpts_mdl(ij)) = 100.0
-         snow_dep_mdl(ipts_mdl(ij),jpts_mdl(ij)) = snow_dep_mdl_tmp(ij)
+         snow_dep_mdl(ipts_mdl(ij),jpts_mdl(ij)) = snow_dep_mdl_tmp(ij,km)
        endif
      endif
    enddo
@@ -510,10 +523,10 @@ subroutine interp
    print*,"- SET DEPTH/COVER FROM NESDIS/IMS DATA IN NH."
    do ij = 1, ijmdl
      if (lats_mdl(ij) >= 0.0) then
-       if (snow_cvr_mdl_1d(ij) >= snow_cvr_threshold) then
+       if (snow_cvr_mdl_1d(ij,km) >= snow_cvr_threshold) then
          snow_dep_mdl(ipts_mdl(ij),jpts_mdl(ij)) = min_snow_depth
        endif   
-       snow_cvr_mdl(ipts_mdl(ij),jpts_mdl(ij)) = snow_cvr_mdl_1d(ij)
+       snow_cvr_mdl(ipts_mdl(ij),jpts_mdl(ij)) = snow_cvr_mdl_1d(ij,km)
      endif
    enddo
    deallocate (snow_cvr_mdl_1d)
@@ -526,6 +539,7 @@ subroutine interp
  AUTOSNOW : if (use_autosnow) then
 
    ipopt = 0
+   ibi = 1
    if (autosnow_res < (0.5*resol_mdl)) then
      print*,"- INTERPOLATE AUTOSNOW DATA TO MODEL GRID USING BUDGET METHOD."
      ipopt(1)=2    ! break model grid cell into 25 points.
@@ -545,16 +559,16 @@ subroutine interp
      no = ijmdl  ! an input when kgds(1) < 0 (subset of grid)
    end if
 
-   allocate (snow_cvr_mdl_1d(ijmdl))
+   allocate (snow_cvr_mdl_1d(ijmdl,km))
    snow_cvr_mdl_1d = 0.0
 
-   allocate (bitmap_mdl(ijmdl))
+   allocate (bitmap_mdl(ijmdl,km))
    bitmap_mdl=.false.  ! if interpolation routine can't find data
                        ! at a point, this flag is false.
 
    call ipolates(int_opt, ipopt, kgds_autosnow, kgds_mdl_tmp,   &
                 (iautosnow*jautosnow), ijmdl,               &
-                 1, 1, bitmap_autosnow, snow_cvr_autosnow,  &
+                 1, ibi, bitmap_autosnow, snow_cvr_autosnow,  &
                  no, lats_mdl, lons_mdl, ibo, bitmap_mdl,     &
                  snow_cvr_mdl_1d, iret)
 
@@ -573,11 +587,11 @@ subroutine interp
 
    do ij = 1, ijmdl
      if (lats_mdl(ij) < 0.0) then
-       if (.not. bitmap_mdl(ij)) then
+       if (.not. bitmap_mdl(ij,km)) then
          if (abs(lats_mdl(ij)) <= lat_threshold) then
-           snow_cvr_mdl_1d(ij) = 0.0
+           snow_cvr_mdl_1d(ij,km) = 0.0
          else 
-           snow_cvr_mdl_1d(ij) = 100.0
+           snow_cvr_mdl_1d(ij,km) = 100.0
          end if
        end if
      end if
@@ -599,6 +613,7 @@ subroutine interp
 !----------------------------------------------------------------------
 
    ipopt = 0
+   ibi = 1
    if (afwa_res < (0.5*resol_mdl)) then
      print*,"- INTERPOLATE SH AFWA DATA TO MODEL GRID USING BUDGET METHOD."
      ipopt(1)=-1  ! break model grid cell into 25 points.
@@ -617,15 +632,15 @@ subroutine interp
      no = ijmdl  ! an input when kgds(1) < 0 (subset of grid)
    end if
 
-   allocate (snow_dep_mdl_tmp(ijmdl))
+   allocate (snow_dep_mdl_tmp(ijmdl,km))
    snow_dep_mdl_tmp = 0.0
 
-   allocate (bitmap_mdl(ijmdl))
+   allocate (bitmap_mdl(ijmdl,km))
    bitmap_mdl = .false.
 
    call ipolates(int_opt, ipopt, kgds_afwa_sh, kgds_mdl_tmp,    &
                 (iafwa*jafwa), ijmdl,  &
-                 1, 1, bitmap_afwa_sh, snow_dep_afwa_sh, &
+                 1, ibi, bitmap_afwa_sh, snow_dep_afwa_sh, &
                  no, lats_mdl, lons_mdl, ibo, bitmap_mdl,     &
                  snow_dep_mdl_tmp, iret)
 
@@ -644,11 +659,11 @@ subroutine interp
 
    do ij = 1, ijmdl
      if (lats_mdl(ij) < 0.) then
-       if (.not. bitmap_mdl(ij)) then
+       if (.not. bitmap_mdl(ij,km)) then
          if (abs(lats_mdl(ij)) >= lat_threshold) then
-           snow_dep_mdl_tmp(ij) = min_snow_depth
+           snow_dep_mdl_tmp(ij,km) = min_snow_depth
          else
-           snow_dep_mdl_tmp(ij) = 0.0
+           snow_dep_mdl_tmp(ij,km) = 0.0
          endif
        endif
      endif
@@ -667,11 +682,11 @@ subroutine interp
    print*,"- BLEND AUTOSNOW AND AFWA DATA IN SH."
    do ij = 1, ijmdl
      if (lats_mdl(ij) < 0.0) then
-       if (snow_cvr_mdl_1d(ij) >= snow_cvr_threshold) then
+       if (snow_cvr_mdl_1d(ij,km) >= snow_cvr_threshold) then
          snow_dep_mdl(ipts_mdl(ij),jpts_mdl(ij)) =   &
-                      max(snow_dep_mdl_tmp(ij), min_snow_depth)
+                      max(snow_dep_mdl_tmp(ij,km), min_snow_depth)
        endif   
-       snow_cvr_mdl(ipts_mdl(ij),jpts_mdl(ij)) = snow_cvr_mdl_1d(ij)
+       snow_cvr_mdl(ipts_mdl(ij),jpts_mdl(ij)) = snow_cvr_mdl_1d(ij,km)
      endif
    enddo
    deallocate (snow_cvr_mdl_1d)
@@ -680,9 +695,9 @@ subroutine interp
    print*,"- SET DEPTH/COVER FROM AFWA DATA IN SH."
    do ij = 1, ijmdl
      if (lats_mdl(ij) < 0.0) then
-       if (snow_dep_mdl_tmp(ij) > 0.0) then
+       if (snow_dep_mdl_tmp(ij,km) > 0.0) then
          snow_cvr_mdl(ipts_mdl(ij),jpts_mdl(ij)) = 100.0
-         snow_dep_mdl(ipts_mdl(ij),jpts_mdl(ij)) = snow_dep_mdl_tmp(ij)
+         snow_dep_mdl(ipts_mdl(ij),jpts_mdl(ij)) = snow_dep_mdl_tmp(ij,km)
        endif
      endif
    enddo
@@ -691,10 +706,10 @@ subroutine interp
    print*,"- SET DEPTH/COVER FROM AUTOSNOW IN SH."
    do ij = 1, ijmdl
      if (lats_mdl(ij) < 0.0) then
-       if (snow_cvr_mdl_1d(ij) >= snow_cvr_threshold) then
+       if (snow_cvr_mdl_1d(ij,km) >= snow_cvr_threshold) then
          snow_dep_mdl(ipts_mdl(ij),jpts_mdl(ij)) = min_snow_depth
        endif   
-       snow_cvr_mdl(ipts_mdl(ij),jpts_mdl(ij)) = snow_cvr_mdl_1d(ij)
+       snow_cvr_mdl(ipts_mdl(ij),jpts_mdl(ij)) = snow_cvr_mdl_1d(ij,km)
      endif
    enddo
    deallocate (snow_cvr_mdl_1d)
@@ -709,9 +724,9 @@ subroutine interp
  if (kgds_mdl(1) == 4 .and. thinned) then
 
    ijmdl2 = sum(lonsperlat_mdl) * 2
-   allocate (snow_cvr_mdl_1d(ijmdl2))
+   allocate (snow_cvr_mdl_1d(ijmdl2,km))
    allocate (lsmask_1d(ijmdl2))
-   allocate (snow_dep_mdl_tmp(ijmdl2))
+   allocate (snow_dep_mdl_tmp(ijmdl2,km))
 
    lsmask_1d = 0.0
    snow_cvr_mdl_1d = 0.0
@@ -748,10 +763,10 @@ subroutine interp
          sumc = sumc + fraction * snow_cvr_mdl(iii,j)
          sumd = sumd + fraction * snow_dep_mdl(iii,j)
        enddo
-       snow_cvr_mdl_1d(ij) = sumc / r
-       snow_dep_mdl_tmp(ij) = 0.0
-       if (snow_cvr_mdl_1d(ij) > snow_cvr_threshold) then
-         snow_dep_mdl_tmp(ij) = max(sumd / r,min_snow_depth)
+       snow_cvr_mdl_1d(ij,km) = sumc / r
+       snow_dep_mdl_tmp(ij,km) = 0.0
+       if (snow_cvr_mdl_1d(ij,km) > snow_cvr_threshold) then
+         snow_dep_mdl_tmp(ij,km) = max(sumd / r,min_snow_depth)
        end if
     enddo
    enddo
@@ -765,9 +780,9 @@ subroutine interp
    allocate (idum(imdl,jmdl))
    idum = 0
    call uninterpred(1, idum, lsmask_1d, lsmask_mdl, imdl, jmdl, ijmdl2, lonsperlat_mdl)
-   call uninterpred(1, idum, snow_cvr_mdl_1d, snow_cvr_mdl, imdl, jmdl, ijmdl2, lonsperlat_mdl)
+   call uninterpred(1, idum, snow_cvr_mdl_1d(:,km), snow_cvr_mdl, imdl, jmdl, ijmdl2, lonsperlat_mdl)
    deallocate(snow_cvr_mdl_1d)
-   call uninterpred(1, idum, snow_dep_mdl_tmp, snow_dep_mdl, imdl, jmdl, ijmdl2, lonsperlat_mdl)
+   call uninterpred(1, idum, snow_dep_mdl_tmp(:,km), snow_dep_mdl, imdl, jmdl, ijmdl2, lonsperlat_mdl)
    deallocate(snow_dep_mdl_tmp)
    deallocate(idum)
 
diff --git a/sorc/orog_mask_tools.fd/orog.fd/CMakeLists.txt b/sorc/orog_mask_tools.fd/orog.fd/CMakeLists.txt
index c292ec530..7cede544b 100644
--- a/sorc/orog_mask_tools.fd/orog.fd/CMakeLists.txt
+++ b/sorc/orog_mask_tools.fd/orog.fd/CMakeLists.txt
@@ -1,5 +1,5 @@
 set(lib_src netcdf_io.F90)
-set(exe_src mtnlm7_oclsm.f)
+set(exe_src mtnlm7_oclsm.F)
 
 if(CMAKE_Fortran_COMPILER_ID MATCHES "^(Intel)$")
   set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -r8 -convert big_endian -assume byterecl")
@@ -9,11 +9,14 @@ elseif(CMAKE_Fortran_COMPILER_ID MATCHES "^(GNU)$")
     set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -fallow-argument-mismatch -fallow-invalid-boz")
   endif()
 endif()
+if(ip_VERSION GREATER_EQUAL 4.0.0)
+  set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -DIP_V4")
+endif()
 
 set(exe_name orog)
 
 add_library(orog_lib STATIC ${lib_src})
-add_executable(${exe_name} mtnlm7_oclsm.f)
+add_executable(${exe_name} mtnlm7_oclsm.F)
 
 set(mod_dir "${CMAKE_CURRENT_BINARY_DIR}/mod")
 set_target_properties(orog_lib PROPERTIES Fortran_MODULE_DIRECTORY ${mod_dir})
diff --git a/sorc/orog_mask_tools.fd/orog.fd/mtnlm7_oclsm.f b/sorc/orog_mask_tools.fd/orog.fd/mtnlm7_oclsm.F
similarity index 99%
rename from sorc/orog_mask_tools.fd/orog.fd/mtnlm7_oclsm.f
rename to sorc/orog_mask_tools.fd/orog.fd/mtnlm7_oclsm.F
index 88d397303..438eefadf 100644
--- a/sorc/orog_mask_tools.fd/orog.fd/mtnlm7_oclsm.f
+++ b/sorc/orog_mask_tools.fd/orog.fd/mtnlm7_oclsm.F
@@ -3365,6 +3365,7 @@ SUBROUTINE MAKEOA2(ZAVG,zslm,VAR,GLAT,OA4,OL,IOA4,ELVMAX,
       WRITE(6,*) "! MAKEOA2 EXIT"
 C
       RETURN
+
       END SUBROUTINE MAKEOA2
  
 !> Compute a great circle distance between two points.
@@ -3553,6 +3554,12 @@ SUBROUTINE MAKEOA3(ZAVG,zslm,VAR,GLAT,OA4,OL,IOA4,ELVMAX,
      2           IM,JM,IMN,JMN,lon_c,lat_c,lon_t,lat_t,
      3           is_south_pole,is_north_pole,IMI,JMI,OA_IN,OL_IN,
      4           slm_in,lon_in,lat_in)
+
+! Required when using iplib v4.0 or higher.
+#ifdef IP_V4
+      use ipolates_mod
+#endif
+
       implicit none
       real, parameter :: MISSING_VALUE = -9999.
       real, parameter :: D2R = 3.14159265358979/180.
@@ -3592,13 +3599,13 @@ SUBROUTINE MAKEOA3(ZAVG,zslm,VAR,GLAT,OA4,OL,IOA4,ELVMAX,
       integer int_opt, ipopt(20), kgds_input(200), kgds_output(200)
       integer count_land_output
       integer ij, ijmdl_output, iret, num_mismatch_land, num
-      integer ibo(1)
+      integer ibo(1), ibi(1)
       logical*1, allocatable :: bitmap_input(:,:)
-      logical*1, allocatable :: bitmap_output(:)
+      logical*1, allocatable :: bitmap_output(:,:)
       integer, allocatable :: ijsav_land_output(:)
       real,    allocatable :: lats_land_output(:)
       real,    allocatable :: lons_land_output(:)
-      real,    allocatable :: output_data_land(:)
+      real,    allocatable :: output_data_land(:,:)
       real,    allocatable :: lons_mismatch_output(:)
       real,    allocatable :: lats_mismatch_output(:)
       real,    allocatable :: data_mismatch_output(:)
@@ -3735,8 +3742,8 @@ SUBROUTINE MAKEOA3(ZAVG,zslm,VAR,GLAT,OA4,OL,IOA4,ELVMAX,
       where(slm_in > 0.0) bitmap_input=.true.   
       print*, "count(bitmap_input)", count(bitmap_input)
       
-      allocate(bitmap_output(count_land_output))
-      allocate(output_data_land(count_land_output))
+      allocate(bitmap_output(count_land_output,1))
+      allocate(output_data_land(count_land_output,1))
       allocate(ijsav_land_output(count_land_output))
       allocate(lats_land_output(count_land_output))
       allocate(lons_land_output(count_land_output))
@@ -3757,13 +3764,14 @@ SUBROUTINE MAKEOA3(ZAVG,zslm,VAR,GLAT,OA4,OL,IOA4,ELVMAX,
 
       oa4 = 0.0
       ol = 0.0
+      ibi = 1
       
       do KWD=1,4
         bitmap_output = .false.
          output_data_land = 0.0
         call ipolates(int_opt, ipopt, kgds_input, kgds_output,   
      &         (IMI*JMI), count_land_output,               
-     &          1, 1, bitmap_input, oa_in(:,:,KWD),  
+     &          1, ibi, bitmap_input, oa_in(:,:,KWD),  
      &          count_land_output, lats_land_output,
      &          lons_land_output, ibo,  
      &          bitmap_output, output_data_land, iret)
@@ -3774,10 +3782,10 @@ SUBROUTINE MAKEOA3(ZAVG,zslm,VAR,GLAT,OA4,OL,IOA4,ELVMAX,
 
         num_mismatch_land = 0     
         do ij = 1, count_land_output
-          if (bitmap_output(ij)) then
+          if (bitmap_output(ij,1)) then
             j = (ijsav_land_output(ij)-1)/IM + 1
             i = mod(ijsav_land_output(ij)-1,IM) + 1
-            oa4(i,j,KWD)=output_data_land(ij)
+            oa4(i,j,KWD)=output_data_land(ij,1)
           else  ! default value
             num_mismatch_land =  num_mismatch_land + 1
           endif
@@ -3793,7 +3801,7 @@ SUBROUTINE MAKEOA3(ZAVG,zslm,VAR,GLAT,OA4,OL,IOA4,ELVMAX,
           
           num = 0     
           do ij = 1, count_land_output
-            if (.not. bitmap_output(ij)) then   
+            if (.not. bitmap_output(ij,1)) then   
               num = num+1
               lons_mismatch_output(num) = lons_land_output(ij)
               lats_mismatch_output(num) = lats_land_output(ij)
@@ -3815,7 +3823,7 @@ SUBROUTINE MAKEOA3(ZAVG,zslm,VAR,GLAT,OA4,OL,IOA4,ELVMAX,
        
         num = 0
         do ij = 1, count_land_output
-          if (.not. bitmap_output(ij)) then   
+          if (.not. bitmap_output(ij,1)) then   
             num = num+1
             j = (ijsav_land_output(ij)-1)/IM + 1
             i = mod(ijsav_land_output(ij)-1,IM) + 1
@@ -3835,7 +3843,7 @@ SUBROUTINE MAKEOA3(ZAVG,zslm,VAR,GLAT,OA4,OL,IOA4,ELVMAX,
         output_data_land = 0.0
         call ipolates(int_opt, ipopt, kgds_input, kgds_output,   
      &         (IMI*JMI), count_land_output,               
-     &          1, 1, bitmap_input, ol_in(:,:,KWD),  
+     &          1, ibi, bitmap_input, ol_in(:,:,KWD),  
      &          count_land_output, lats_land_output,
      &          lons_land_output, ibo,  
      &          bitmap_output, output_data_land, iret)
@@ -3846,10 +3854,10 @@ SUBROUTINE MAKEOA3(ZAVG,zslm,VAR,GLAT,OA4,OL,IOA4,ELVMAX,
 
         num_mismatch_land = 0     
         do ij = 1, count_land_output
-          if (bitmap_output(ij)) then
+          if (bitmap_output(ij,1)) then
             j = (ijsav_land_output(ij)-1)/IM + 1
             i = mod(ijsav_land_output(ij)-1,IM) + 1
-            ol(i,j,KWD)=output_data_land(ij)
+            ol(i,j,KWD)=output_data_land(ij,1)
           else  ! default value
             num_mismatch_land =  num_mismatch_land + 1
           endif
@@ -3862,7 +3870,7 @@ SUBROUTINE MAKEOA3(ZAVG,zslm,VAR,GLAT,OA4,OL,IOA4,ELVMAX,
        
         num = 0
         do ij = 1, count_land_output
-          if (.not. bitmap_output(ij)) then   
+          if (.not. bitmap_output(ij,1)) then   
             num = num+1
             j = (ijsav_land_output(ij)-1)/IM + 1
             i = mod(ijsav_land_output(ij)-1,IM) + 1