diff --git a/src/core_atmosphere/Registry.xml b/src/core_atmosphere/Registry.xml
index 4281c40bba..f5d74ef379 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"/>
+
+                <nml_option name="config_epssm_minimum" type="real" default_value="0.1"
+                     units="-"
+                     description="Value of epssm below transition zone"
+                     possible_values="Positive real values between 0 and 1"/>
+
+                <nml_option name="config_epssm_maximum" type="real" default_value="0.5"
+                     units="-"
+                     description="Value of epssm above transition zone"
+                     possible_values="Positive real values between 0 and 1"/>
+
+                <nml_option name="config_epssm_transition_bottom_z" type="real" default_value="30000."
+                     units="m"
+                     description="Height MSL of bottom of transition zone for epssm"
+                     possible_values="Positive real values"/>
+
+                <nml_option name="config_epssm_transition_top_z" type="real" default_value="50000."
+                     units="m"
+                     description="Height MSL of top of transition zone for epssm"
+                     possible_values="Positive real values"/>
+
         </nml_record>
 
         <nml_record name="limited_area" in_defaults="true">
@@ -495,6 +516,10 @@
 			<var name="defc_a"/>
 			<var name="defc_b"/>
 			<var name="coeffs_reconstruct"/>
+			<var name="etp"/>
+			<var name="etm"/>
+			<var name="ewp"/>
+			<var name="ewm"/>
 #ifdef MPAS_CAM_DYCORE
 			<var name="cell_gradient_coef_x"/>
 			<var name="cell_gradient_coef_y"/>
@@ -1484,6 +1509,18 @@
                 <var name="rdzw" type="real" dimensions="nVertLevels" units="unitless"
                      description="Reciprocal dzw"/>
 
+                <var name="etp" type="real" dimensions="nVertLevels" units="unitless"
+                     description="etp = (1+epssm(z))/2 at theta point"/>
+
+                <var name="etm" type="real" dimensions="nVertLevels" units="unitless"
+                     description="etm = (1-epssm(z))/2 at theta point"/>
+
+                <var name="ewp" type="real" dimensions="nVertLevelsP1" units="unitless"
+                     description="ewp = (1+epssm(z))/2 at w point"/>
+
+                <var name="ewm" type="real" dimensions="nVertLevelsP1" units="unitless"
+                     description="ewm = (1-epssm(z))/2 at w point"/>
+
                 <var name="dzu" type="real" dimensions="nVertLevels" units="unitless"
                      description="d(zeta) at w levels"/>
 
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
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