From d635aed972f498545ce6723c7b9eb2adfa437baa Mon Sep 17 00:00:00 2001 From: Bill Skamarock Date: Tue, 30 Sep 2025 14:54:23 -0600 Subject: [PATCH 1/4] Variable epssm(zeta) introduction: epssm is an offcentering coefficient for the vertically (column-wise) semi-implicit acoustic/gravity-wave integration. Presently the coefficent is a single constant. We now allow this coefficient to vary with model and interface levels as a function of height zeta between two values with a transition layer between the two. This commit, to the Registry only, adds namelist control variables for specifying the lower and upper values of epssm and the lower and upper heights of the transition region for these values. The commit also includes Registry-defined arrays storing the values (1 +- epssm(z))/2 at the levels and interfaces needed in the solution procedure. --- src/core_atmosphere/Registry.xml | 37 ++++++++++++++++++++++++++++++++ 1 file changed, 37 insertions(+) diff --git a/src/core_atmosphere/Registry.xml b/src/core_atmosphere/Registry.xml index 4281c40bba..fc599be2be 100644 --- a/src/core_atmosphere/Registry.xml +++ b/src/core_atmosphere/Registry.xml @@ -292,6 +292,27 @@ units="-" description="Number of layers in which to apply Rayleigh damping on horizontal velocity at top of model; damping linearly ramps to zero by layer number from the top" possible_values="Positive integer values"/> + + + + + + + + + @@ -495,6 +516,10 @@ + + + + #ifdef MPAS_CAM_DYCORE @@ -1484,6 +1509,18 @@ + + + + + + + + From 14ae8f8a37cfd72202e0502bd68c038868bf891e Mon Sep 17 00:00:00 2001 From: Bill Skamarock Date: Wed, 1 Oct 2025 13:48:53 -0600 Subject: [PATCH 2/4] Variable epssm(zeta) introduction, step (2): Added code into mpas_atm_core.F to initialize the Registry-defined arrays storing the values (1 +- epssm(z))/2 at the levels and interfaces needed in the solution procedure. The initialization occurs at model integration start up. The vertically varying epssm is NOT active at this point; the simulation still uses the constant value from config_epssm. --- src/core_atmosphere/mpas_atm_core.F | 72 +++++++++++++++++++++++++++++ 1 file changed, 72 insertions(+) diff --git a/src/core_atmosphere/mpas_atm_core.F b/src/core_atmosphere/mpas_atm_core.F index f7d04a1f0c..03cc02d8f6 100644 --- a/src/core_atmosphere/mpas_atm_core.F +++ b/src/core_atmosphere/mpas_atm_core.F @@ -1253,6 +1253,11 @@ subroutine atm_compute_damping_coefs(mesh, configs) real (kind=RKIND), dimension(:), pointer :: meshDensity real (kind=RKIND) :: dx_scale_power + real (kind=RKIND), dimension(:), pointer :: rdzw, etp, etm, ewp, ewm + real (kind=RKIND), pointer :: max_coeff, min_coeff, transition_lower_bound, transition_upper_bound + real (kind=RKIND) :: transition_width + real (kind=RKIND), allocatable, dimension(:) :: height_u_levels, epssm_coeff_u, height_w_levels, epssm_coeff_w + m1 = -1.0 pii = acos(m1) @@ -1279,6 +1284,73 @@ subroutine atm_compute_damping_coefs(mesh, configs) end do end do + call mpas_pool_get_array(mesh, 'rdzw', rdzw) + call mpas_pool_get_array(mesh, 'etp', etp) + call mpas_pool_get_array(mesh, 'etm', etm) + call mpas_pool_get_array(mesh, 'ewp', ewp) + call mpas_pool_get_array(mesh, 'ewm', ewm) + call mpas_pool_get_config(configs, 'config_epssm_minimum', min_coeff) + call mpas_pool_get_config(configs, 'config_epssm_maximum', max_coeff) + call mpas_pool_get_config(configs, 'config_epssm_transition_bottom_z', transition_lower_bound) + call mpas_pool_get_config(configs, 'config_epssm_transition_top_z', transition_upper_bound) + + allocate(height_u_levels(nVertLevels)) + allocate(epssm_coeff_u(nVertLevels)) + allocate(height_w_levels(nVertLevels+1)) + allocate(epssm_coeff_w(nVertLevels+1)) + + transition_width = transition_upper_bound - transition_lower_bound + +! initialization for heights of coordinate at u and w levels +! These are the heights of the computational coordinate zeta + + height_w_levels(:) = 0.0_RKIND + do k =1, nVertLevels + height_w_levels(k+1) = height_w_levels(k) + 1.0_RKIND/rdzw(k) + height_u_levels(k) = 0.5*(height_w_levels(k+1) + height_w_levels(k)) + enddo + +! Height dependent values of epssm; profiles stored in etp, etm, ewp, and ewm, + + call mpas_log_write(' setting epssm coefficients ') + call mpas_log_write(' minimum epssm: $r ',realArgs=(/min_coeff/)) + call mpas_log_write(' maximum epssm: $r ',realArgs=(/max_coeff/)) + call mpas_log_write(' transition lower bound (m): $r ',realArgs=(/transition_lower_bound/)) + call mpas_log_write(' transition upper bound (m): $r ',realArgs=(/transition_upper_bound/)) + call mpas_log_write(' ') + + do k = 1,nVertLevels + if(height_u_levels(k).le.transition_lower_bound) then + epssm_coeff_u(k) = min_coeff + else if(height_u_levels(k).ge.transition_upper_bound) then + epssm_coeff_u(k) = max_coeff + else + z = (height_u_levels(k)-transition_lower_bound)/transition_width + epssm_coeff_u(k) = min_coeff + sin(0.5_RKIND*pii*z)**2*(max_coeff-min_coeff) + end if + etp(k) = 0.5*(1.0 + epssm_coeff_u(k)) + etm(k) = 0.5*(1.0 - epssm_coeff_u(k)) + call mpas_log_write('k, etp, etm $i $r $r ',intArgs=(/k/),realArgs=(/etp(k),etm(k)/)) + end do + do k= 1,nVertlevels+1 + if(height_w_levels(k).le.transition_lower_bound) then + epssm_coeff_w(k) = min_coeff + else if(height_w_levels(k).ge.transition_upper_bound) then + epssm_coeff_w(k) = max_coeff + else + z = (height_w_levels(k)-transition_lower_bound)/transition_width + epssm_coeff_w(k) = min_coeff + sin(0.5_RKIND*pii*z)**2*(max_coeff-min_coeff) + end if + ewp(k) = 0.5*(1.0 + epssm_coeff_w(k)) + ewm(k) = 0.5*(1.0 - epssm_coeff_w(k)) + call mpas_log_write('k, ewp, ewm $i $r $r ',intArgs=(/k/),realArgs=(/ewp(k),ewm(k)/)) + end do + + deallocate(height_u_levels) + deallocate(epssm_coeff_u) + deallocate(height_w_levels) + deallocate(epssm_coeff_w) + end subroutine atm_compute_damping_coefs From 097f9fedf391eaec3cf429172b2dc152019a054b Mon Sep 17 00:00:00 2001 From: Bill Skamarock Date: Tue, 14 Oct 2025 11:06:28 -0600 Subject: [PATCH 3/4] Variable epssm(zeta) introduction, step (3): modified code to set the coefficients for the vertically-implicit acoustic solution to use the variable epssm formulation, and modified the acoustic step solver to use the variable epssm coefficients. This commit changes the solution at roundoff level even if a constant epssm is configured because we have changed the order of the computations in some places and we have incorporated the acoustic timestep delta tau in a different manner than in the previous formulation. --- .../dynamics/mpas_atm_time_integration.F | 157 +++++++++++++----- 1 file changed, 120 insertions(+), 37 deletions(-) diff --git a/src/core_atmosphere/dynamics/mpas_atm_time_integration.F b/src/core_atmosphere/dynamics/mpas_atm_time_integration.F index 4fe2faefc4..353c3574cc 100644 --- a/src/core_atmosphere/dynamics/mpas_atm_time_integration.F +++ b/src/core_atmosphere/dynamics/mpas_atm_time_integration.F @@ -2177,10 +2177,16 @@ subroutine atm_compute_vert_imp_coefs(state, mesh, diag, configs, nVertLevels, d real (kind=RKIND), pointer :: epssm - integer, pointer :: nCells, moist_start, moist_end + ! variable epssm arrays + real (kind=RKIND), dimension(:), pointer :: etp, etm, ewp, ewm + integer, pointer :: nCells, moist_start, moist_end call mpas_pool_get_config(configs, 'config_epssm', epssm) + call mpas_pool_get_array(mesh, 'etp', etp) + call mpas_pool_get_array(mesh, 'etm', etm) + call mpas_pool_get_array(mesh, 'ewp', ewp) + call mpas_pool_get_array(mesh, 'ewm', ewm) call mpas_pool_get_array(mesh, 'rdzu', rdzu) call mpas_pool_get_array(mesh, 'rdzw', rdzw) @@ -2215,6 +2221,7 @@ subroutine atm_compute_vert_imp_coefs(state, mesh, diag, configs, nVertLevels, d call atm_compute_vert_imp_coefs_work(nCells, moist_start, moist_end, dts, epssm, & zz, cqw, p, t, rb, rtb, pb, rt, cofwr, cofwz, coftz, cofwt, & a_tri, alpha_tri, gamma_tri, cofrz, rdzw, fzm, fzp, rdzu, scalars, & + etp, etm, ewp, ewm, & cellStart, cellEnd, edgeStart, edgeEnd, & cellSolveStart, cellSolveEnd, edgeSolveStart, edgeSolveEnd) @@ -2225,6 +2232,7 @@ end subroutine atm_compute_vert_imp_coefs subroutine atm_compute_vert_imp_coefs_work(nCells, moist_start, moist_end, dts, epssm, & zz, cqw, p, t, rb, rtb, pb, rt, cofwr, cofwz, coftz, cofwt, & a_tri, alpha_tri, gamma_tri, cofrz, rdzw, fzm, fzp, rdzu, scalars, & + etp, etm, ewp, ewm, & cellStart, cellEnd, edgeStart, edgeEnd, & cellSolveStart, cellSolveEnd, edgeSolveStart, edgeSolveEnd) @@ -2260,6 +2268,8 @@ subroutine atm_compute_vert_imp_coefs_work(nCells, moist_start, moist_end, dts, real (kind=RKIND), dimension(nVertLevels) :: fzm real (kind=RKIND), dimension(nVertLevels) :: fzp real (kind=RKIND), dimension(nVertLevels) :: rdzu + real (kind=RKIND), dimension(nVertLevels ) :: etp,etm + real (kind=RKIND), dimension(nVertLevels+1) :: ewp,ewm real (kind=RKIND), dimension(num_scalars,nVertLevels,nCells+1) :: scalars integer, intent(in) :: cellStart, cellEnd, edgeStart, edgeEnd @@ -2280,7 +2290,7 @@ subroutine atm_compute_vert_imp_coefs_work(nCells, moist_start, moist_end, dts, MPAS_ACC_TIMER_STOP('atm_compute_vert_imp_coefs_work [ACC_data_xfer]') ! set coefficients - dtseps = .5*dts*(1.+epssm) + ! dtseps = .5*dts*(1.+epssm) ! not needed for epssm_z rcv = rgas/(cp-rgas) c2 = cp*rcv @@ -2288,7 +2298,8 @@ subroutine atm_compute_vert_imp_coefs_work(nCells, moist_start, moist_end, dts, !$acc loop gang worker ! MGD bad to have all threads setting this variable? do k=1,nVertLevels - cofrz(k) = dtseps*rdzw(k) + ! cofrz(k) = dtseps*rdzw(k) + cofrz(k) = rdzw(k) ! epssm_z_change end do !$acc end parallel @@ -2299,15 +2310,21 @@ subroutine atm_compute_vert_imp_coefs_work(nCells, moist_start, moist_end, dts, !DIR$ IVDEP !$acc loop vector do k=2,nVertLevels - cofwr(k,iCell) =.5*dtseps*gravity*(fzm(k)*zz(k,iCell)+fzp(k)*zz(k-1,iCell)) + ! cofwr(k,iCell) =.5*dtseps*gravity*(fzm(k)*zz(k,iCell)+fzp(k)*zz(k-1,iCell)) + cofwr(k,iCell) =.5*gravity*(fzm(k)*zz(k,iCell)+fzp(k)*zz(k-1,iCell)) ! epssm_z_change end do coftz(1,iCell) = 0.0 !DIR$ IVDEP !$acc loop vector do k=2,nVertLevels - cofwz(k,iCell) = dtseps*c2*(fzm(k)*zz(k,iCell)+fzp(k)*zz(k-1,iCell)) & + ! cofwz(k,iCell) = dtseps*c2*(fzm(k)*zz(k,iCell)+fzp(k)*zz(k-1,iCell)) & + ! *rdzu(k)*cqw(k,iCell)*(fzm(k)*p (k,iCell)+fzp(k)*p (k-1,iCell)) + ! coftz(k,iCell) = dtseps* (fzm(k)*t (k,iCell)+fzp(k)*t (k-1,iCell)) + + ! ! epssm_z_change + cofwz(k,iCell) = c2*(fzm(k)*zz(k,iCell)+fzp(k)*zz(k-1,iCell)) & *rdzu(k)*cqw(k,iCell)*(fzm(k)*p (k,iCell)+fzp(k)*p (k-1,iCell)) - coftz(k,iCell) = dtseps* (fzm(k)*t (k,iCell)+fzp(k)*t (k-1,iCell)) + coftz(k,iCell) = (fzm(k)*t (k,iCell)+fzp(k)*t (k-1,iCell)) end do coftz(nVertLevels+1,iCell) = 0.0 !DIR$ IVDEP @@ -2320,9 +2337,14 @@ subroutine atm_compute_vert_imp_coefs_work(nCells, moist_start, moist_end, dts, ! end do qtotal = qtot(k,iCell) - cofwt(k,iCell) = .5*dtseps*rcv*zz(k,iCell)*gravity*rb(k,iCell)/(1.+qtotal) & + ! cofwt(k,iCell) = .5*dtseps*rcv*zz(k,iCell)*gravity*rb(k,iCell)/(1.+qtotal) & + ! *p(k,iCell)/((rtb(k,iCell)+rt(k,iCell))*pb(k,iCell)) + + ! epssm_z_change + cofwt(k,iCell) = .5*rcv*zz(k,iCell)*gravity*rb(k,iCell)/(1.+qtotal) & *p(k,iCell)/((rtb(k,iCell)+rt(k,iCell))*pb(k,iCell)) -! cofwt(k,iCell) = 0. + ! cofwt(k,iCell) = .5*rcv*zz(k,iCell)*gravity/t(k,iCell) ! zero base state option + ! cofwt(k,iCell) = 0. end do a_tri(1,iCell) = 0. ! note, this value is never used @@ -2337,21 +2359,38 @@ subroutine atm_compute_vert_imp_coefs_work(nCells, moist_start, moist_end, dts, a_tri(k,iCell) = -cofwz(k ,iCell)* coftz(k-1,iCell)*rdzw(k-1)*zz(k-1,iCell) & +cofwr(k ,iCell)* cofrz(k-1 ) & -cofwt(k-1,iCell)* coftz(k-1,iCell)*rdzw(k-1) - b_tri(k) = 1. & - +cofwz(k ,iCell)*(coftz(k ,iCell)*rdzw(k )*zz(k ,iCell) & - +coftz(k ,iCell)*rdzw(k-1)*zz(k-1,iCell)) & - -coftz(k ,iCell)*(cofwt(k ,iCell)*rdzw(k ) & - -cofwt(k-1,iCell)*rdzw(k-1)) & - +cofwr(k, iCell)*(cofrz(k )-cofrz(k-1)) + a_tri(k,iCell) = a_tri(k,iCell)*etp(k-1)*ewp(k-1) ! epssm_z_change (addition) + + ! b_tri(k) = 1. & + ! +cofwz(k ,iCell)*(coftz(k ,iCell)*rdzw(k )*zz(k ,iCell) & + ! +coftz(k ,iCell)*rdzw(k-1)*zz(k-1,iCell)) & + ! -coftz(k ,iCell)*(cofwt(k ,iCell)*rdzw(k ) & + ! -cofwt(k-1,iCell)*rdzw(k-1)) & + ! +cofwr(k, iCell)*(cofrz(k )-cofrz(k-1)) + + ! epssm_z_change + b_tri(k) = +cofwz(k ,iCell)*coftz(k,iCell)* & + ( etp(k )*rdzw(k )*zz(k ,iCell) & + +etp(k-1)*rdzw(k-1)*zz(k-1,iCell)) & + -coftz(k ,iCell)*( etp(k )*cofwt(k ,iCell)*rdzw(k ) & + -etp(k-1)*cofwt(k-1,iCell)*rdzw(k-1)) & + +cofwr(k, iCell)*(etp(k)*cofrz(k )-etp(k-1)*cofrz(k-1)) + b_tri(k) = b_tri(k)*ewp(k) + c_tri(k) = -cofwz(k ,iCell)* coftz(k+1,iCell)*rdzw(k )*zz(k ,iCell) & -cofwr(k ,iCell)* cofrz(k ) & +cofwt(k ,iCell)* coftz(k+1,iCell)*rdzw(k ) + c_tri(k) = c_tri(k)*etp(k)*ewp(k+1) ! epssm_z_change (addition) end do + c_tri(nVertLevels) = 0.0 ! epssm_z_change (addition) !MGD VECTOR DEPENDENCE !$acc loop seq do k=2,nVertLevels - alpha_tri(k,iCell) = 1./(b_tri(k)-a_tri(k,iCell)*gamma_tri(k-1,iCell)) - gamma_tri(k,iCell) = c_tri(k)*alpha_tri(k,iCell) + ! alpha_tri(k,iCell) = 1./(b_tri(k)-a_tri(k,iCell)*gamma_tri(k-1,iCell)) + ! gamma_tri(k,iCell) = c_tri(k)*alpha_tri(k,iCell) + ! epssm_z_change + alpha_tri(k,iCell) = 1./(1.0+(dts**2)*(b_tri(k)-a_tri(k,iCell)*gamma_tri(k-1,iCell))) + gamma_tri(k,iCell) = (dts**2)*c_tri(k)*alpha_tri(k,iCell) end do end do ! loop over cells @@ -2561,6 +2600,8 @@ subroutine atm_advance_acoustic_step( state, diag, tend, mesh, configs, nCells, real (kind=RKIND), pointer :: cf1, cf2, cf3 + real (kind=RKIND), dimension(:), pointer :: etp, etm, ewp, ewm + integer, pointer :: nEdges, nCellsSolve call mpas_pool_get_array(mesh, 'cellsOnEdge', cellsOnEdge) @@ -2570,6 +2611,11 @@ subroutine atm_advance_acoustic_step( state, diag, tend, mesh, configs, nCells, call mpas_pool_get_array(mesh, 'specZoneMaskEdge', specZoneMaskEdge) call mpas_pool_get_array(mesh, 'specZoneMaskCell', specZoneMaskCell) + call mpas_pool_get_array(mesh, 'etp', etp) + call mpas_pool_get_array(mesh, 'etm', etm) + call mpas_pool_get_array(mesh, 'ewp', ewp) + call mpas_pool_get_array(mesh, 'ewm', ewm) + call mpas_pool_get_array(state, 'rho_zz', rho_zz, 2) ! call mpas_pool_get_array(state, 'theta_m', theta_m, 2) call mpas_pool_get_array(state, 'theta_m', theta_m, 1) @@ -2637,6 +2683,7 @@ subroutine atm_advance_acoustic_step( state, diag, tend, mesh, configs, nCells, tend_rw, zgrid, cofwr, cofwz, w, ru, ru_save, rw, rw_save, fzm, fzp, rdzw, dcEdge, invDcEdge, & invAreaCell, cofrz, dvEdge, nEdgesOnCell, cellsOnEdge, edgesOnCell, edgesOnCell_sign, & dts, small_step, epssm, cf1, cf2, cf3, & + etp, etm, ewp, ewm, & specZoneMaskEdge, specZoneMaskCell & ) @@ -2650,6 +2697,7 @@ subroutine atm_advance_acoustic_step_work(nCells, nEdges, nCellsSolve, cellStart tend_rw, zgrid, cofwr, cofwz, w, ru, ru_save, rw, rw_save, fzm, fzp, rdzw, dcEdge, invDcEdge, & invAreaCell, cofrz, dvEdge, nEdgesOnCell, cellsOnEdge, edgesOnCell, edgesOnCell_sign, & dts, small_step, epssm, cf1, cf2, cf3, & + etp, etm, ewp, ewm, & specZoneMaskEdge, specZoneMaskCell & ) @@ -2703,6 +2751,10 @@ subroutine atm_advance_acoustic_step_work(nCells, nEdges, nCellsSolve, cellStart real (kind=RKIND), dimension(nVertLevels) :: fzm real (kind=RKIND), dimension(nVertLevels) :: fzp real (kind=RKIND), dimension(nVertLevels) :: rdzw + real (kind=RKIND), dimension(nVertLevels ) :: etp + real (kind=RKIND), dimension(nVertLevels ) :: etm + real (kind=RKIND), dimension(nVertLevels+1) :: ewp + real (kind=RKIND), dimension(nVertLevels+1) :: ewm real (kind=RKIND), dimension(nEdges+1) :: dcEdge real (kind=RKIND), dimension(nEdges+1) :: invDcEdge real (kind=RKIND), dimension(nCells+1) :: invAreaCell @@ -2888,31 +2940,51 @@ subroutine atm_advance_acoustic_step_work(nCells, nEdges, nCellsSolve, cellStart !DIR$ IVDEP !$acc loop vector do k=1, nVertLevels - rs(k) = rho_pp(k,iCell) + dts*tend_rho(k,iCell) + rs(k) & - - cofrz(k)*resm*(rw_p(k+1,iCell)-rw_p(k,iCell)) - ts(k) = rtheta_pp(k,iCell) + dts*tend_rt(k,iCell) + ts(k) & - - resm*rdzw(k)*( coftz(k+1,iCell)*rw_p(k+1,iCell) & - -coftz(k,iCell)*rw_p(k,iCell)) + ! rs(k) = rho_pp(k,iCell) + dts*tend_rho(k,iCell) + rs(k) & + ! - cofrz(k)*resm*(rw_p(k+1,iCell)-rw_p(k,iCell)) + ! ts(k) = rtheta_pp(k,iCell) + dts*tend_rt(k,iCell) + ts(k) & + ! - resm*rdzw(k)*( coftz(k+1,iCell)*rw_p(k+1,iCell) & + ! -coftz(k,iCell)*rw_p(k,iCell)) + + ! epssm_z change + rs(k) = rho_pp(k,iCell) + dts*tend_rho(k,iCell) + rs(k) & + - dts*cofrz(k)*(ewm(k+1)*rw_p(k+1,iCell)-ewm(k)*rw_p(k,iCell)) + ts(k) = rtheta_pp(k,iCell) + dts*tend_rt(k,iCell) + ts(k) & + - dts*rdzw(k)*( ewm(k+1)*coftz(k+1,iCell)*rw_p(k+1,iCell) & + -ewm(k )*coftz(k,iCell)*rw_p(k,iCell)) end do !DIR$ IVDEP !$acc loop vector do k=2, nVertLevels - wwAvg(k,iCell) = wwAvg(k,iCell) + 0.5*(1.0-epssm)*rw_p(k,iCell) + ! wwAvg(k,iCell) = wwAvg(k,iCell) + 0.5*(1.0-epssm)*rw_p(k,iCell) + wwavg(k,iCell) = wwavg(k,iCell) + ewm(k)*rw_p(k,iCell) ! epssm_z change end do !DIR$ IVDEP !$acc loop vector do k=2, nVertLevels - rw_p(k,iCell) = rw_p(k,iCell) + dts*tend_rw(k,iCell) & - - cofwz(k,iCell)*((zz(k ,iCell)*ts(k) & - -zz(k-1,iCell)*ts(k-1)) & - +resm*(zz(k ,iCell)*rtheta_pp(k ,iCell) & - -zz(k-1,iCell)*rtheta_pp(k-1,iCell))) & - - cofwr(k,iCell)*((rs(k)+rs(k-1)) & - +resm*(rho_pp(k,iCell)+rho_pp(k-1,iCell))) & - + cofwt(k ,iCell)*(ts(k )+resm*rtheta_pp(k ,iCell)) & - + cofwt(k-1,iCell)*(ts(k-1)+resm*rtheta_pp(k-1,iCell)) + ! rw_p(k,iCell) = rw_p(k,iCell) + dts*tend_rw(k,iCell) & + ! - cofwz(k,iCell)*((zz(k ,iCell)*ts(k) & + ! -zz(k-1,iCell)*ts(k-1)) & + ! +resm*(zz(k ,iCell)*rtheta_pp(k ,iCell) & + ! -zz(k-1,iCell)*rtheta_pp(k-1,iCell))) & + ! - cofwr(k,iCell)*((rs(k)+rs(k-1)) & + ! +resm*(rho_pp(k,iCell)+rho_pp(k-1,iCell))) & + ! + cofwt(k ,iCell)*(ts(k )+resm*rtheta_pp(k ,iCell)) & + ! + cofwt(k-1,iCell)*(ts(k-1)+resm*rtheta_pp(k-1,iCell)) + + ! epssm_z change + rw_p(k,iCell) = rw_p(k,iCell) + dts*(tend_rw(k,iCell) & + - cofwz(k,iCell)*(( etp(k )*zz(k ,iCell)*ts(k) & + -etp(k-1)*zz(k-1,iCell)*ts(k-1)) & + + ( etm(k )*zz(k ,iCell)*rtheta_pp(k ,iCell) & + -etm(k-1)*zz(k-1,iCell)*rtheta_pp(k-1,iCell))) & + - cofwr(k,iCell)*((etp(k)*rs(k)+etp(k-1)*rs(k-1)) & + +( etm(k )*rho_pp(k ,iCell) & + +etm(k-1)*rho_pp(k-1,iCell))) & + + cofwt(k ,iCell)*(etp(k )*ts(k )+etm(k )*rtheta_pp(k ,iCell)) & + + cofwt(k-1,iCell)*(etp(k-1)*ts(k-1)+etm(k-1)*rtheta_pp(k-1,iCell))) end do ! tridiagonal solve sweeping up and then down the column @@ -2920,7 +2992,10 @@ subroutine atm_advance_acoustic_step_work(nCells, nEdges, nCellsSolve, cellStart !MGD VECTOR DEPENDENCE !$acc loop seq do k=2,nVertLevels - rw_p(k,iCell) = (rw_p(k,iCell)-a_tri(k,iCell)*rw_p(k-1,iCell))*alpha_tri(k,iCell) + ! rw_p(k,iCell) = (rw_p(k,iCell)-a_tri(k,iCell)*rw_p(k-1,iCell))*alpha_tri(k,iCell) + + ! epssm_z change + rw_p(k,iCell) = (rw_p(k,iCell)-(dts**2)*a_tri(k,iCell)*rw_p(k-1,iCell))*alpha_tri(k,iCell) end do !MGD VECTOR DEPENDENCE @@ -2945,7 +3020,8 @@ subroutine atm_advance_acoustic_step_work(nCells, nEdges, nCellsSolve, cellStart !DIR$ IVDEP !$acc loop vector do k=2,nVertLevels - wwAvg(k,iCell) = wwAvg(k,iCell) + 0.5*(1.0+epssm)*rw_p(k,iCell) + ! wwAvg(k,iCell) = wwAvg(k,iCell) + 0.5*(1.0+epssm)*rw_p(k,iCell) + wwAvg(k,iCell) = wwAvg(k,iCell) + ewp(k)*rw_p(k,iCell) ! epssm_z change end do ! update rho_pp and theta_pp given updated rw_p @@ -2953,9 +3029,15 @@ subroutine atm_advance_acoustic_step_work(nCells, nEdges, nCellsSolve, cellStart !DIR$ IVDEP !$acc loop vector do k=1,nVertLevels - rho_pp(k,iCell) = rs(k) - cofrz(k) *(rw_p(k+1,iCell)-rw_p(k ,iCell)) - rtheta_pp(k,iCell) = ts(k) - rdzw(k)*(coftz(k+1,iCell)*rw_p(k+1,iCell) & - -coftz(k ,iCell)*rw_p(k ,iCell)) + ! rho_pp(k,iCell) = rs(k) - cofrz(k) *(rw_p(k+1,iCell)-rw_p(k ,iCell)) + ! rtheta_pp(k,iCell) = ts(k) - rdzw(k)*(coftz(k+1,iCell)*rw_p(k+1,iCell) & + ! -coftz(k ,iCell)*rw_p(k ,iCell)) + ! + ! epssm_z change + rho_pp(k,iCell) = rs(k) - dts*cofrz(k) *( ewp(k+1)*rw_p(k+1,iCell) & + -ewp(k )*rw_p(k ,iCell)) + rtheta_pp(k,iCell) = ts(k) - dts*rdzw(k)*( ewp(k+1)*coftz(k+1,iCell)*rw_p(k+1,iCell) & + -ewp(k )*coftz(k ,iCell)*rw_p(k ,iCell)) end do else ! specifed zone in regional_MPAS @@ -2965,7 +3047,8 @@ subroutine atm_advance_acoustic_step_work(nCells, nEdges, nCellsSolve, cellStart rho_pp(k,iCell) = rho_pp(k,iCell) + dts*tend_rho(k,iCell) rtheta_pp(k,iCell) = rtheta_pp(k,iCell) + dts*tend_rt(k,iCell) rw_p(k,iCell) = rw_p(k,iCell) + dts*tend_rw(k,iCell) - wwAvg(k,iCell) = wwAvg(k,iCell) + 0.5*(1.0+epssm)*rw_p(k,iCell) + ! wwAvg(k,iCell) = wwAvg(k,iCell) + 0.5*(1.0+epssm)*rw_p(k,iCell) + wwAvg(k,iCell) = wwAvg(k,iCell) + ewp(k)*rw_p(k,iCell) ! epssm_z change end do end if From 51d8389060519b1e97472d34239515b3ab549963 Mon Sep 17 00:00:00 2001 From: Bill Skamarock Date: Tue, 14 Oct 2025 11:31:53 -0600 Subject: [PATCH 4/4] Changed the default values of the variable_epssm configuration variables to the intended general configuration for real-data cases. --- src/core_atmosphere/Registry.xml | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/core_atmosphere/Registry.xml b/src/core_atmosphere/Registry.xml index fc599be2be..f5d74ef379 100644 --- a/src/core_atmosphere/Registry.xml +++ b/src/core_atmosphere/Registry.xml @@ -298,17 +298,17 @@ description="Value of epssm below transition zone" possible_values="Positive real values between 0 and 1"/> - - -