diff --git a/src/core_init_atmosphere/Registry.xml b/src/core_init_atmosphere/Registry.xml
index d45a115e33..a06165ab64 100644
--- a/src/core_init_atmosphere/Registry.xml
+++ b/src/core_init_atmosphere/Registry.xml
@@ -66,8 +66,10 @@
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"/>
-
diff --git a/src/core_init_atmosphere/mpas_init_atm_cases.F b/src/core_init_atmosphere/mpas_init_atm_cases.F
index fb67a5164b..6a224a68c9 100644
--- a/src/core_init_atmosphere/mpas_init_atm_cases.F
+++ b/src/core_init_atmosphere/mpas_init_atm_cases.F
@@ -16,7 +16,7 @@ module init_atm_cases
use atm_advection
use mpas_RBF_interpolation
use mpas_vector_reconstruction
- use init_atm_vinterp, only: vertical_interp
+ use init_atm_vinterp, only: vertical_interp, vertical_interp_ecmwf
use mpas_timer
use mpas_init_atm_static
use mpas_init_atm_surface
@@ -389,6 +389,155 @@ subroutine init_atm_setup_case(domain, stream_manager)
block_ptr => block_ptr % next
end do
+ else if (config_init_case == 11) then
+
+ call mpas_log_write('Real-data ECMWF model-level test case (ERA5/IFS)')
+ block_ptr => domain % blocklist
+
+ do while (associated(block_ptr))
+
+ call mpas_pool_get_config(block_ptr % configs, 'config_static_interp', config_static_interp)
+ call mpas_pool_get_config(block_ptr % configs, 'config_native_gwd_static', config_native_gwd_static)
+ call mpas_pool_get_config(block_ptr % configs, 'config_native_gwd_gsl_static', config_native_gwd_gsl_static)
+ call mpas_pool_get_config(block_ptr % configs, 'config_met_interp', config_met_interp)
+ call mpas_pool_get_config(block_ptr % configs, 'config_blend_bdy_terrain', config_blend_bdy_terrain)
+
+ call mpas_pool_get_subpool(block_ptr % structs, 'mesh', mesh)
+ call mpas_pool_get_subpool(block_ptr % structs, 'fg', fg)
+ call mpas_pool_get_subpool(block_ptr % structs, 'state', state)
+ call mpas_pool_get_subpool(block_ptr % structs, 'diag', diag)
+ call mpas_pool_get_subpool(block_ptr % structs, 'diag_physics', diag_physics)
+
+ call mpas_pool_get_dimension(block_ptr % dimensions, 'nCells', nCells)
+ call mpas_pool_get_dimension(block_ptr % dimensions, 'nEdges', nEdges)
+ call mpas_pool_get_dimension(block_ptr % dimensions, 'nVertLevels', nVertLevels)
+
+ if (config_blend_bdy_terrain) then
+ call mpas_pool_get_config(block_ptr % configs, 'config_start_time', config_start_time)
+ call mpas_pool_get_config(block_ptr % configs, 'config_met_prefix', config_met_prefix)
+
+ call blend_bdy_terrain(config_met_prefix, config_start_time, &
+ nCells, latCell, lonCell, bdyMaskCell, ter, .true., ierr)
+ if (ierr /= 0) then
+ call mpas_log_write('*************************************************************', messageType=MPAS_LOG_ERR)
+ call mpas_log_write('Blending of terrain along domain boundaries would fail, and', messageType=MPAS_LOG_ERR)
+ call mpas_log_write('config_blend_bdy_terrain = true in the namelist.init_atmosphere file.', &
+ messageType=MPAS_LOG_ERR)
+ call mpas_log_write('*************************************************************', messageType=MPAS_LOG_CRIT)
+ end if
+ end if
+
+ if (config_static_interp) then
+ call init_atm_static(mesh, block_ptr % dimensions, block_ptr % configs)
+ end if
+
+ if (config_native_gwd_static) then
+ call mpas_log_write(' ')
+ call mpas_log_write('Computing GWDO static fields on the native MPAS mesh')
+ call mpas_log_write(' ')
+ ierr = compute_gwd_fields(domain)
+ if (ierr /= 0) then
+ call mpas_log_write('****************************************************************', messageType=MPAS_LOG_ERR)
+ call mpas_log_write('Error while trying to compute sub-grid-scale orography', messageType=MPAS_LOG_ERR)
+ call mpas_log_write('statistics for use with the GWDO scheme.', messageType=MPAS_LOG_ERR)
+ call mpas_log_write('****************************************************************', messageType=MPAS_LOG_CRIT)
+ end if
+ end if
+
+ if (config_native_gwd_gsl_static) then
+ call mpas_log_write(' ')
+ call mpas_log_write('Computing GWDO static fields for UGWP orog drag on the native MPAS mesh')
+ call mpas_log_write(' ')
+ call calc_gsl_oro_data(domain,ierr)
+ if (ierr /= 0) then
+ call mpas_log_write('****************************************************************', messageType=MPAS_LOG_ERR)
+ call mpas_log_write('Error while trying to compute sub-grid-scale GSL orography', messageType=MPAS_LOG_ERR)
+ call mpas_log_write('statistics for use with the GWDO scheme.', messageType=MPAS_LOG_ERR)
+ call mpas_log_write('****************************************************************', messageType=MPAS_LOG_CRIT)
+ end if
+ end if
+
+ call mpas_pool_get_array(mesh, 'mminlu', mminlu)
+ if (len_trim(mminlu) == 0) then
+ call mpas_log_write('****************************************************************')
+ call mpas_log_write('No information on land use dataset is available.')
+ call mpas_log_write('Assume that we are using ''USGS''.')
+ call mpas_log_write('****************************************************************')
+ write(mminlu,'(a)') 'USGS'
+ end if
+
+ call init_atm_case_ecmwf(block_ptr, mesh, nCells, nEdges, nVertLevels, fg, state, &
+ diag, diag_physics, block_ptr % dimensions, block_ptr % configs)
+
+ if (config_met_interp) then
+ call init_atm_thompson_aerosols(block_ptr, mesh, block_ptr % configs, diag, state)
+ call physics_initialize_real(mesh, fg, domain % dminfo, block_ptr % dimensions, block_ptr % configs)
+ end if
+
+ block_ptr => block_ptr % next
+ end do
+
+ else if (config_init_case == 12) then
+
+ call mpas_log_write('ECMWF model-level LBC case (ERA5/IFS)')
+
+ clock_interval = mpas_get_clock_timestep(domain % clock, ierr=ierr)
+ lbc_stream_interval = MPAS_stream_mgr_get_stream_interval(stream_manager, 'lbc', MPAS_STREAM_OUTPUT, ierr)
+ if (clock_interval /= lbc_stream_interval) then
+ call mpas_log_write('****************************************************************', messageType=MPAS_LOG_ERR)
+ call mpas_log_write('The intermediate file interval specified by ''config_fg_interval''', messageType=MPAS_LOG_ERR)
+ call mpas_log_write('does not match the output_interval for the ''lbc'' stream.', messageType=MPAS_LOG_ERR)
+ call mpas_log_write('Please correct the namelist.init_atmosphere and/or', messageType=MPAS_LOG_ERR)
+ call mpas_log_write('streams.init_atmosphere files.', messageType=MPAS_LOG_ERR)
+ call mpas_log_write('****************************************************************', messageType=MPAS_LOG_CRIT)
+ end if
+
+ curr_time = mpas_get_clock_time(domain % clock, MPAS_NOW)
+ stop_time = mpas_get_clock_time(domain % clock, MPAS_STOP_TIME)
+ start_time = mpas_get_clock_time(domain % clock, MPAS_START_TIME)
+
+ do while (curr_time <= stop_time)
+
+ block_ptr => domain % blocklist
+ do while (associated(block_ptr))
+ call mpas_pool_get_subpool(block_ptr % structs, 'mesh', mesh)
+ call mpas_pool_get_subpool(block_ptr % structs, 'fg', fg)
+ call mpas_pool_get_subpool(block_ptr % structs, 'state', state)
+ call mpas_pool_get_subpool(block_ptr % structs, 'diag', diag)
+ call mpas_pool_get_subpool(block_ptr % structs, 'lbc_state', lbc_state)
+
+ call mpas_pool_get_array(state, 'xtime', xtime)
+ call mpas_pool_get_array(state, 'Time', Time)
+
+ call mpas_pool_get_dimension(block_ptr % dimensions, 'nCells', nCells)
+ call mpas_pool_get_dimension(block_ptr % dimensions, 'nEdges', nEdges)
+ call mpas_pool_get_dimension(block_ptr % dimensions, 'nVertLevels', nVertLevels)
+
+ call mpas_get_time(curr_time, dateTimeString=timeString)
+ xtime = timeString
+ time_since_start = curr_time - start_time
+ call mpas_get_TimeInterval(time_since_start, dt=dt)
+ Time = dt
+
+ call init_atm_case_lbc_ecmwf(timeString, block_ptr, mesh, nCells, nEdges, nVertLevels, fg, state, &
+ diag, lbc_state, block_ptr % dimensions, block_ptr % configs)
+
+ call mpas_get_time(start_time, dateTimeString=timeStart)
+ call init_atm_thompson_aerosols_lbc(timeString, timeStart, block_ptr, mesh, diag, lbc_state)
+
+ block_ptr => block_ptr % next
+ end do
+
+ call mpas_stream_mgr_write(stream_manager, streamID='lbc', ierr=ierr)
+ call mpas_stream_mgr_reset_alarms(stream_manager, streamID='lbc', direction=MPAS_STREAM_OUTPUT, ierr=ierr)
+
+ call mpas_advance_clock(domain % clock)
+ curr_time = mpas_get_clock_time(domain % clock, MPAS_NOW)
+
+ end do
+
+ call mpas_stream_mgr_reset_alarms(stream_manager, streamID='lbc', direction=MPAS_STREAM_OUTPUT, ierr=ierr)
+
else if (config_init_case == 13 ) then
call mpas_log_write(' CAM-MPAS grid ')
@@ -2652,9 +2801,13 @@ subroutine init_atm_case_mtn_wave(dminfo, mesh, nCells, nVertLevels, state, diag
end subroutine init_atm_case_mtn_wave
- subroutine init_atm_case_gfs(block, mesh, nCells, nEdges, nVertLevels, fg, state, diag, diag_physics, dims, configs)
+ subroutine init_atm_case(block, mesh, nCells, nEdges, nVertLevels, fg, state, diag, diag_physics, dims, configs, use_ecmwf)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Real-data test case using GFS data
+ ! Shared core for real-data initialization.
+ ! use_ecmwf=.false. → GFS/standard behaviour (3-level LS window, extrap_airtemp)
+ ! use_ecmwf=.true. → ECMWF model-level behaviour (physical-depth LS window)
+ ! Callers should use init_atm_case_gfs or init_atm_case_ecmwf;
+ ! do not call this subroutine directly from the dispatch block.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
use mpas_dmpar
@@ -2677,6 +2830,7 @@ subroutine init_atm_case_gfs(block, mesh, nCells, nEdges, nVertLevels, fg, state
type (mpas_pool_type), intent(inout):: diag_physics
type (mpas_pool_type), intent(inout):: dims
type (mpas_pool_type), intent(inout):: configs
+ logical, intent(in) :: use_ecmwf ! .true. selects physical-depth LS extrapolation
type (parallel_info), pointer :: parinfo
type (dm_info), pointer :: dminfo
@@ -2835,7 +2989,7 @@ subroutine init_atm_case_gfs(block, mesh, nCells, nEdges, nVertLevels, fg, state
real (kind=RKIND), dimension(:,:), pointer :: t_fg
real (kind=RKIND), dimension(:,:), pointer :: rh_fg
real (kind=RKIND), dimension(:,:), pointer :: sh_fg
- real (kind=RKIND), dimension(:,:), pointer :: gfs_z
+ real (kind=RKIND), dimension(:,:), pointer :: ght_fg
real (kind=RKIND), dimension(:,:), pointer :: p_fg
real (kind=RKIND), dimension(:,:), pointer :: st_fg
real (kind=RKIND), dimension(:,:), pointer :: sm_fg
@@ -2975,7 +3129,7 @@ subroutine init_atm_case_gfs(block, mesh, nCells, nEdges, nVertLevels, fg, state
call mpas_pool_get_array(fg, 't', t_fg)
call mpas_pool_get_array(fg, 'rh', rh_fg)
call mpas_pool_get_array(fg, 'sh', sh_fg)
- call mpas_pool_get_array(fg, 'gfs_z', gfs_z)
+ call mpas_pool_get_array(fg, 'ght_fg', ght_fg)
call mpas_pool_get_array(fg, 'p', p_fg)
call mpas_pool_get_dimension(dims, 'nVertLevels', nz1)
@@ -5023,91 +5177,145 @@ subroutine init_atm_case_gfs(block, mesh, nCells, nEdges, nVertLevels, fg, state
do iCell=1,nCells
- ! T
+ ! Temperature
+
sorted_arr(:,:) = -999.0
+
do k = 1, nfglevels_actual
sorted_arr(1,k) = z_fg(k,iCell)
if (vert_level(k) == 200100.0) sorted_arr(1,k) = 99999.0
sorted_arr(2,k) = t_fg(k,iCell)
end do
+
call mpas_quicksort(nfglevels_actual, sorted_arr)
+
do k = 1, nVertLevels
target_z = 0.5 * (zgrid(k,iCell) + zgrid(k+1,iCell))
- t(k,iCell) = vertical_interp(target_z, nfglevels_actual-1, &
- sorted_arr(:,1:nfglevels_actual-1), order=1, &
- extrap=extrap_airtemp, ierr=istatus)
+ if (use_ecmwf) then
+ t(k,iCell) = vertical_interp_ecmwf(target_z, nfglevels_actual-1, &
+ sorted_arr(:,1:nfglevels_actual-1), ierr=istatus)
+ else
+ t(k,iCell) = vertical_interp(target_z, nfglevels_actual-1, &
+ sorted_arr(:,1:nfglevels_actual-1), order=1, &
+ extrap=extrap_airtemp, ierr=istatus)
+ end if
+
if (istatus /= 0) then
write(errstring,'(a,i4,a,i10)') 'Error in interpolation of t(k,iCell) for k=', k, ', iCell=', iCell
call mpas_log_write('*****************************************************************', messageType=MPAS_LOG_ERR)
call mpas_log_write(trim(errstring), messageType=MPAS_LOG_ERR)
call mpas_log_write('*****************************************************************', messageType=MPAS_LOG_CRIT)
end if
+
end do
- ! RH
+ ! Relative humidity: if first-guess values are negative, set those values to zero before
+ ! vertical interpolation.
+
sorted_arr(:,:) = -999.0
relhum(:,iCell) = 0._RKIND
+
do k = 1, nfglevels_actual
sorted_arr(1,k) = z_fg(k,iCell)
if (vert_level(k) == 200100.0) sorted_arr(1,k) = 99999.0
sorted_arr(2,k) = rh_fg(k,iCell)
end do
+
call mpas_quicksort(nfglevels_actual, sorted_arr)
+
do k = nVertLevels, 1, -1
target_z = 0.5 * (zgrid(k,iCell) + zgrid(k+1,iCell))
- relhum(k,iCell) = vertical_interp(target_z, nfglevels_actual-1, &
- sorted_arr(:,1:nfglevels_actual-1), order=1, extrap=0)
+ if (use_ecmwf) then
+ relhum(k,iCell) = vertical_interp_ecmwf(target_z, nfglevels_actual-1, &
+ sorted_arr(:,1:nfglevels_actual-1))
+ else
+ 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
! SPECHUM: if first-guess values are negative, set those values to zero before
! vertical interpolation.
+
sorted_arr(:,:) = -999.0
spechum(:,iCell) = 0._RKIND
+
do k = 1, nfglevels_actual
sorted_arr(1,k) = z_fg(k,iCell)
if (vert_level(k) == 200100.0) sorted_arr(1,k) = 99999.0
sorted_arr(2,k) = max(0._RKIND,sh_fg(k,iCell))
end do
+
call mpas_quicksort(nfglevels_actual, sorted_arr)
+
do k = nVertLevels, 1, -1
target_z = 0.5 * (zgrid(k,iCell) + zgrid(k+1,iCell))
- spechum(k,iCell) = vertical_interp(target_z, nfglevels_actual-1, &
- sorted_arr(:,1:nfglevels_actual-1), order=1, extrap=0)
+ if (use_ecmwf) then
+ spechum(k,iCell) = vertical_interp_ecmwf(target_z, nfglevels_actual-1, &
+ sorted_arr(:,1:nfglevels_actual-1))
+ else
+ spechum(k,iCell) = vertical_interp(target_z, nfglevels_actual-1, &
+ sorted_arr(:,1:nfglevels_actual-1), order=1, extrap=2)
+ end if
+ ! Specific humidity must be non-negative; clamp after LS extrapolation
+ spechum(k,iCell) = max(0.0_RKIND, spechum(k,iCell))
end do
- ! GHT
+ ! Geopotential height
+
sorted_arr(:,:) = -999.0
+
do k = 1, nfglevels_actual
sorted_arr(1,k) = z_fg(k,iCell)
if (vert_level(k) == 200100.0) sorted_arr(1,k) = 99999.0
sorted_arr(2,k) = z_fg(k,iCell)
end do
+
call mpas_quicksort(nfglevels_actual, sorted_arr)
+
do k = 1, nVertLevels
target_z = 0.5 * (zgrid(k,iCell) + zgrid(k+1,iCell))
- gfs_z(k,iCell) = vertical_interp(target_z, nfglevels_actual-1, &
- sorted_arr(:,1:nfglevels_actual-1), order=1, extrap=1)
+ if (use_ecmwf) then
+ ght_fg(k,iCell) = vertical_interp_ecmwf(target_z, nfglevels_actual-1, &
+ 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)
+ end if
end do
- ! PRESSURE
+ ! PRESSURE: interpolate log(p) so exp() recovers a strictly positive pressure.
+ ! LS lapse-rate extrapolation of ln(p) is equivalent to fitting a local scale
+ ! height from the bottom three levels.
+
sorted_arr(:,:) = -999.0
+
do k = 1, nfglevels_actual
sorted_arr(1,k) = z_fg(k,iCell)
- if (vert_level(k) == 200100.0) then
+ if (vert_level(k) == 200100.0) then
sorted_arr(1,k) = 99999.0
sfc_k = k
end if
sorted_arr(2,k) = log(p_fg(k,iCell))
end do
+
call mpas_quicksort(nfglevels_actual, sorted_arr)
+
do k = 1, nVertLevels
target_z = 0.5 * (zgrid(k,iCell) + zgrid(k+1,iCell))
- pressure(k,iCell) = exp(vertical_interp(target_z, nfglevels_actual-1, &
- sorted_arr(:,1:nfglevels_actual-1), order=1, extrap=1))
+ if (use_ecmwf) then
+ pressure(k,iCell) = exp(vertical_interp_ecmwf(target_z, nfglevels_actual-1, &
+ 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))
+ end if
end do
end do
@@ -5115,19 +5323,26 @@ subroutine init_atm_case_gfs(block, mesh, nCells, nEdges, nVertLevels, fg, state
do iEdge=1,nEdges
- ! U
+ ! Wind: interpolate normal wind component to edges. We will reconstruct zonal and meridional winds later.
+
sorted_arr(:,:) = -999.0
+
do k=1,nfglevels_actual
sorted_arr(1,k) = 0.5 * (z_fg(k,cellsOnEdge(1,iEdge)) + z_fg(k,cellsOnEdge(2,iEdge)))
!NOSFC if (vert_level(k) == 200100.0) sorted_arr(1,k) = 0.5 * (fg % soilz % array(cellsOnEdge(1,iEdge)) + fg % soilz % array(cellsOnEdge(2,iEdge)))
if (vert_level(k) == 200100.0) sorted_arr(1,k) = 99999.0
sorted_arr(2,k) = u_fg(k,iEdge)
end do
+
call mpas_quicksort(nfglevels_actual, sorted_arr)
+
do k=1,nVertLevels
target_z = 0.25 * (zgrid(k,cellsOnEdge(1,iEdge)) + zgrid(k+1,cellsOnEdge(1,iEdge)) + zgrid(k,cellsOnEdge(2,iEdge)) + zgrid(k+1,cellsOnEdge(2,iEdge)))
-! u(k,iEdge) = vertical_interp(target_z, nfglevels_actual, sorted_arr, order=1, extrap=0)
- u(k,iEdge) = vertical_interp(target_z, nfglevels_actual-1, sorted_arr(:,1:nfglevels_actual-1), order=1, extrap=0)
+ 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)
+ end if
end do
end do
@@ -5167,7 +5382,11 @@ subroutine init_atm_case_gfs(block, mesh, nCells, nEdges, nVertLevels, fg, state
end do
call mpas_quicksort(nfglevels_actual, sorted_arr)
target_z = zgrid(1,iCell)
- psfc(iCell) = exp(vertical_interp(target_z, nfglevels_actual, sorted_arr, order=1, extrap=1))
+ 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))
+ end if
end do
@@ -5414,9 +5633,61 @@ subroutine init_atm_case_gfs(block, mesh, nCells, nEdges, nVertLevels, fg, state
end if ! config_met_interp
+ end subroutine init_atm_case
+
+
+ !------------------------------------------------------------------------------
+ ! Routine: init_atm_case_gfs
+ !
+ ! Purpose: thin wrapper around init_atm_case to select the GFS/standard pathway
+ ! (i.e., no ECMWF-specific handling of interpolation/extrapolation).
+ !
+ ! Author: Rubaiat Islam, NCAR/MMM, 2026-05-01
+ !------------------------------------------------------------------------------
+ subroutine init_atm_case_gfs(block, mesh, nCells, nEdges, nVertLevels, fg, state, diag, diag_physics, dims, configs)
+ implicit none
+ type (block_type), intent(inout), target :: block
+ type (mpas_pool_type), intent(inout) :: mesh
+ integer, intent(in) :: nCells
+ integer, intent(in) :: nEdges
+ integer, intent(in) :: nVertLevels
+ type (mpas_pool_type), intent(inout) :: fg
+ type (mpas_pool_type), intent(inout) :: state
+ type (mpas_pool_type), intent(inout) :: diag
+ type (mpas_pool_type), intent(inout) :: diag_physics
+ type (mpas_pool_type), intent(inout) :: dims
+ type (mpas_pool_type), intent(inout) :: configs
+ call init_atm_case(block, mesh, nCells, nEdges, nVertLevels, fg, state, diag, &
+ diag_physics, dims, configs, use_ecmwf=.false.)
end subroutine init_atm_case_gfs
+ !-----------------------------------------------------------------------
+ ! Routine: init_atm_case_ecmwf
+ !
+ ! Purpose: thin wrapper around init_atm_case to select the ECMWF-specific pathway
+ ! (i.e., use of physical-depth-based vertical interpolation/extrapolation).
+ !
+ ! Author: Rubaiat Islam, NCAR/MMM, 2026-05-01
+ !-----------------------------------------------------------------------
+ subroutine init_atm_case_ecmwf(block, mesh, nCells, nEdges, nVertLevels, fg, state, diag, diag_physics, dims, configs)
+ implicit none
+ type (block_type), intent(inout), target :: block
+ type (mpas_pool_type), intent(inout) :: mesh
+ integer, intent(in) :: nCells
+ integer, intent(in) :: nEdges
+ integer, intent(in) :: nVertLevels
+ type (mpas_pool_type), intent(inout) :: fg
+ type (mpas_pool_type), intent(inout) :: state
+ type (mpas_pool_type), intent(inout) :: diag
+ type (mpas_pool_type), intent(inout) :: diag_physics
+ type (mpas_pool_type), intent(inout) :: dims
+ type (mpas_pool_type), intent(inout) :: configs
+ call init_atm_case(block, mesh, nCells, nEdges, nVertLevels, fg, state, diag, &
+ diag_physics, dims, configs, use_ecmwf=.true.)
+ end subroutine init_atm_case_ecmwf
+
+
!-----------------------------------------------------------------------
! routine init_atm_case_lbc
!
@@ -5431,7 +5702,7 @@ end subroutine init_atm_case_gfs
!> fields that are needed as model lateral boundary conditions.
!
!-----------------------------------------------------------------------
- subroutine init_atm_case_lbc(timestamp, block, mesh, nCells, nEdges, nVertLevels, fg, state, diag, lbc_state, dims, configs)
+ subroutine init_atm_lbc_core(timestamp, block, mesh, nCells, nEdges, nVertLevels, fg, state, diag, lbc_state, dims, configs, use_ecmwf)
use mpas_dmpar, only : mpas_dmpar_min_real, mpas_dmpar_max_real
use init_atm_read_met, only : met_data, read_met_init, read_met_close, read_next_met_field
@@ -5453,6 +5724,7 @@ subroutine init_atm_case_lbc(timestamp, block, mesh, nCells, nEdges, nVertLevels
type (mpas_pool_type), intent(inout) :: lbc_state
type (mpas_pool_type), intent(inout):: dims
type (mpas_pool_type), intent(inout):: configs
+ logical, intent(in) :: use_ecmwf ! .true. selects physical-depth LS extrapolation
type (dm_info), pointer :: dminfo
@@ -5932,9 +6204,14 @@ subroutine init_atm_case_lbc(timestamp, block, mesh, nCells, nEdges, nVertLevels
call mpas_quicksort(nfglevels_actual, sorted_arr)
do k = 1, nVertLevels
target_z = 0.5 * (zgrid(k,iCell) + zgrid(k+1,iCell))
- t(k,iCell) = vertical_interp(target_z, nfglevels_actual-1, &
- sorted_arr(:,1:nfglevels_actual-1), order=1, &
- extrap=extrap_airtemp, ierr=istatus)
+ if (use_ecmwf) then
+ t(k,iCell) = vertical_interp_ecmwf(target_z, nfglevels_actual-1, &
+ sorted_arr(:,1:nfglevels_actual-1), ierr=istatus)
+ else
+ t(k,iCell) = vertical_interp(target_z, nfglevels_actual-1, &
+ sorted_arr(:,1:nfglevels_actual-1), order=1, &
+ extrap=extrap_airtemp, ierr=istatus)
+ end if
if (istatus /= 0) then
write(errstring,'(a,i4,a,i10)') 'Error in interpolation of t(k,iCell) for k=', k, ', iCell=', iCell
call mpas_log_write('*****************************************************************', messageType=MPAS_LOG_ERR)
@@ -5955,8 +6232,15 @@ subroutine init_atm_case_lbc(timestamp, block, mesh, nCells, nEdges, nVertLevels
call mpas_quicksort(nfglevels_actual, sorted_arr)
do k = nVertLevels, 1, -1
target_z = 0.5 * (zgrid(k,iCell) + zgrid(k+1,iCell))
- relhum(k,iCell) = vertical_interp(target_z, nfglevels_actual-1, &
- sorted_arr(:,1:nfglevels_actual-1), order=1, extrap=0)
+ if (use_ecmwf) then
+ relhum(k,iCell) = vertical_interp_ecmwf(target_z, nfglevels_actual-1, &
+ sorted_arr(:,1:nfglevels_actual-1))
+ else
+ 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
@@ -5972,8 +6256,15 @@ subroutine init_atm_case_lbc(timestamp, block, mesh, nCells, nEdges, nVertLevels
call mpas_quicksort(nfglevels_actual, sorted_arr)
do k = nVertLevels, 1, -1
target_z = 0.5 * (zgrid(k,iCell) + zgrid(k+1,iCell))
- spechum(k,iCell) = vertical_interp(target_z, nfglevels_actual-1, &
- sorted_arr(:,1:nfglevels_actual-1), order=1, extrap=0)
+ if (use_ecmwf) then
+ spechum(k,iCell) = vertical_interp_ecmwf(target_z, nfglevels_actual-1, &
+ sorted_arr(:,1:nfglevels_actual-1))
+ else
+ spechum(k,iCell) = vertical_interp(target_z, nfglevels_actual-1, &
+ sorted_arr(:,1:nfglevels_actual-1), order=1, extrap=2)
+ end if
+ ! Specific humidity must be non-negative; clamp after LS extrapolation
+ spechum(k,iCell) = max(0.0_RKIND, spechum(k,iCell))
end do
@@ -5991,8 +6282,13 @@ subroutine init_atm_case_lbc(timestamp, block, mesh, nCells, nEdges, nVertLevels
call mpas_quicksort(nfglevels_actual, sorted_arr)
do k = 1, nVertLevels
target_z = 0.5 * (zgrid(k,iCell) + zgrid(k+1,iCell))
- pressure(k,iCell) = exp(vertical_interp(target_z, nfglevels_actual-1, &
- sorted_arr(:,1:nfglevels_actual-1), order=1, extrap=1))
+ if (use_ecmwf) then
+ pressure(k,iCell) = exp(vertical_interp_ecmwf(target_z, nfglevels_actual-1, &
+ 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))
+ end if
end do
end do
@@ -6011,7 +6307,11 @@ subroutine init_atm_case_lbc(timestamp, block, mesh, nCells, nEdges, nVertLevels
do k=1,nVertLevels
target_z = 0.25 * (zgrid(k,cellsOnEdge(1,iEdge)) + zgrid(k+1,cellsOnEdge(1,iEdge)) &
+ zgrid(k,cellsOnEdge(2,iEdge)) + zgrid(k+1,cellsOnEdge(2,iEdge)))
- u(k,iEdge) = vertical_interp(target_z, nfglevels_actual-1, sorted_arr(:,1:nfglevels_actual-1), order=1, extrap=1)
+ 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)
+ end if
end do
end do
@@ -6215,8 +6515,64 @@ subroutine init_atm_case_lbc(timestamp, block, mesh, nCells, nEdges, nVertLevels
end do
end do
+ end subroutine init_atm_lbc_core
+
+
+ !-----------------------------------------------------------------------
+ ! Routine: init_atm_case_lbc
+ !
+ ! Purpose: This is a thin wrapper around init_atm_lbc_core for the standard
+ ! LBC pathway, which does not use ECMWF/ERA5/IFS-specific vertical
+ ! interpolation.
+ !
+ ! Author: Rubaiat Islam, NCAR/MMM, 2026-05-01
+ !-----------------------------------------------------------------------
+ subroutine init_atm_case_lbc(timestamp, block, mesh, nCells, nEdges, nVertLevels, fg, state, diag, lbc_state, dims, configs)
+ implicit none
+ character(len=*), intent(in) :: timestamp
+ type (block_type), intent(inout), target :: block
+ type (mpas_pool_type), intent(inout) :: mesh
+ integer, intent(in) :: nCells
+ integer, intent(in) :: nEdges
+ integer, intent(in) :: nVertLevels
+ type (mpas_pool_type), intent(inout) :: fg
+ type (mpas_pool_type), intent(inout) :: state
+ type (mpas_pool_type), intent(inout) :: diag
+ type (mpas_pool_type), intent(inout) :: lbc_state
+ type (mpas_pool_type), intent(inout) :: dims
+ type (mpas_pool_type), intent(inout) :: configs
+ call init_atm_lbc_core(timestamp, block, mesh, nCells, nEdges, nVertLevels, fg, state, &
+ diag, lbc_state, dims, configs, use_ecmwf=.false.)
end subroutine init_atm_case_lbc
+
+ !-----------------------------------------------------------------------
+ ! Routine: init_atm_case_lbc_ecmwf
+ !
+ ! Purpose: This is a thin wrapper around init_atm_lbc_core for the ERA5/IFS
+ ! LBC pathway, which uses ECMWF-specific vertical interpolation.
+ !
+ ! Author: Rubaiat Islam, NCAR/MMM, 2026-05-01
+ !-----------------------------------------------------------------------
+ subroutine init_atm_case_lbc_ecmwf(timestamp, block, mesh, nCells, nEdges, nVertLevels, fg, state, diag, lbc_state, dims, configs)
+ implicit none
+ character(len=*), intent(in) :: timestamp
+ type (block_type), intent(inout), target :: block
+ type (mpas_pool_type), intent(inout) :: mesh
+ integer, intent(in) :: nCells
+ integer, intent(in) :: nEdges
+ integer, intent(in) :: nVertLevels
+ type (mpas_pool_type), intent(inout) :: fg
+ type (mpas_pool_type), intent(inout) :: state
+ type (mpas_pool_type), intent(inout) :: diag
+ type (mpas_pool_type), intent(inout) :: lbc_state
+ type (mpas_pool_type), intent(inout) :: dims
+ type (mpas_pool_type), intent(inout) :: configs
+ call init_atm_lbc_core(timestamp, block, mesh, nCells, nEdges, nVertLevels, fg, state, &
+ diag, lbc_state, dims, configs, use_ecmwf=.true.)
+ end subroutine init_atm_case_lbc_ecmwf
+
+
!---------------------
subroutine init_atm_case_les(dminfo, mesh, fg, nCells, nVertLevels, state, diag, test_case, configs)
diff --git a/src/core_init_atmosphere/mpas_init_atm_core_interface.F b/src/core_init_atmosphere/mpas_init_atm_core_interface.F
index f277a4a72f..45a5b42433 100644
--- a/src/core_init_atmosphere/mpas_init_atm_core_interface.F
+++ b/src/core_init_atmosphere/mpas_init_atm_core_interface.F
@@ -195,7 +195,7 @@ function init_atm_setup_packages(configs, streamInfo, packages, iocontext) resul
sfc_update = .false.
end if
- if (config_init_case == 9) then
+ 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)
@@ -258,6 +258,40 @@ function init_atm_setup_packages(configs, streamInfo, packages, iocontext) resul
initial_conds = .false. ! Also, turn off the initial_conds package to avoid writing the IC "output" stream
+ else if (config_init_case == 11) then
+ gwd_stage_in = config_native_gwd_static .and. &
+ (.not. config_static_interp)
+ gwd_stage_out = config_native_gwd_static
+ gwd_gsl_stage_out = config_native_gwd_gsl_static
+ vertical_stage_in = config_vertical_grid .and. &
+ (.not. config_native_gwd_static) .and. &
+ (.not. config_static_interp)
+ vertical_stage_out = config_vertical_grid
+ met_stage_in = config_met_interp .and. &
+ (.not. config_native_gwd_static) .and. &
+ (.not. config_static_interp) .and. &
+ (.not. config_vertical_grid)
+ met_stage_out = config_met_interp
+
+ mp_thompson_aers_in = .false.
+ inquire(file="QNWFA_QNIFA_SIGMA_MONTHLY.dat",exist=lexist)
+ if((lexist .and. met_stage_out) .or. (lexist .and. met_stage_in)) mp_thompson_aers_in = .true.
+
+ else if (config_init_case == 12) then
+ gwd_stage_in = .false.
+ gwd_stage_out = .false.
+ gwd_gsl_stage_out = .false.
+ vertical_stage_in = .false.
+ vertical_stage_out = .false.
+ met_stage_in = .true.
+ met_stage_out = .true.
+
+ mp_thompson_aers_in = .false.
+ inquire(file="QNWFA_QNIFA_SIGMA_MONTHLY.dat",exist=lexist)
+ if((lexist .and. met_stage_out) .or. (lexist .and. met_stage_in)) mp_thompson_aers_in = .true.
+
+ initial_conds = .false.
+
else if (config_init_case == 13) then
gwd_stage_in = .false.
gwd_stage_out = .false.
@@ -285,7 +319,9 @@ function init_atm_setup_packages(configs, streamInfo, packages, iocontext) resul
call mpas_pool_get_package(packages, 'first_guess_fieldActive', first_guess_field)
first_guess_field = .false.
- if ((config_init_case == 7 .and. config_met_interp) .or. config_init_case == 9) then
+
+ if ((config_init_case == 7 .and. config_met_interp) .or. config_init_case == 9 .or. &
+ (config_init_case == 11 .and. config_met_interp) .or. config_init_case == 12) then
first_guess_field = .true.
end if
diff --git a/src/core_init_atmosphere/mpas_init_atm_vinterp.F b/src/core_init_atmosphere/mpas_init_atm_vinterp.F
index 164cfac546..5f198a7d90 100644
--- a/src/core_init_atmosphere/mpas_init_atm_vinterp.F
+++ b/src/core_init_atmosphere/mpas_init_atm_vinterp.F
@@ -38,9 +38,10 @@ 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
+ logical :: found
integer :: interp_order, extrap_type
real (kind=RKIND) :: surface, sealevel
@@ -81,7 +82,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 +95,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
@@ -103,20 +106,212 @@ real (kind=RKIND) function vertical_interp(target_z, nz, zf, order, extrap, surf
!
! No extrapolation required
!
+
+ found = .false.
+
do k=1,nz-1
if (target_z >= zf(1,k) .and. target_z < zf(1,k+1)) then
lm = k
lp = k+1
wm = (zf(1,k+1) - target_z) / (zf(1,k+1) - zf(1,k))
wp = (target_z - zf(1,k)) / (zf(1,k+1) - zf(1,k))
+ found = .true.
exit
end if
end do
+ 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 = zf(2,lm)
+ return
+ end if
+
vertical_interp = wm*zf(2,lm) + wp*zf(2,lp)
return
end function vertical_interp
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Routine: ls_lapse_slope
+ !
+ ! Purpose: Compute an ordinary least-squares slope over levels k_start..k_end
+ ! of zf. x is shifted by the first selected level before accumulating
+ ! to prevent catastrophic cancellation in n*sum_x2 - sum_x**2 when z
+ ! values are large relative to their spacing (e.g. z~50 km, dz~1 km
+ ! loses ~3 digits in single precision without the shift). After
+ ! centering, the denominator is a sum of squares and therefore always
+ ! non-negative; a relative epsilon guard catches degenerate z-spacing
+ ! before a divide-by-zero or Inf/NaN can propagate. Returns 0 in all
+ ! degenerate cases (constant extrapolation fallback).
+ !
+ ! Author: Rubaiat Islam, NCAR/MMM, 2026-04-17
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ 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) :: x_ref, dx
+ real (kind=RKIND) :: sum_x, sum_y, sum_x2, sum_xy, n_real, denom
+
+ n_pts = k_end - k_start + 1
+
+ if (n_pts <= 1) then
+ ls_lapse_slope = 0.0_RKIND
+ return
+ end if
+
+ x_ref = zf(1, k_start) ! shift origin: dx values are O(dz), not O(z)
+
+ 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
+ dx = zf(1,k) - x_ref
+ sum_x = sum_x + dx
+ sum_y = sum_y + zf(2,k)
+ sum_x2 = sum_x2 + dx**2
+ sum_xy = sum_xy + dx * zf(2,k)
+ end do
+
+ n_real = real(n_pts, RKIND)
+ denom = n_real * sum_x2 - sum_x * sum_x
+
+ ! denom >= 0 after centering; it is ~0 only when the selected levels
+ ! share the same (or nearly the same) vertical coordinate.
+ if (denom > epsilon(1.0_RKIND) * n_real * sum_x2) then
+ ls_lapse_slope = (n_real*sum_xy - sum_x*sum_y) / denom
+ else
+ ls_lapse_slope = 0.0_RKIND
+ end if
+
+ end function ls_lapse_slope
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Routine: vertical_interp_ecmwf
+ !
+ ! Purpose: ECMWF model-level aware vertical interpolation.
+ !
+ ! Interior targets: identical linear interpolation to vertical_interp.
+ !
+ ! Below/above-data extrapolation: OLS slope fitted over ALL levels whose
+ ! height lies within a fixed physical depth from the boundary:
+ !
+ ! Below data: levels k where zf(1,k) - zf(1,1) <= DEPTH_BOT (3000 m)
+ ! Above data: levels k where zf(1,nz) - zf(1,k) <= DEPTH_TOP (5000 m)
+ !
+ ! With ERA5/IFS 137-level data this selects ~20-25 levels near the surface
+ ! (spanning the well-mixed boundary layer and lower free troposphere) and
+ ! ~5-10 levels near the model top, compared with only 3 levels for the
+ ! GFS-oriented vertical_interp. The result is a slope that represents the
+ ! actual lower-tropospheric lapse rate rather than the shallow skin layer
+ ! captured by the first three fine-resolution levels (~100 m depth).
+ !
+ ! N is clamped to [N_MIN=3, N_MAX=30] for numerical stability.
+ ! The existing ls_lapse_slope helper (x-centred OLS, epsilon-guarded
+ ! denominator) is reused unchanged.
+ !
+ ! Author: Rubaiat Islam, NCAR/MMM, 2026-04-17
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ real (kind=RKIND) function vertical_interp_ecmwf(target_z, nz, zf, ierr)
+
+ implicit none
+
+ real (kind=RKIND), intent(in) :: target_z
+ integer, intent(in) :: nz
+ real (kind=RKIND), dimension(2,nz), intent(in) :: zf
+ integer, intent(out), optional :: ierr
+
+ real (kind=RKIND), parameter :: DEPTH_BOT = 3000.0_RKIND ! m: below-boundary window
+ real (kind=RKIND), parameter :: DEPTH_TOP = 5000.0_RKIND ! m: above-boundary window
+ integer, parameter :: N_MIN = 3 ! minimum OLS levels (stability)
+ integer, parameter :: N_MAX = 30 ! maximum OLS levels (cap distant levels)
+
+ integer :: k, lm, lp, n_depth
+ real (kind=RKIND) :: wm, wp, slope
+ logical :: found
+
+ if (present(ierr)) ierr = 0
+
+ !
+ ! Below lowest level: physical-depth window from the bottom
+ !
+
+ if (target_z < zf(1,1)) then
+ n_depth = 1
+ do k = 2, nz
+ if (zf(1,k) - zf(1,1) <= DEPTH_BOT) then
+ n_depth = k
+ else
+ exit ! zf is sorted ascending; levels beyond are out of window
+ end if
+ end do
+ n_depth = min(max(n_depth, N_MIN), min(N_MAX, nz))
+ slope = ls_lapse_slope(zf, nz, 1, n_depth)
+ vertical_interp_ecmwf = zf(2,1) + slope * (target_z - zf(1,1))
+ return
+ end if
+
+ !
+ ! Above highest level: physical-depth window from the top
+ !
+
+ if (target_z >= zf(1,nz)) then
+ n_depth = 1
+ do k = nz-1, 1, -1
+ if (zf(1,nz) - zf(1,k) <= DEPTH_TOP) then
+ n_depth = nz - k + 1
+ else
+ exit ! zf is sorted ascending; levels beyond are out of window
+ end if
+ end do
+ n_depth = min(max(n_depth, N_MIN), min(N_MAX, nz))
+ slope = ls_lapse_slope(zf, nz, nz - n_depth + 1, nz)
+ vertical_interp_ecmwf = zf(2,nz) + slope * (target_z - zf(1,nz))
+ return
+ end if
+
+ !
+ ! Interior: linear interpolation
+ !
+
+ found = .false.
+
+ do k = 1, nz-1
+ if (target_z >= zf(1,k) .and. target_z < zf(1,k+1)) then
+ lm = k
+ lp = k+1
+ wm = (zf(1,k+1) - target_z) / (zf(1,k+1) - zf(1,k))
+ wp = (target_z - zf(1,k)) / (zf(1,k+1) - zf(1,k))
+ found = .true.
+ exit
+ end if
+ end do
+
+ 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
+ end if
+
+ vertical_interp_ecmwf = wm * zf(2,lm) + wp * zf(2,lp)
+
+ return
+
+ end function vertical_interp_ecmwf
+
end module init_atm_vinterp