Implement MPAS init-atmosphere case 11 + 12 for model-level IFS and ERA5 with lapse-rate vinterp 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. For model-level ECMWF data e.g. ERA5 and IFS, also implemented init-atmosphere case 11 (ECMWF initial condition) and case 12 (ECMWF LBC).

Files Modified:

src/core_init_atmosphere/mpas_init_atm_cases.F (+392 lines, -36 lines)
src/core_init_atmosphere/mpas_init_atm_core_interface.F (+38 lines, -2 lines)
src/core_init_atmosphere/mpas_init_atm_vinterp.F (+200 lines, -5 lines)
src/core_init_atmosphere/Registry.xml (+4 lines, -2 lines)

Testing:

Successfully built MPAS model (both init-atmosphere and atmosphere model) with modified code. Test simulations were performed using the built init_atmosphere_model and atmosphere_model executables at 240km uniform mesh. 48h forecasts were successfully generated using era5 model level initial conditions. ERA5 intermediate files were generated with ERA5toInt tool.

Init-atmosphere namelist specifications for test u240k experiment:

config_init_case = 11
config_extrap_airtemp = 'lapse-rate'
config_nvertlevels = 66
config_ztop = 40000.0

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.

Copilot

Pull request overview

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

Comments suppressed due to low confidence (1)

src/core_init_atmosphere/mpas_init_atm_core_interface.F:206

The PR description/metadata focuses on implementing extrap_type==2 for vertical interpolation, but this PR also introduces new init cases (11/12), adds a new vertical_interp_ecmwf API, changes package activation logic, and renames a Registry field (gfs_z → ght_fg). Please update the PR description (and any user-facing documentation) to reflect these additional behavioral/API changes so reviewers and users can assess impact and backward compatibility.

      if (config_init_case == 9 .or. config_init_case == 12) then
         lbcs = .true.
         mp_thompson_aers_in = .false.
         inquire(file="QNWFA_QNIFA_SIGMA_MONTHLY.dat",exist=lexist)
         if(lexist) mp_thompson_aers_in = .true.
      else
         lbcs = .false.
         mp_thompson_aers_in = .false.
      end if

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

+            if (use_ecmwf) then
+               u(k,iEdge) = vertical_interp_ecmwf(target_z, nfglevels_actual-1, sorted_arr(:,1:nfglevels_actual-1))
+            else
+               u(k,iEdge) = vertical_interp(target_z, nfglevels_actual-1, sorted_arr(:,1:nfglevels_actual-1), order=1, extrap=2)


+      if (.not. found) then
+         ! Degenerate sorted data (duplicate adjacent heights); return nearest level.
+         if (present(ierr)) ierr = 1
+         lm = 1
+         do k = 2, nz
+            if (abs(zf(1,k) - target_z) < abs(zf(1,lm) - target_z)) lm = k
+         end do
+         vertical_interp_ecmwf = zf(2,lm)
+         return


+         ! Relative humidity: if first-guess values are negative, set those values to zero before
+         ! vertical interpolation.


+                                          sorted_arr(:,1:nfglevels_actual-1), order=1, extrap=2)
+            end if
+            ! RH is physically bounded to [0, 100]; clamp after LS extrapolation


+                                                        sorted_arr(:,1:nfglevels_actual-1))
+            else
+               ght_fg(k,iCell) = vertical_interp(target_z, nfglevels_actual-1, &
+                                          sorted_arr(:,1:nfglevels_actual-1), order=1, extrap=2)


+            if (use_ecmwf) then
+               psfc(iCell) = exp(vertical_interp_ecmwf(target_z, nfglevels_actual, sorted_arr))
+            else
+               psfc(iCell) = exp(vertical_interp(target_z, nfglevels_actual, sorted_arr, order=1, extrap=2))


+               relhum(k,iCell) = vertical_interp(target_z, nfglevels_actual-1, &
+                                          sorted_arr(:,1:nfglevels_actual-1), order=1, extrap=2)
+            end if
+            ! RH is physically bounded to [0, 100]; clamp after LS extrapolation
+            relhum(k,iCell) = max(0.0_RKIND, min(100.0_RKIND, relhum(k,iCell)))
         end do


+                                        sorted_arr(:,1:nfglevels_actual-1)))
+            else
+               pressure(k,iCell) = exp(vertical_interp(target_z, nfglevels_actual-1, &
+                                       sorted_arr(:,1:nfglevels_actual-1), order=1, extrap=2))


                                        8 = surface field (SST, sea-ice) update file for use with real-data simulations \newline
                                        9 = lateral boundary conditions update file for use with real-data simulations \newline
                                        10 = idealized LES case \newline
+                                        11 = real-data initial conditions from ECMWF model levels (ERA5/IFS) \newline
+                                        12 = lateral boundary conditions from ECMWF model levels (ERA5/IFS) \newline
                                        13 = CAM-MPAS 3-d grid with specified topography and zeta levels"
-                     possible_values="1 -- 10, or 13"/>
+                     possible_values="1 -- 13"/>


                     packages="met_stage_out"/>

-                <var name="gfs_z" type="real" dimensions="nVertLevels nCells Time" units="m"
+                <var name="ght_fg" type="real" dimensions="nVertLevels nCells Time" units="m"


mrixlam · 2026-05-04T03:50:51Z

Hi @mgduda , I think code changes are in a stage where we can request your review. MPAS model was successfully built with modified code and 48h forecasts were generated with the built model executables. As for the vinterp, I am attaching three panels with temperature, wind, humidity and pressure interpolated with the least square based method. If you would like to check the MPAS build and run directory, it can be found here: /glade/derecho/scratch/mrislam/work/pandac/mpas/gifsim/gifsml

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

Add ECMWF model level initialization

34306d3

mrixlam requested a review from Copilot May 3, 2026 07:15

Copilot started reviewing on behalf of mrixlam May 3, 2026 07:15 View session

Copilot AI reviewed May 3, 2026

View reviewed changes

mrixlam changed the title ~~Implement lapse-rate extrapolation for vertical interpolation (extrap_type == 2)~~ Implement Init-atmosphere case 11 + 12 for model-level IFS and ERA5 with lapse-rate vinterp May 4, 2026

mrixlam changed the title ~~Implement Init-atmosphere case 11 + 12 for model-level IFS and ERA5 with lapse-rate vinterp~~ Implement MPAS init-atmosphere case 11 + 12 for model-level IFS and ERA5 with lapse-rate vinterp May 4, 2026

		! Relative humidity: if first-guess values are negative, set those values to zero before
		! vertical interpolation.

Conversation

mrixlam commented Feb 12, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Summary:

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!

Copilot AI left a comment

Choose a reason for hiding this comment

Pull request overview

Uh oh!

mrixlam commented May 4, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

mrixlam commented Feb 12, 2026 •

edited

Loading