diff --git a/src/core_atmosphere/Registry.xml b/src/core_atmosphere/Registry.xml
index 4e50bea204..a8f08f9c96 100644
--- a/src/core_atmosphere/Registry.xml
+++ b/src/core_atmosphere/Registry.xml
@@ -251,7 +251,8 @@
units="-"
description="Coefficient for the divergent component of the Laplacian filter of momentum in the relaxation zone"
possible_values="Positive real values"/>
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
@@ -487,6 +512,7 @@
+
#ifdef MPAS_CAM_DYCORE
@@ -607,6 +633,10 @@
+
+
+
+
@@ -623,6 +653,7 @@
+
#ifdef MPAS_CAM_DYCORE
@@ -1353,6 +1384,18 @@
+
+
+
+
+
+
+
+
@@ -1400,6 +1443,9 @@
+
+
diff --git a/src/core_atmosphere/dynamics/mpas_atm_time_integration.F b/src/core_atmosphere/dynamics/mpas_atm_time_integration.F
index 63b2ddd25e..337489f9c1 100644
--- a/src/core_atmosphere/dynamics/mpas_atm_time_integration.F
+++ b/src/core_atmosphere/dynamics/mpas_atm_time_integration.F
@@ -94,7 +94,6 @@ end subroutine halo_exchange_routine
! if damping in the model changed with the simulation calendar
real (kind=RKIND), parameter, private :: seconds_per_day = 86400.0_RKIND
-
contains
@@ -376,6 +375,11 @@ subroutine atm_srk3(domain, dt, itimestep, exchange_halo_group)
real (kind=RKIND), dimension(:,:), pointer :: rqvdynten, rthdynten, theta_m
real (kind=RKIND) :: theta_local, fac_m
+!--------------------------
+!SK_epssm: wrote this to retrieve rdzw
+ real (kind=RKIND), dimension(:), pointer:: rdzw
+!--------------------------
+
#ifndef MPAS_CAM_DYCORE
! Used in allocating scratch fields for physics tendencies
type (field2DReal), pointer :: tend_ru_physicsField, tend_rtheta_physicsField, tend_rho_physicsField
@@ -415,7 +419,6 @@ subroutine atm_srk3(domain, dt, itimestep, exchange_halo_group)
#ifdef DO_PHYSICS
call mpas_pool_get_subpool(block % structs, 'diag_physics', diag_physics)
#endif
-
!
! Retrieve dimensions
! Note: nCellsSolve and nVerticesSolve are not currently used in this function
@@ -1603,9 +1606,16 @@ subroutine atm_compute_vert_imp_coefs(state, mesh, diag, configs, nVertLevels, d
integer, pointer :: nCells, moist_start, moist_end
-
+!!!JBK-vepssm and reference sounding
+ real (kind=RKIND), dimension(:), pointer :: etp, etm, ewp, ewm
+ integer, pointer :: reference_sounding_switch
+ 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, 'reference_sounding_switch', reference_sounding_switch)
call mpas_pool_get_config(configs, 'config_epssm', epssm)
-
+!!!JBK
call mpas_pool_get_array(mesh, 'rdzu', rdzu)
call mpas_pool_get_array(mesh, 'rdzw', rdzw)
call mpas_pool_get_array(mesh, 'fzm', fzm)
@@ -1635,11 +1645,13 @@ subroutine atm_compute_vert_imp_coefs(state, mesh, diag, configs, nVertLevels, d
call mpas_pool_get_dimension(state, 'moist_start', moist_start)
call mpas_pool_get_dimension(state, 'moist_end', moist_end)
-
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, &
cellStart, cellEnd, edgeStart, edgeEnd, &
+!!!JBK
+ etp, etm, ewp, ewm, reference_sounding_switch, &
+!!!JBK
cellSolveStart, cellSolveEnd, edgeSolveStart, edgeSolveEnd)
@@ -1650,6 +1662,9 @@ subroutine atm_compute_vert_imp_coefs_work(nCells, moist_start, moist_end, dts,
zz, cqw, p, t, rb, rtb, pb, rt, cofwr, cofwz, coftz, cofwt, &
a_tri, alpha_tri, gamma_tri, cofrz, rdzw, fzm, fzp, rdzu, scalars, &
cellStart, cellEnd, edgeStart, edgeEnd, &
+!!!JBK
+ etp, etm, ewp, ewm, reference_sounding_switch, &
+!!!JBK
cellSolveStart, cellSolveEnd, edgeSolveStart, edgeSolveEnd)
use mpas_atm_dimensions
@@ -1688,6 +1703,12 @@ subroutine atm_compute_vert_imp_coefs_work(nCells, moist_start, moist_end, dts,
integer, intent(in) :: cellStart, cellEnd, edgeStart, edgeEnd
integer, intent(in) :: cellSolveStart, cellSolveEnd, edgeSolveStart, edgeSolveEnd
+!!!JBK
+ real (kind=RKIND), dimension(nVertLevels ) :: etp,etm
+ real (kind=RKIND), dimension(nVertLevels+1) :: ewp,ewm
+!!!JBK Swith to use reference sounding and perturbation thermodynamic variables
+ integer :: reference_sounding_switch
+!!!JBK
!
@@ -1697,7 +1718,6 @@ subroutine atm_compute_vert_imp_coefs_work(nCells, moist_start, moist_end, dts,
real (kind=RKIND) :: dtseps, c2, qtotal, rcv
real (kind=RKIND), dimension( nVertLevels ) :: b_tri, c_tri
-
! set coefficients
dtseps = .5*dts*(1.+epssm)
rcv = rgas/(cp-rgas)
@@ -1705,24 +1725,30 @@ subroutine atm_compute_vert_imp_coefs_work(nCells, moist_start, moist_end, dts,
! MGD bad to have all threads setting this variable?
do k=1,nVertLevels
- cofrz(k) = dtseps*rdzw(k)
+ ! WCS_epssm
+ ! cofrz(k) = dtseps*rdzw(k)
+ cofrz(k) = rdzw(k)
end do
do iCell = cellSolveStart,cellSolveEnd ! we only need to do cells we are solving for, not halo cells
!DIR$ IVDEP
do k=2,nVertLevels
- cofwr(k,iCell) =.5*dtseps*gravity*(fzm(k)*zz(k,iCell)+fzp(k)*zz(k-1,iCell))
+ cofwr(k,iCell) = 0.5*gravity*(fzm(k)*zz(k,iCell)+fzp(k)*zz(k-1,iCell))
end do
coftz(1,iCell) = 0.0
!DIR$ IVDEP
do k=2,nVertLevels
- 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))
+ 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) = (fzm(k)*t (k,iCell)+fzp(k)*t (k-1,iCell))
end do
coftz(nVertLevels+1,iCell) = 0.0
!DIR$ IVDEP
+
+!!!JBK Mod to use either reference sounding or full variables and log(p)
+ if(reference_sounding_switch.eq.1) then
+!!!JBK
do k=1,nVertLevels
! qtotal = 0.
@@ -1731,10 +1757,19 @@ 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) &
- *p(k,iCell)/((rtb(k,iCell)+rt(k,iCell))*pb(k,iCell))
+ 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.
end do
+!!!JBK Mod to use full variables and log(p)
+ else
+ do k=1,nVertLevels
+ cofwt(k,iCell) = .5*rcv*zz(k,iCell)*gravity/t(k,iCell)
+ cofwt(k,iCell) = 0.
+ end do
+ end if
+!!!JBK end
a_tri(1,iCell) = 0. ! note, this value is never used
b_tri(1) = 1. ! note, this value is never used
@@ -1744,23 +1779,29 @@ subroutine atm_compute_vert_imp_coefs_work(nCells, moist_start, moist_end, dts,
!DIR$ IVDEP
do k=2,nVertLevels
- a_tri(k,iCell) = -cofwz(k ,iCell)* coftz(k-1,iCell)*rdzw(k-1)*zz(k-1,iCell) &
- +cofwr(k ,iCell)* cofrz(k-1 ) &
+ 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))
- c_tri(k) = -cofwz(k ,iCell)* coftz(k+1,iCell)*rdzw(k )*zz(k ,iCell) &
- -cofwr(k ,iCell)* cofrz(k ) &
+ a_tri(k,iCell) = a_tri(k,iCell)*etp(k-1)*ewp(k-1)
+
+ 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)
end do
+ c_tri(nVertLevels) = 0.0
!MGD VECTOR DEPENDENCE
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./(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
@@ -2001,6 +2042,14 @@ subroutine atm_advance_acoustic_step( state, diag, tend, mesh, configs, nCells,
integer, pointer :: nEdges, nCellsSolve
+!!!JBK-vepssm
+ real (kind=RKIND), dimension(:), pointer :: etp, etm, ewp, ewm
+ 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)
+!!!JBK
+
call mpas_pool_get_array(mesh, 'cellsOnEdge', cellsOnEdge)
call mpas_pool_get_array(mesh, 'edgesOnCell', edgesOnCell)
call mpas_pool_get_array(mesh, 'nEdgesOnCell', nEdgesOnCell)
@@ -2075,6 +2124,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 &
)
@@ -2088,6 +2138,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 &
)
@@ -2159,8 +2210,10 @@ subroutine atm_advance_acoustic_step_work(nCells, nEdges, nCellsSolve, cellStart
integer, intent(in) :: small_step
real (kind=RKIND), intent(in) :: dts, epssm,cf1, cf2, cf3
real (kind=RKIND), dimension(nVertLevels) :: ts, rs
-
-
+!!!JBK
+ real (kind=RKIND), dimension(nVertLevels ) :: etp,etm
+ real (kind=RKIND), dimension(nVertLevels+1) :: ewp,ewm
+!!!JBK
!
! Local variables
!
@@ -2281,36 +2334,39 @@ subroutine atm_advance_acoustic_step_work(nCells, nEdges, nCellsSolve, cellStart
!DIR$ IVDEP
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) &
+ - 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
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) + ewm(k)*rw_p(k,iCell)
end do
!DIR$ IVDEP
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)*(( 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
!MGD VECTOR DEPENDENCE
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)-(dts**2)*a_tri(k,iCell)*rw_p(k-1,iCell))*alpha_tri(k,iCell)
+
end do
!MGD VECTOR DEPENDENCE
@@ -2332,16 +2388,17 @@ subroutine atm_advance_acoustic_step_work(nCells, nEdges, nCellsSolve, cellStart
! accumulate (rho*omega)' for use later in scalar transport
!DIR$ IVDEP
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) + ewp(k)*rw_p(k,iCell)
end do
! update rho_pp and theta_pp given updated rw_p
!DIR$ IVDEP
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) - 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
@@ -2350,7 +2407,7 @@ 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) + ewp(k)*rw_p(k,iCell)
end do
end if
@@ -3925,10 +3982,12 @@ subroutine atm_compute_dyn_tend(tend, tend_physics, state, diag, mesh, configs,
real (kind=RKIND), pointer :: config_h_theta_eddy_visc2, config_v_theta_eddy_visc2
real (kind=RKIND), pointer :: config_mpas_cam_coef
- logical, pointer :: config_rayleigh_damp_u
+ logical, pointer :: config_rayleigh_damp_u
real (kind=RKIND), pointer :: config_rayleigh_damp_u_timescale_days
- integer, pointer :: config_number_rayleigh_damp_u_levels, config_number_cam_damping_levels
-
+ integer, pointer :: config_number_rayleigh_damp_u_levels, config_number_cam_damping_levels
+!!!JBK - Mod for reference sounding
+ integer, pointer :: reference_sounding_switch
+!!!JBK end
call mpas_pool_get_config(mesh, 'sphere_radius', r_earth)
call mpas_pool_get_config(configs, 'config_coef_3rd_order', coef_3rd_order)
@@ -3979,7 +4038,9 @@ subroutine atm_compute_dyn_tend(tend, tend_physics, state, diag, mesh, configs,
call mpas_pool_get_array(diag, 'exner', exner)
call mpas_pool_get_array(tend_physics, 'rthdynten', rthdynten)
-
+!!!JBK Mod for reference_sounding_switch
+ call mpas_pool_get_array(mesh, 'reference_sounding_switch', reference_sounding_switch)
+!!!
call mpas_pool_get_array(mesh, 'weightsOnEdge', weightsOnEdge)
call mpas_pool_get_array(mesh, 'cellsOnEdge', cellsOnEdge)
call mpas_pool_get_array(mesh, 'cellsOnCell', cellsOnCell)
@@ -4072,6 +4133,9 @@ subroutine atm_compute_dyn_tend(tend, tend_physics, state, diag, mesh, configs,
config_mpas_cam_coef, &
config_rayleigh_damp_u, config_rayleigh_damp_u_timescale_days, &
config_number_rayleigh_damp_u_levels, config_number_cam_damping_levels, &
+!!!JBK - Mod for reference sounding
+ reference_sounding_switch, &
+!!!JBK end
rthdynten, &
cellStart, cellEnd, vertexStart, vertexEnd, edgeStart, edgeEnd, &
cellSolveStart, cellSolveEnd, vertexSolveStart, vertexSolveEnd, edgeSolveStart, edgeSolveEnd)
@@ -4096,6 +4160,10 @@ subroutine atm_compute_dyn_tend_work(nCells, nEdges, nVertices, nVertLevels_dumm
config_mpas_cam_coef, &
config_rayleigh_damp_u, config_rayleigh_damp_u_timescale_days, &
config_number_rayleigh_damp_u_levels, config_number_cam_damping_levels, &
+
+!!!JBK - Mod for reference sounding
+ reference_sounding_switch, &
+!!!JBK end
rthdynten, &
cellStart, cellEnd, vertexStart, vertexEnd, edgeStart, edgeEnd, &
cellSolveStart, cellSolveEnd, vertexSolveStart, vertexSolveEnd, edgeSolveStart, edgeSolveEnd)
@@ -4253,6 +4321,12 @@ subroutine atm_compute_dyn_tend_work(nCells, nEdges, nVertices, nVertLevels_dumm
real (kind=RKIND) :: flux3, flux4
real (kind=RKIND) :: q_im2, q_im1, q_i, q_ip1, ua, coef3
+!!!JBK - mod for reference sounding switch
+ integer :: reference_sounding_switch
+! Declaration for plog
+ real (kind=RKIND) :: p0
+ real (kind=RKIND), dimension(nVertLevels,nCells+1) :: plog
+!!!JBK end
flux4(q_im2, q_im1, q_i, q_ip1, ua) = &
ua*( 7.*(q_i + q_im1) - (q_ip1 + q_im2) )/12.0
@@ -4269,6 +4343,12 @@ subroutine atm_compute_dyn_tend_work(nCells, nEdges, nVertices, nVertLevels_dumm
if (rk_step == 1) then
+ p0 = 1.0e+05 ! this should come from somewhere else...
+
+!!!JBK - Temporary mod to print reference_sounding
+ call mpas_log_write('reference_sounding_switch = $i ', intArgs=(/reference_sounding_switch/))
+!!!JBK end
+
! tend_u_euler(1:nVertLevels,edgeStart:edgeEnd) = 0.0
! Smagorinsky eddy viscosity, based on horizontal deformation (in this case on model coordinate surfaces).
@@ -4360,10 +4440,27 @@ subroutine atm_compute_dyn_tend_work(nCells, nEdges, nVertices, nVertLevels_dumm
do iCell = cellStart,cellEnd
!DIR$ IVDEP
+
+!!!JBK - Mods to use log(p)
+! do k = 1,nVertLevels
+! tend_rho(k,iCell) = -h_divergence(k,iCell)-rdzw(k)*(rw(k+1,iCell)-rw(k,iCell)) + tend_rho_physics(k,iCell)
+! dpdz(k,iCell) = -gravity*(rb(k,iCell)*(qtot(k,iCell)) + rr_save(k,iCell)*(1.+qtot(k,iCell)))
+! end do
+
do k = 1,nVertLevels
- tend_rho(k,iCell) = -h_divergence(k,iCell)-rdzw(k)*(rw(k+1,iCell)-rw(k,iCell)) + tend_rho_physics(k,iCell)
- dpdz(k,iCell) = -gravity*(rb(k,iCell)*(qtot(k,iCell)) + rr_save(k,iCell)*(1.+qtot(k,iCell)))
+ tend_rho(k,iCell) = -h_divergence(k,iCell)-rdzw(k)*(rw(k+1,iCell)-rw(k,iCell)) + tend_rho_physics(k,iCell)
end do
+ if(reference_sounding_switch.eq.1) then
+ do k = 1,nVertLevels
+ dpdz(k,iCell) = -gravity*(rb(k,iCell)*(qtot(k,iCell)) + rr_save(k,iCell)*(1.+qtot(k,iCell)))
+ end do
+ else
+ do k = 1,nVertLevels
+ dpdz(k,iCell) = -gravity*rr_save(k,iCell)/pp(k,iCell)
+ plog(k,iCell) = alog(pp(k,iCell)/p0)
+ end do
+ end if
+!!!JBK end
end do
end if
@@ -4382,11 +4479,27 @@ subroutine atm_compute_dyn_tend_work(nCells, nEdges, nVertices, nVertLevels_dumm
if(rk_step == 1) then
!DIR$ IVDEP
- do k=1,nVertLevels
- tend_u_euler(k,iEdge) = - cqu(k,iEdge)*( (pp(k,cell2)-pp(k,cell1))*invDcEdge(iEdge)/(.5*(zz(k,cell2)+zz(k,cell1))) &
- -0.5*zxu(k,iEdge)*(dpdz(k,cell1)+dpdz(k,cell2)) )
- end do
+!!!JBK - Mods to use log(p)
+! do k=1,nVertLevels
+! tend_u_euler(k,iEdge) = - cqu(k,iEdge)*( (pp(k,cell2)-pp(k,cell1))*invDcEdge(iEdge)/(.5*(zz(k,cell2)+zz(k,cell1))) &
+! -0.5*zxu(k,iEdge)*(dpdz(k,cell1)+dpdz(k,cell2)) )
+! end do
+ if(reference_sounding_switch.eq.1) then
+ do k=1,nVertLevels
+ tend_u_euler(k,iEdge) = - cqu(k,iEdge)*((pp(k,cell2)-pp(k,cell1))*invDcEdge(iEdge) &
+ /(.5*(zz(k,cell2)+zz(k,cell1))) &
+ -0.5*zxu(k,iEdge)*(dpdz(k,cell2)+dpdz(k,cell1)) )
+ end do
+ else
+ do k=1,nVertLevels
+ tend_u_euler(k,iEdge) = - .5*(pp (k,cell2)+pp (k,cell1))* &
+ (cqu(k,iEdge)*(plog(k,cell2)-plog(k,cell1))*invDcEdge(iEdge) &
+ /(.5*(zz (k,cell2)+zz (k,cell1))) &
+ -0.5*zxu(k,iEdge)*(dpdz(k,cell2)+dpdz(k,cell1)) )
+ end do
+ end if
+!!!JBK end
end if
! vertical transport of u
@@ -4789,11 +4902,30 @@ subroutine atm_compute_dyn_tend_work(nCells, nEdges, nVertices, nVertLevels_dumm
if(rk_step == 1) then
!DIR$ IVDEP
- do k=2,nVertLevels
- tend_w_euler(k,iCell) = tend_w_euler(k,iCell) - cqw(k,iCell)*( &
- rdzu(k)*(pp(k,iCell)-pp(k-1,iCell)) &
- - (fzm(k)*dpdz(k,iCell) + fzp(k)*dpdz(k-1,iCell)) ) ! dpdz is the buoyancy term here.
- end do
+
+!!!JBK - Mods for using log(p)
+
+! do k=2,nVertLevels
+! tend_w_euler(k,iCell) = tend_w_euler(k,iCell) - cqw(k,iCell)*( &
+! rdzu(k)*(pp(k,iCell)-pp(k-1,iCell)) &
+! - (fzm(k)*dpdz(k,iCell) + fzp(k)*dpdz(k-1,iCell)) ) ! dpdz is the buoyancy term here.
+! end do
+
+ if(reference_sounding_switch.eq.1) then
+ do k=2,nVertLevels
+ tend_w_euler(k,iCell) = tend_w_euler(k,iCell) - cqw(k,iCell) &
+ *(rdzu(k)*(pp(k,iCell)-pp(k-1,iCell)) &
+ - (fzm(k)*dpdz(k,iCell) + fzp(k)*dpdz(k-1,iCell)) )
+! dpdz is the buoyancy term here.
+ end do
+ else
+ do k=2,nVertLevels
+ tend_w_euler(k,iCell) = tend_w_euler(k,iCell) - (fzm(k)*pp(k,iCell)+fzp(k)*pp(k-1,iCell)) &
+ *(cqw(k,iCell)*rdzu(k)*(plog(k,iCell)-plog(k-1,iCell)) &
+ -(fzm(k)*dpdz(k,iCell) + fzp(k)*dpdz(k-1,iCell)) )
+ end do
+ end if
+!!!JBK end
end if
end do
diff --git a/src/core_atmosphere/mpas_atm_core.F b/src/core_atmosphere/mpas_atm_core.F
index 065e74c547..542e9e71b1 100644
--- a/src/core_atmosphere/mpas_atm_core.F
+++ b/src/core_atmosphere/mpas_atm_core.F
@@ -551,7 +551,7 @@ subroutine atm_mpas_init_block(dminfo, stream_manager, block, mesh, dt)
call atm_compute_mesh_scaling(mesh, block % configs)
- call atm_compute_damping_coefs(mesh, block % configs)
+ call atm_compute_damping_coefs(mesh, block % configs, nVertLevels)
!
! Set up mask fields used in limited-area simulations
@@ -1205,31 +1205,50 @@ subroutine atm_compute_signs(mesh)
end subroutine atm_compute_signs
- subroutine atm_compute_damping_coefs(mesh, configs)
+ subroutine atm_compute_damping_coefs(mesh, configs, nVertLevels)
implicit none
type (mpas_pool_type), intent(inout) :: mesh
type (mpas_pool_type), intent(in) :: configs
- integer :: iCell, k
- integer, pointer :: nCells, nVertLevels
+ integer :: iCell, k, nVertLevels
+ integer, pointer :: nCells
real (kind=RKIND), pointer :: config_xnutr, config_zd
real (kind=RKIND) :: z, zt, m1, pii
real (kind=RKIND), dimension(:,:), pointer :: dss, zgrid
real (kind=RKIND), dimension(:), pointer :: meshDensity
real (kind=RKIND) :: dx_scale_power
+!!!JBK-vepssm
+ real (kind=RKIND), dimension(:), pointer :: rdzw, etp, etm, ewp, ewm
+ real (kind=RKIND), dimension(nVertLevels) :: height_u_levels, epssm_coeff_u
+ real (kind=RKIND), dimension(nVertLevels+1) :: height_w_levels, epssm_coeff_w
+ real (kind=RKIND), pointer :: max_coeff, min_coeff, transition_lower_bound, transition_upper_bound
+ real (kind=RKIND) :: transition_width
+!!!JBK
+
m1 = -1.0
pii = acos(m1)
call mpas_pool_get_dimension(mesh, 'nCells', nCells)
- call mpas_pool_get_dimension(mesh, 'nVertLevels', nVertLevels)
+! call mpas_pool_get_dimension(mesh, 'nVertLevels', nVertLevels)
call mpas_pool_get_array(mesh, 'meshDensity', meshDensity)
call mpas_pool_get_array(mesh, 'dss', dss)
call mpas_pool_get_array(mesh, 'zgrid', zgrid)
+!!!JBK-vepssm
+ 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)
+!!!JBK
call mpas_pool_get_config(configs, 'config_zd', config_zd)
call mpas_pool_get_config(configs, 'config_xnutr', config_xnutr)
@@ -1246,6 +1265,44 @@ subroutine atm_compute_damping_coefs(mesh, configs)
end do
end do
+ transition_width = transition_upper_bound - transition_lower_bound
+
+! initialization for heights of coordinate at u and w levels
+
+ 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,
+
+ 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))
+ 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))
+ end do
+!!!JBK
+
end subroutine atm_compute_damping_coefs
diff --git a/src/core_init_atmosphere/Registry.xml b/src/core_init_atmosphere/Registry.xml
index e40f5c3722..da251e2395 100644
--- a/src/core_init_atmosphere/Registry.xml
+++ b/src/core_init_atmosphere/Registry.xml
@@ -99,6 +99,11 @@
description="projecting layer values to the interface, linear vertical interpolation or integral"
possible_values="'linear_interpolation' or 'layer_integral' (layer average value)"/>
+
+
@@ -540,6 +545,7 @@
+
@@ -739,6 +745,9 @@
+
+
diff --git a/src/core_init_atmosphere/mpas_init_atm_cases.F b/src/core_init_atmosphere/mpas_init_atm_cases.F
index e3e1ba56ea..b156dd6f2c 100644
--- a/src/core_init_atmosphere/mpas_init_atm_cases.F
+++ b/src/core_init_atmosphere/mpas_init_atm_cases.F
@@ -2005,7 +2005,10 @@ subroutine init_atm_case_mtn_wave(mesh, nCells, nVertLevels, state, diag, config
integer, dimension(nCells, 2) :: next_cell
logical, parameter :: terrain_smooth = .false.
-
+!!!JBK reference sounding
+ logical, pointer :: config_reference_sounding
+ integer, pointer :: reference_sounding_switch
+!!!
real (kind=RKIND), dimension(:), pointer :: xCell, yCell, zCell
real (kind=RKIND), dimension(:), pointer :: xEdge, yEdge, zEdge
real (kind=RKIND), dimension(:), pointer :: xVertex, yVertex, zVertex
@@ -2046,7 +2049,10 @@ subroutine init_atm_case_mtn_wave(mesh, nCells, nVertLevels, state, diag, config
call mpas_pool_get_config(configs, 'config_coef_3rd_order', config_coef_3rd_order)
call mpas_pool_get_config(configs, 'config_theta_adv_order', config_theta_adv_order)
call mpas_pool_get_config(configs, 'config_interface_projection', config_interface_projection)
-
+!!!JBK reference sounding switch
+ call mpas_pool_get_config(configs, 'config_reference_sounding', config_reference_sounding)
+ call mpas_pool_get_array(mesh, 'reference_sounding_switch', reference_sounding_switch)
+!!!
call mpas_pool_get_array(mesh, 'weightsOnEdge', weightsOnEdge)
call mpas_pool_get_array(mesh, 'nEdgesOnEdge', nEdgesOnEdge)
call mpas_pool_get_array(mesh, 'edgesOnEdge', edgesOnEdge)
@@ -2135,7 +2141,13 @@ subroutine init_atm_case_mtn_wave(mesh, nCells, nVertLevels, state, diag, config
call atm_initialize_advection_rk(mesh, nCells, nEdges, maxEdges, on_a_sphere, sphere_radius )
call atm_initialize_deformation_weights(mesh, nCells, on_a_sphere, sphere_radius)
-
+!!!JBK
+ if(config_reference_sounding) then
+ reference_sounding_switch = 1
+ else
+ reference_sounding_switch = 0
+ end if
+!!!
xnutr = 0.1
zd = 10500.
@@ -2380,6 +2392,21 @@ subroutine init_atm_case_mtn_wave(mesh, nCells, nVertLevels, state, diag, config
rr (k,i) = p (k,i)**(1./rcv)/((rgas/p0)*t (k,i)*zz(k,i))-rb(k,i)
cqw(k,i) = 1.
end do
+!
+!!!JBK Mods to remove reference state
+!
+ if(reference_sounding_switch.eq.0) then
+ do k=1,nz1
+ rr (k,i) = rb(k,i)+rr(k,i)
+ rb (k,i) = 0.
+ rtb(k,i) = 0.
+ pb (k,i) = 0.
+ tb (k,i) = 0.
+ end do
+ end if
+!
+!!!JBK
+!
end do
call mpas_log_write(' ***** base state sounding ***** ')