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 1/3] 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 + From a2fd267e234f23e74d97edcc103562c7426be6c5 Mon Sep 17 00:00:00 2001 From: Rubaiat Islam Date: Wed, 29 Apr 2026 01:30:53 -0600 Subject: [PATCH 2/3] bugfix vinterp --- .../mpas_init_atm_vinterp.F | 65 ++++++++++++------- 1 file changed, 42 insertions(+), 23 deletions(-) diff --git a/src/core_init_atmosphere/mpas_init_atm_vinterp.F b/src/core_init_atmosphere/mpas_init_atm_vinterp.F index e39a496066..72bbd20d7b 100644 --- a/src/core_init_atmosphere/mpas_init_atm_vinterp.F +++ b/src/core_init_atmosphere/mpas_init_atm_vinterp.F @@ -24,26 +24,30 @@ module init_atm_vinterp !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 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 ! error code: 0 = success, 1 = invalid extrapolation type - - integer :: k, lm, lp + integer, intent(out), optional :: ierr + + integer :: k, lm, lp, n_lapse real (kind=RKIND) :: wm, wp real (kind=RKIND) :: slope - real (kind=RKIND) :: lapse_rate + real (kind=RKIND) :: sum_x, sum_y, sum_x2, sum_xy, n_real integer :: interp_order, extrap_type real (kind=RKIND) :: surface, sealevel - + if (present(ierr)) ierr = 0 + if (present(order)) then interp_order = order else @@ -68,11 +72,6 @@ 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 ! @@ -83,12 +82,22 @@ 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 - ! 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 + 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) + 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_lapse, RKIND) + if (n_lapse > 1) then + slope = (n_real*sum_xy - sum_x*sum_y) / (n_real*sum_x2 - sum_x**2) + else + slope = 0.0_RKIND + end if + vertical_interp = zf(2,1) + slope * (target_z - zf(1,1)) end if return end if @@ -99,12 +108,22 @@ 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 - ! 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 + 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) + 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_lapse, RKIND) + if (n_lapse > 1) then + slope = (n_real*sum_xy - sum_x*sum_y) / (n_real*sum_x2 - sum_x**2) + else + slope = 0.0_RKIND + end if + vertical_interp = zf(2,nz) + slope * (target_z - zf(1,nz)) end if return end if From 924b456c20e146d47b227a8140816dc8985260ec Mon Sep 17 00:00:00 2001 From: Rubaiat Islam Date: Wed, 29 Apr 2026 02:36:24 -0600 Subject: [PATCH 3/3] Add function for least square slope estimation --- .../mpas_init_atm_vinterp.F | 69 +++++++++++-------- 1 file changed, 40 insertions(+), 29 deletions(-) diff --git a/src/core_init_atmosphere/mpas_init_atm_vinterp.F b/src/core_init_atmosphere/mpas_init_atm_vinterp.F index 72bbd20d7b..a88505cb4a 100644 --- a/src/core_init_atmosphere/mpas_init_atm_vinterp.F +++ b/src/core_init_atmosphere/mpas_init_atm_vinterp.F @@ -41,7 +41,6 @@ real (kind=RKIND) function vertical_interp(target_z, nz, zf, order, extrap, surf integer :: k, lm, lp, n_lapse real (kind=RKIND) :: wm, wp real (kind=RKIND) :: slope - real (kind=RKIND) :: sum_x, sum_y, sum_x2, sum_xy, n_real integer :: interp_order, extrap_type real (kind=RKIND) :: surface, sealevel @@ -83,20 +82,7 @@ real (kind=RKIND) function vertical_interp(target_z, nz, zf, order, extrap, surf vertical_interp = zf(2,1) + slope * (target_z - zf(1,1)) else if (extrap_type == 2) then 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) - 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_lapse, RKIND) - if (n_lapse > 1) then - slope = (n_real*sum_xy - sum_x*sum_y) / (n_real*sum_x2 - sum_x**2) - else - slope = 0.0_RKIND - end if + slope = ls_lapse_slope(zf, nz, 1, n_lapse) vertical_interp = zf(2,1) + slope * (target_z - zf(1,1)) end if return @@ -109,20 +95,7 @@ real (kind=RKIND) function vertical_interp(target_z, nz, zf, order, extrap, surf vertical_interp = zf(2,nz) + slope * (target_z - zf(1,nz)) else if (extrap_type == 2) then 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) - 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_lapse, RKIND) - if (n_lapse > 1) then - slope = (n_real*sum_xy - sum_x*sum_y) / (n_real*sum_x2 - sum_x**2) - else - slope = 0.0_RKIND - end if + 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 @@ -148,5 +121,43 @@ real (kind=RKIND) function vertical_interp(target_z, nz, zf, order, extrap, surf end function vertical_interp + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! Name: ls_lapse_slope + ! + ! Purpose: Compute an ordinary least-squares slope over levels k_start..k_end of + ! zf. Returns 0 when only one level is supplied (constant fallback). + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + 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) :: 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 + end if + + end function ls_lapse_slope + end module init_atm_vinterp