Implement lapse-rate extrapolation for vertical interpolation (extrap_type == 2) by mrixlam · Pull Request #1413 · MPAS-Dev/MPAS-Model

mrixlam · 2026-02-12T07:02:27Z

Summary:

This PR implements the missing lapse-rate extrapolation functionality (extrap_type == 2) for vertical interpolation in the init_atmosphere module. Previously, this extrapolation type was documented but not implemented for the upper boundary condition, causing fatal errors when model top exceeded input data vertical range.

Problem Statement:

When initializing MPAS-Atmosphere with input data (e.g., IFS/ECMWF) that has a top level below or near the model's configured vertical extent (config_ztop), vertical interpolation may require extrapolation beyond available data levels. While the code supported three extrapolation types:

Type 0: Constant (use nearest value)
Type 1: Linear (extend using slope between last two points)
Type 2: Lapse-rate (use atmospheric temperature gradient)

Type 2 (lapse-rate) was not implemented for the upper boundary (target_z >= zf(1,nz)), resulting in this fatal error:

ERROR: extrap_type == 2 not implemented for target_z >= zf(1,nz)
CRITICAL ERROR: Error in interpolation of t(k,iCell) for k=1, iCell=1

This prevented users from using the physically more appropriate lapse-rate extrapolation method via config_extrap_airtemp = 'lapse-rate' in namelist.init_atmosphere.

Files Modified:

src/core_init_atmosphere/mpas_init_atm_vinterp.F (+5 lines, -0 lines)
src/core_init_atmosphere/mpas_init_atm_cases.F (+4 lines, -3 lines)

Testing:

Successfully tested MPAS model initialization with quasi-uniform global 15km quasi-uniform mesh with following specifications:

IFS 137-level input data (model top ~80 km)
Model configuration: 66 vertical levels, 40 km top
config_extrap_airtemp = 'lapse-rate' in namelist

Copilot

Pull request overview

Implements the previously missing extrap_type == 2 (“lapse-rate”) behavior for vertical extrapolation at the upper boundary during init-atmosphere vertical interpolation, preventing fatal errors when the model top exceeds the input data’s vertical range.

Changes:

Extend vertical_interp to handle extrap_type == 2 for target_z >= zf(1,nz) (upper-boundary extrapolation).
Add/standardize an optional ierr output on the init_atm_vinterp implementation to signal invalid extrapolation types.
Remove the prior “not implemented” fatal-path for extrap_type == 2 in mpas_init_atm_cases.F.

Reviewed changes

Copilot reviewed 2 out of 2 changed files in this pull request and generated 3 comments.

File	Description
`src/core_init_atmosphere/mpas_init_atm_vinterp.F`	Adds optional `ierr` and implements `extrap_type==2` for upper/lower boundary extrapolation.
`src/core_init_atmosphere/mpas_init_atm_cases.F`	Replaces the previous fatal error for upper-boundary `extrap_type==2` with an extrapolation calculation.

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

Copilot · 2026-02-12T07:07:00Z

+         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))


extrap_type == 2 (lapse-rate) is implemented here using the same formula as extrap_type == 1 (linear), i.e., a slope computed from the bottom two levels. That makes the new lapse-rate option behaviorally identical to linear extrapolation. If lapse-rate is intended to be a fixed temperature lapse rate (as in the other vertical_interp implementation in mpas_init_atm_cases.F), please apply that constant lapse rate here (or otherwise document/rename the option so it’s not indistinguishable from linear).

Copilot · 2026-02-12T07:07:01Z

+         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))


Upper-boundary extrap_type == 2 uses a lapse rate computed from the top two levels, which is mathematically equivalent to the existing extrap_type == 1 linear extrapolation. If lapse-rate is meant to behave differently from linear (e.g., fixed -6.5 K/km), this should be adjusted so the two extrapolation modes are not duplicates.

Copilot · 2026-02-12T07:07:01Z

+            ! 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))


For extrap_type == 2 (lapse-rate) at the upper boundary, this change computes lapse_rate from the top two levels, which makes it effectively the same as extrap_type == 1 linear extrapolation. In this same function, the lower-boundary lapse-rate branch uses a fixed 0.0065 K/m (line 6926), so the upper and lower implementations are now inconsistent. Please align the two (either both fixed-lapse-rate or both gradient-derived) and ensure linear vs lapse-rate remain meaningfully different options.

Suggested change

! 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))

! Lapse-rate extrapolation: use fixed lapse rate (0.0065 K/m)

vertical_interp = zf(2,nz) - (target_z - zf(1,nz))*0.0065

mgduda · 2026-03-24T16:04:33Z

@mrixlam While looking through the changes in this PR, I was reminded that we in fact have two vertical_interp functions in the init_atmosphere core: one in mpas_init_atm_vinterp.F and the other in mpas_init_atm_cases.F. Before considering the changes in this PR, we may try to refactor the init_atmosphere code so that there is only a single vertical_interp routine. After we've done this, I'll follow up here so that we can figure out next steps for this PR.

