diff --git a/src/core_init_atmosphere/mpas_init_atm_cases.F b/src/core_init_atmosphere/mpas_init_atm_cases.F
index fb67a5164b..9f485804f2 100644
--- a/src/core_init_atmosphere/mpas_init_atm_cases.F
+++ b/src/core_init_atmosphere/mpas_init_atm_cases.F
@@ -8154,3 +8154,4 @@ end function read_text_array
 
 
 end module init_atm_cases
+
diff --git a/src/core_init_atmosphere/mpas_init_atm_vinterp.F b/src/core_init_atmosphere/mpas_init_atm_vinterp.F
index 164cfac546..a88505cb4a 100644
--- a/src/core_init_atmosphere/mpas_init_atm_vinterp.F
+++ b/src/core_init_atmosphere/mpas_init_atm_vinterp.F
@@ -38,7 +38,7 @@ real (kind=RKIND) function vertical_interp(target_z, nz, zf, order, extrap, surf
       real (kind=RKIND), intent(in), optional :: sealev_val
       integer, intent(out), optional :: ierr
       
-      integer :: k, lm, lp
+      integer :: k, lm, lp, n_lapse
       real (kind=RKIND) :: wm, wp
       real (kind=RKIND) :: slope
 
@@ -81,7 +81,9 @@ real (kind=RKIND) function vertical_interp(target_z, nz, zf, order, extrap, surf
             slope = (zf(2,2) - zf(2,1)) / (zf(1,2) - zf(1,1))
             vertical_interp = zf(2,1) + slope * (target_z - zf(1,1))
          else if (extrap_type == 2) then
-            vertical_interp = zf(2,1) - (target_z - zf(1,1))*0.0065
+            n_lapse = min(nz, 3)
+            slope = ls_lapse_slope(zf, nz, 1, n_lapse)
+            vertical_interp = zf(2,1) + slope * (target_z - zf(1,1))
          end if
          return
       end if
@@ -92,9 +94,9 @@ real (kind=RKIND) function vertical_interp(target_z, nz, zf, order, extrap, surf
             slope = (zf(2,nz) - zf(2,nz-1)) / (zf(1,nz) - zf(1,nz-1))
             vertical_interp = zf(2,nz) + slope * (target_z - zf(1,nz))
          else if (extrap_type == 2) then
-             call mpas_log_write('extrap_type == 2 not implemented for target_z >= zf(1,nz)', messageType=MPAS_LOG_ERR)
-             if (present(ierr)) ierr = 1
-             return
+            n_lapse = min(nz, 3)
+            slope = ls_lapse_slope(zf, nz, nz - n_lapse + 1, nz)
+            vertical_interp = zf(2,nz) + slope * (target_z - zf(1,nz))
          end if
          return
       end if
@@ -119,4 +121,43 @@ real (kind=RKIND) function vertical_interp(target_z, nz, zf, order, extrap, surf
 
    end function vertical_interp
 
+   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+   ! Name: ls_lapse_slope
+   !
+   ! Purpose: Compute an ordinary least-squares slope over levels k_start..k_end of
+   !          zf.  Returns 0 when only one level is supplied (constant fallback).
+   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+   real (kind=RKIND) function ls_lapse_slope(zf, nz, k_start, k_end)
+
+      implicit none
+
+      integer, intent(in) :: nz, k_start, k_end
+      real (kind=RKIND), dimension(2,nz), intent(in) :: zf
+
+      integer :: k, n_pts
+      real (kind=RKIND) :: sum_x, sum_y, sum_x2, sum_xy, n_real
+
+      n_pts  = k_end - k_start + 1
+
+      sum_x  = 0.0_RKIND;  sum_y  = 0.0_RKIND
+      sum_x2 = 0.0_RKIND;  sum_xy = 0.0_RKIND
+
+      do k = k_start, k_end
+         sum_x  = sum_x  + zf(1,k)
+         sum_y  = sum_y  + zf(2,k)
+         sum_x2 = sum_x2 + zf(1,k)**2
+         sum_xy = sum_xy + zf(1,k)*zf(2,k)
+      end do
+
+      n_real = real(n_pts, RKIND)
+      
+      if (n_pts > 1) then
+         ls_lapse_slope = (n_real*sum_xy - sum_x*sum_y) / (n_real*sum_x2 - sum_x**2)
+      else
+         ls_lapse_slope = 0.0_RKIND
+      end if
+
+   end function ls_lapse_slope
+
 end module init_atm_vinterp
+