Skip to content

Consolidating multiple definitions of vertical_interp into init_atm_vinterp #1429

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
mgduda merged 1 commit into from
Mar 25, 2026
+16 −99
Merged

Consolidating multiple definitions of vertical_interp into init_atm_vinterp #1429

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
96 changes: 1 addition & 95 deletions src/core_init_atmosphere/mpas_init_atm_cases.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +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 mpas_timer
use mpas_init_atm_static
use mpas_init_atm_surface
Expand Down Expand Up @@ -7595,101 +7596,6 @@ integer function nearest_edge(target_lat, target_lon, &
end function nearest_edge


real (kind=RKIND) function vertical_interp(target_z, nz, zf, order, extrap, surface_val, sealev_val, ierr)

implicit none

real (kind=RKIND), intent(in) :: target_z
integer, intent(in) :: nz
real (kind=RKIND), dimension(2,nz), intent(in) :: zf ! zf(1,:) is column of vertical coordinate values, zf(2,:) is column of field values
integer, intent(in), optional :: order
integer, intent(in), optional :: extrap ! can take values 0 = constant, 1 = linear (default), 2 = lapse-rate
real (kind=RKIND), intent(in), optional :: surface_val
real (kind=RKIND), intent(in), optional :: sealev_val
integer, intent(out), optional :: ierr

integer :: k, lm, lp
real (kind=RKIND) :: wm, wp
real (kind=RKIND) :: slope

integer :: interp_order, extrap_type
real (kind=RKIND) :: surface, sealevel

if (present(ierr)) ierr = 0

if (present(order)) then
interp_order = order
else
interp_order = 2
end if

if (present(extrap)) then
extrap_type = extrap
else
extrap_type = 1
end if

if (present(surface_val)) then
surface = surface_val
else
surface = 200100.0
end if

if (present(sealev_val)) then
sealevel = sealev_val
else
sealevel = 201300.0
end if

!
! Extrapolation required
!
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
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
end if
return
end if


!
! No extrapolation required
!
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))
exit
end if
end do

vertical_interp = wm*zf(2,lm) + wp*zf(2,lp)

return

end function vertical_interp


!----------------------------------------------------------------------------------------------------------

real (kind=RKIND) function env_qv( temperature, pressure, rh_max )
Expand Down
19 changes: 15 additions & 4 deletions src/core_init_atmosphere/mpas_init_atm_vinterp.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,26 +22,31 @@ 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)

use mpas_log, only : mpas_log_write
use mpas_derived_types, only : MPAS_LOG_ERR

implicit none

real (kind=RKIND), intent(in) :: target_z
integer, intent(in) :: nz
real (kind=RKIND), dimension(2,nz), intent(in) :: zf ! zf(1,:) is column of vertical coordinate values, zf(2,:) is column of field values
integer, intent(in), optional :: order
integer, intent(in), optional :: extrap
integer, intent(in), optional :: extrap ! can take values 0 = constant, 1 = linear (default), 2 = lapse-rate
real (kind=RKIND), intent(in), optional :: surface_val
real (kind=RKIND), intent(in), optional :: sealev_val

integer, intent(out), optional :: ierr

integer :: k, lm, lp
real (kind=RKIND) :: wm, wp
real (kind=RKIND) :: slope

integer :: interp_order, extrap_type
real (kind=RKIND) :: surface, sealevel


if (present(ierr)) ierr = 0

if (present(order)) then
interp_order = order
else
Expand Down Expand Up @@ -75,6 +80,8 @@ 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
vertical_interp = zf(2,1) - (target_z - zf(1,1))*0.0065
end if
return
end if
Expand All @@ -84,6 +91,10 @@ 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
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
end if
return
end if
Expand Down