mgduda · 2026-03-25T22:24:42Z

@mrixlam We've just merged PR #1429, which consolidates the two different definitions of the vertical_interp function. Would you be able to rework your changes based on merge commit 88ae401, modifying the vertical_interp function in src/core_init_atmosphere/mpas_init_atm_vinterp.F?

Copilot

Pull request overview

Copilot reviewed 2 out of 2 changed files in this pull request and generated 2 comments.

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

Copilot · 2026-04-29T08:02:26Z

         else if (extrap_type == 2) then
-            vertical_interp = zf(2,1) - (target_z - zf(1,1))*0.0065
+            n_lapse = min(nz, 3)
+            sum_x = 0.0_RKIND; sum_y = 0.0_RKIND
+            sum_x2 = 0.0_RKIND; sum_xy = 0.0_RKIND
+            do k = 1, n_lapse
+               sum_x  = sum_x  + zf(1,k)


The lapse-rate extrapolation for the lower boundary (target_z < zf(1,1)) used to apply a fixed standard lapse rate (0.0065 K/m). This change switches extrap_type==2 to a least-squares slope fit from the lowest 1–3 input levels, which is a behavioral change beyond the PR description (which focuses on missing upper boundary support). Please confirm this is intended; if not, keep the fixed lapse-rate behavior for the lower boundary (and implement the same for the upper boundary) or update the PR description/namelist documentation accordingly.

The hard-code 0.0065 K/m threshold is a field-specific constant that only produce a reasonable extrapolation when it is done for temperature in K. However, this is a field-agnostic function that is also used to extrapolate humidity, for example, which might produce unphysical extrapolated values. We will have a change in PR description before merging it.

Copilot · 2026-04-29T08:02:26Z

+            n_lapse = min(nz, 3)
+            sum_x = 0.0_RKIND; sum_y = 0.0_RKIND
+            sum_x2 = 0.0_RKIND; sum_xy = 0.0_RKIND
+            do k = nz - n_lapse + 1, nz
+               sum_x  = sum_x  + zf(1,k)


The least-squares slope calculation for extrap_type==2 is duplicated in both the lower- and upper-boundary extrapolation branches. Consider extracting the slope computation into a small helper routine/function (or a shared local block) so future adjustments (e.g., number of points, guarding against degenerate z spacing) don’t need to be made in two places.

…m/MPAS-Model into bugfix-vinterp-lapse-rate

Copilot

Pull request overview

Copilot reviewed 2 out of 2 changed files in this pull request and generated 4 comments.

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

Copilot · 2026-04-29T09:26:43Z

+      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


ls_lapse_slope divides by (n_real*sum_x2 - sum_x**2) without guarding against a zero/near-zero denominator. If the selected z-levels contain duplicate (or nearly duplicate) zf(1,k) values, this will trigger a divide-by-zero (or overflow to Inf/NaN) and silently propagate into the extrapolated field. Please add a tolerance check on the denominator (and return a safe fallback slope, e.g., 0 or a 2-point difference), and consider computing the variance term using centered values to avoid catastrophic cancellation when RKIND is single precision.

Suggested change

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

real (kind=RKIND) :: sum_x, sum_y, mean_x, mean_y

real (kind=RKIND) :: dx, dy, var_x, cov_xy, tol, scale_x, denom_2pt

real (kind=RKIND) :: n_real

n_pts = k_end - k_start + 1

if (n_pts <= 1) then

ls_lapse_slope = 0.0_RKIND

return

end if

sum_x = 0.0_RKIND

sum_y = 0.0_RKIND

do k = k_start, k_end

sum_x = sum_x + zf(1,k)

sum_y = sum_y + zf(2,k)

end do

n_real = real(n_pts, RKIND)

mean_x = sum_x / n_real

mean_y = sum_y / n_real

var_x = 0.0_RKIND

cov_xy = 0.0_RKIND

scale_x = 0.0_RKIND

do k = k_start, k_end

dx = zf(1,k) - mean_x

dy = zf(2,k) - mean_y

var_x = var_x + dx*dx

cov_xy = cov_xy + dx*dy

scale_x = max(scale_x, abs(zf(1,k)))

end do

tol = max(1.0_RKIND, scale_x)**2 * epsilon(1.0_RKIND) * n_real

if (var_x > tol) then

ls_lapse_slope = cov_xy / var_x

else

denom_2pt = zf(1,k_end) - zf(1,k_start)

if (abs(denom_2pt) > tol) then

ls_lapse_slope = (zf(2,k_end) - zf(2,k_start)) / denom_2pt

else

ls_lapse_slope = 0.0_RKIND

end if

