From fa48c91a4d832fa4cca3ff318bc4a06cb319d03c Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?=E2=80=9Cmrixlam=E2=80=9D?= <“mrislam386@gmail.com”>
Date: Wed, 11 Feb 2026 23:33:07 -0700
Subject: [PATCH] Fix lapse-rate vinterp for EC

---
 .../mpas_init_atm_cases.F                     |  8 ++++---
 .../mpas_init_atm_vinterp.F                   | 24 ++++++++++++++++++-
 2 files changed, 28 insertions(+), 4 deletions(-)

diff --git a/src/core_init_atmosphere/mpas_init_atm_cases.F b/src/core_init_atmosphere/mpas_init_atm_cases.F
index 673ebfc525..36dccd8adb 100644
--- a/src/core_init_atmosphere/mpas_init_atm_cases.F
+++ b/src/core_init_atmosphere/mpas_init_atm_cases.F
@@ -6882,6 +6882,7 @@ real (kind=RKIND) function vertical_interp(target_z, nz, zf, order, extrap, surf
       integer :: k, lm, lp
       real (kind=RKIND) :: wm, wp
       real (kind=RKIND) :: slope
+      real (kind=RKIND) :: lapse_rate
 
       integer :: interp_order, extrap_type
       real (kind=RKIND) :: surface, sealevel
@@ -6933,9 +6934,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
+            ! Lapse-rate extrapolation: calculate lapse rate from top two levels
+            lapse_rate = (zf(2,nz) - zf(2,nz-1)) / (zf(1,nz) - zf(1,nz-1))
+            vertical_interp = zf(2,nz) + lapse_rate * (target_z - zf(1,nz))
          end if
          return
       end if
@@ -7522,3 +7523,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 c5163d2708..e39a496066 100644
--- a/src/core_init_atmosphere/mpas_init_atm_vinterp.F
+++ b/src/core_init_atmosphere/mpas_init_atm_vinterp.F
@@ -22,7 +22,7 @@ module init_atm_vinterp
    !
    ! Purpose:
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-   real (kind=RKIND) function vertical_interp(target_z, nz, zf, order, extrap, surface_val, sealev_val)
+   real (kind=RKIND) function vertical_interp(target_z, nz, zf, order, extrap, surface_val, sealev_val, ierr)
 
       implicit none
 
@@ -33,10 +33,12 @@ real (kind=RKIND) function vertical_interp(target_z, nz, zf, order, extrap, surf
       integer, intent(in), optional :: extrap
       real (kind=RKIND), intent(in), optional :: surface_val
       real (kind=RKIND), intent(in), optional :: sealev_val
+      integer, intent(out), optional :: ierr                    ! error code: 0 = success, 1 = invalid extrapolation type
 
       integer :: k, lm, lp
       real (kind=RKIND) :: wm, wp
       real (kind=RKIND) :: slope
+      real (kind=RKIND) :: lapse_rate
 
       integer :: interp_order, extrap_type
       real (kind=RKIND) :: surface, sealevel
@@ -66,6 +68,11 @@ real (kind=RKIND) function vertical_interp(target_z, nz, zf, order, extrap, surf
          sealevel = 201300.0
       end if
 
+      ! Initialize ierr to 0 (success)
+      if (present(ierr)) then
+         ierr = 0
+      end if
+
       !
       ! Extrapolation required
       !
@@ -75,6 +82,13 @@ real (kind=RKIND) function vertical_interp(target_z, nz, zf, order, extrap, surf
          else if (extrap_type == 1) then
             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
+            ! Lapse-rate extrapolation: calculate lapse rate from bottom two levels
+            lapse_rate = (zf(2,2) - zf(2,1)) / (zf(1,2) - zf(1,1))
+            vertical_interp = zf(2,1) + lapse_rate * (target_z - zf(1,1))
+         else
+            vertical_interp = zf(2,1)
+            if (present(ierr)) ierr = 1
          end if
          return
       end if
@@ -84,6 +98,13 @@ real (kind=RKIND) function vertical_interp(target_z, nz, zf, order, extrap, surf
          else if (extrap_type == 1) then
             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
+            ! Lapse-rate extrapolation: calculate lapse rate from top two levels
+            lapse_rate = (zf(2,nz) - zf(2,nz-1)) / (zf(1,nz) - zf(1,nz-1))
+            vertical_interp = zf(2,nz) + lapse_rate * (target_z - zf(1,nz))
+         else
+            vertical_interp = zf(2,nz)
+            if (present(ierr)) ierr = 1
          end if
          return
       end if
@@ -109,3 +130,4 @@ real (kind=RKIND) function vertical_interp(target_z, nz, zf, order, extrap, surf
    end function vertical_interp
 
 end module init_atm_vinterp
+