Copilot · 2026-04-29T09:26:43Z

      if (target_z < zf(1,1)) then
         if (extrap_type == 0) then
            vertical_interp = zf(2,1)
         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
-            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
      if (target_z >= zf(1,nz)) then
         if (extrap_type == 0) then
            vertical_interp = zf(2,nz)
         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))


In the extrapolation branches, the linear extrapolation path assumes nz >= 2 (accesses zf(:,2) / zf(:,nz-1)). If the caller provides only a single vertical level (e.g., nz==1), this becomes an out-of-bounds access. Please add an explicit nz < 2 fallback (e.g., treat as constant extrapolation and/or set ierr + log an error) before computing the slope for extrap_type == 1.

Copilot · 2026-04-29T09:26:44Z

      if (target_z < zf(1,1)) then
         if (extrap_type == 0) then
            vertical_interp = zf(2,1)
         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
-            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
      if (target_z >= zf(1,nz)) then
         if (extrap_type == 0) then
            vertical_interp = zf(2,nz)
         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
-             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


The extrapolation if/else if chain has no final else to handle unexpected extrap_type values. If a caller passes something other than {0,1,2}, vertical_interp can be returned uninitialized while ierr remains 0. Please add a default branch that logs an error and sets ierr (and/or falls back to constant extrapolation).

Copilot · 2026-04-29T09:26:44Z

         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


The PR description focuses on adding missing upper-boundary support for extrap_type==2, but this diff also changes the lower-boundary extrap_type==2 behavior from a fixed 0.0065 K/m to a least-squares fitted slope from the lowest levels. Please update the PR description (and any user-facing docs for config_extrap_airtemp='lapse-rate') to reflect this broader behavior change.

Fix lapse-rate vinterp for EC

fa48c91

Copilot AI review requested due to automatic review settings February 12, 2026 07:02

Copilot started reviewing on behalf of mrixlam February 12, 2026 07:02 View session

Copilot AI reviewed Feb 12, 2026

View reviewed changes

abishekg7 mentioned this pull request Mar 25, 2026

Consolidating multiple definitions of vertical_interp into init_atm_vinterp #1429

Merged

mgduda changed the base branch from master to develop March 25, 2026 22:22

mrixlam and others added 2 commits April 29, 2026 01:30

bugfix vinterp

a2fd267

Merge branch 'develop' into bugfix-vinterp-lapse-rate

86845cb

mrixlam requested a review from Copilot April 29, 2026 07:57

Copilot started reviewing on behalf of mrixlam April 29, 2026 07:58 View session

Copilot AI reviewed Apr 29, 2026

View reviewed changes

mrixlam added 2 commits April 29, 2026 02:36

Add function for least square slope estimation

924b456

Merge branch 'bugfix-vinterp-lapse-rate' of https://github.com/mrixla…

ef43c50

…m/MPAS-Model into bugfix-vinterp-lapse-rate

mrixlam closed this Apr 29, 2026

mrixlam deleted the bugfix-vinterp-lapse-rate branch April 29, 2026 08:42

mrixlam restored the bugfix-vinterp-lapse-rate branch April 29, 2026 09:02

mrixlam reopened this Apr 29, 2026

mrixlam requested a review from Copilot April 29, 2026 09:22

Copilot started reviewing on behalf of mrixlam April 29, 2026 09:23 View session

Copilot AI reviewed Apr 29, 2026

View reviewed changes

Conversation

mrixlam commented Feb 12, 2026

Summary:

Problem Statement:

Files Modified:

Testing:

Uh oh!

Copilot AI left a comment

Choose a reason for hiding this comment

Pull request overview

Reviewed changes

Uh oh!

Copilot AI Feb 12, 2026

Choose a reason for hiding this comment

Uh oh!

Copilot AI Feb 12, 2026

Choose a reason for hiding this comment

Uh oh!

Copilot AI Feb 12, 2026

Choose a reason for hiding this comment

Uh oh!

mgduda commented Mar 24, 2026

Uh oh!

mgduda commented Mar 25, 2026

Uh oh!

Copilot AI left a comment

Choose a reason for hiding this comment

Pull request overview

Uh oh!

Copilot AI Apr 29, 2026

Choose a reason for hiding this comment

Uh oh!

mrixlam Apr 29, 2026

Choose a reason for hiding this comment

Uh oh!

Copilot AI Apr 29, 2026

Choose a reason for hiding this comment

Uh oh!

Copilot AI left a comment

Choose a reason for hiding this comment

Pull request overview

Uh oh!

Copilot AI Apr 29, 2026

Choose a reason for hiding this comment

Uh oh!

Copilot AI Apr 29, 2026

Choose a reason for hiding this comment

Uh oh!

Copilot AI Apr 29, 2026

Choose a reason for hiding this comment

Uh oh!

Copilot AI Apr 29, 2026

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants