From f0b26e6ad63c678182f6c1bc0e2f51eb10389edd Mon Sep 17 00:00:00 2001 From: Joe Klemp Date: Mon, 27 Nov 2023 17:21:53 -0700 Subject: [PATCH 1/3] 1) Added code to allow epssm to vary with height according to a specified profile 2) Added a logical switch (reference_sounding) to specify whether the thermodynamic variables are computed as a perturbation from a specified reference sounding. If false, the full thermodynamic variables are used and the pressure gradients on the large time steps are computed using log(p). --- src/core_atmosphere/Registry.xml | 45 ++ .../dynamics/mpas_atm_time_integration.F | 489 ++++++++++++++++-- src/core_atmosphere/mpas_atm_core.F | 76 ++- .../mpas_init_atm_cases.F | 13 + 4 files changed, 564 insertions(+), 59 deletions(-) diff --git a/src/core_atmosphere/Registry.xml b/src/core_atmosphere/Registry.xml index 4e50bea204..009893c627 100644 --- a/src/core_atmosphere/Registry.xml +++ b/src/core_atmosphere/Registry.xml @@ -197,6 +197,11 @@ description="Whether to advect scalar fields" possible_values=".true. or .false."/> + + + + + + + + + + + + + + @@ -607,6 +636,10 @@ + + + + @@ -1353,6 +1386,18 @@ + + + + + + + + diff --git a/src/core_atmosphere/dynamics/mpas_atm_time_integration.F b/src/core_atmosphere/dynamics/mpas_atm_time_integration.F index 63b2ddd25e..e01e4e8555 100644 --- a/src/core_atmosphere/dynamics/mpas_atm_time_integration.F +++ b/src/core_atmosphere/dynamics/mpas_atm_time_integration.F @@ -94,6 +94,8 @@ 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 + !WCM_epssm +! real (kind=RKIND), dimension(:), allocatable :: ewp, ewm, etp, etm ! variable_epssm contains @@ -376,6 +378,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,6 +422,15 @@ subroutine atm_srk3(domain, dt, itimestep, exchange_halo_group) #ifdef DO_PHYSICS call mpas_pool_get_subpool(block % structs, 'diag_physics', diag_physics) #endif +!---------------------- +!SK_epssm: setting the 'mesh' pool pointer + call mpas_pool_get_subpool(domain % blocklist % structs, 'mesh', mesh) +!---------------------------- + +!------------------------ +!SK_epssm: added this to retrieve the pointer for rdzw + call mpas_pool_get_array(mesh, 'rdzw', rdzw) +!------------------ ! ! Retrieve dimensions @@ -487,6 +503,14 @@ subroutine atm_srk3(domain, dt, itimestep, exchange_halo_group) call mpas_pool_get_array(tend_physics, 'tend_ru_physics', tend_ru_physics) tend_ru_physics(:,nEdges+1) = 0.0_RKIND + ! WCS_epssm + ! allocate storage for epssm column arrays and set values +! allocate(ewp(nVertLevels+1)) +! allocate(ewm(nVertLevels+1)) +! allocate(etp(nVertLevels)) +! allocate(etm(nVertLevels)) +! call set_epssm(ewp,ewm,etp,etm,nVertLevels,rdzw) + ! ! Initialize RK weights ! @@ -1257,6 +1281,12 @@ subroutine atm_srk3(domain, dt, itimestep, exchange_halo_group) end if ! regional_MPAS addition + ! WCS_epssm +! deallocate(ewp) +! deallocate(ewm) +! deallocate(etp) +! deallocate(etm) + call summarize_timestep(domain) end subroutine atm_srk3 @@ -1603,6 +1633,15 @@ 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 + logical, pointer :: reference_sounding + 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_reference_sounding', reference_sounding) +!!!JBK call mpas_pool_get_config(configs, 'config_epssm', epssm) @@ -1640,6 +1679,9 @@ subroutine atm_compute_vert_imp_coefs(state, mesh, diag, configs, nVertLevels, d 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, & +!!!JBK cellSolveStart, cellSolveEnd, edgeSolveStart, edgeSolveEnd) @@ -1650,6 +1692,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, & +!!!JBK cellSolveStart, cellSolveEnd, edgeSolveStart, edgeSolveEnd) use mpas_atm_dimensions @@ -1688,6 +1733,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 + logical :: reference_sounding +!!!JBK ! @@ -1697,6 +1748,15 @@ 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 +! temporary print to check etm and ewm + + integer :: iprint=0 + if(iprint.eq.0) then + do k=1,nVertLevels + call mpas_log_write(' k, etp(k), ewp(k) = $i $r $r', intArgs=(/k/), realArgs=(/etp(k),ewp(k)/)) + end do + iprint = 1 + end if ! set coefficients dtseps = .5*dts*(1.+epssm) @@ -1705,24 +1765,37 @@ 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)) + + ! WCS_epssm + ! 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)) + ! WCS_epssm + ! 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 full variables and log(p) + if(reference_sounding) then +!!!JBK do k=1,nVertLevels ! qtotal = 0. @@ -1731,10 +1804,22 @@ 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)) + ! WCS_epssm + ! 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 +1829,45 @@ 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 ) & + + ! WCS_epssm + a_tri(k,iCell) = a_tri(k,iCell)*etp(k-1)*ewp(k-1) + + ! WiCS_epssm + ! 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)) + 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)) + ! WCS_epssm + 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 ) + + ! WCS_epssm + c_tri(k) = c_tri(k)*etp(k)*ewp(k+1) end do + ! WCS_epssm + 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) + ! WCS_epssm + ! 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 +2108,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 +2190,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 +2204,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 +2276,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 ! @@ -2169,6 +2288,16 @@ subroutine atm_advance_acoustic_step_work(nCells, nEdges, nCellsSolve, cellStart real (kind=RKIND) :: pgrad, flux, resm, rdts +! temporary print to check etm and ewm + + integer :: iprint=0 + if(iprint.eq.0) then + do k=1,nVertLevels + call mpas_log_write(' k, etm(k), ewm(k) = $i $r $r', intArgs=(/k/), realArgs=(/etm(k),ewm(k)/)) + end do + iprint = 1 + end if + rcv = rgas / (cp - rgas) c2 = cp * rcv resm = (1.0 - epssm) / (1.0 + epssm) @@ -2281,36 +2410,60 @@ 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)) + ! WCS_epssm + ! 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)) + 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)) + ! WCS_epssm + ! 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)) + 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) + ! WCS_epssm + ! 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)) + ! WCS_epssm + ! 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) + ! WCS_epssm + ! 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 +2485,23 @@ 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) + ! WCS_epssm + ! 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)) + ! WCS_epssm + ! 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 +2510,9 @@ 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) + ! WCS_epssm + ! 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 @@ -3928,6 +4090,9 @@ subroutine atm_compute_dyn_tend(tend, tend_physics, state, diag, mesh, configs, 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 +!!!JBK - Mod for reference sounding switch + logical, pointer :: reference_sounding +!!!JBK end call mpas_pool_get_config(mesh, 'sphere_radius', r_earth) @@ -3949,7 +4114,9 @@ subroutine atm_compute_dyn_tend(tend, tend_physics, state, diag, mesh, configs, call mpas_pool_get_config(configs, 'config_rayleigh_damp_u_timescale_days', config_rayleigh_damp_u_timescale_days) call mpas_pool_get_config(configs, 'config_number_rayleigh_damp_u_levels', config_number_rayleigh_damp_u_levels) call mpas_pool_get_config(configs, 'config_number_cam_damping_levels', config_number_cam_damping_levels) - +!!!JBK Mod foe reference sounding switch + call mpas_pool_get_config(configs, 'config_reference_sounding', reference_sounding) +!!!JBK end call mpas_pool_get_array(state, 'rho_zz', rho_zz, 2) call mpas_pool_get_array(state, 'u', u, 2) call mpas_pool_get_array(state, 'w', w, 2) @@ -4072,6 +4239,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 switch + reference_sounding, & +!!!JBK end rthdynten, & cellStart, cellEnd, vertexStart, vertexEnd, edgeStart, edgeEnd, & cellSolveStart, cellSolveEnd, vertexSolveStart, vertexSolveEnd, edgeSolveStart, edgeSolveEnd) @@ -4096,6 +4266,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 switch + reference_sounding, & +!!!JBK end rthdynten, & cellStart, cellEnd, vertexStart, vertexEnd, edgeStart, edgeEnd, & cellSolveStart, cellSolveEnd, vertexSolveStart, vertexSolveEnd, edgeSolveStart, edgeSolveEnd) @@ -4253,6 +4427,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 + logical :: reference_sounding +! 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 +4449,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 = $l ', logicArgs=(/reference_sounding/)) +!!!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 +4546,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) 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 +4585,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) 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 +5008,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) 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 @@ -6604,4 +6842,147 @@ subroutine summarize_timestep(domain) end subroutine summarize_timestep +!---------------------- +!----------------------- + subroutine set_epssm(ewp,ewm,etp,etm,n,rdzw) + implicit none + integer :: n + real (kind=RKIND), dimension(n+1) :: ewp, ewm + real (kind=RKIND), dimension(n ) :: etp, etm +!----------------------- +!SK_epssm: added rdzw to the function call + real (kind=RKIND), dimension(n) :: rdzw + +!---------------------------------- +!SK_epssm: created local variables + real(kind=RKIND), dimension(n):: height_u_levels + real(kind=RKIND), dimension(n+1):: height_w_levels + real (kind=RKIND), dimension(n) :: epssm_coeff_u + real (kind=RKIND), dimension(n+1):: epssm_coeff_w + real (kind=RKIND) :: max_coeff, min_coeff, transition_width, transition_point + integer :: ilevel +!-------------------- +! SK_epssm: commented this out since we are transitiong to height based epssm +! real (kind=RKIND), parameter :: epssm_b = 0.1 + +!------------------ +!SK_epssm: intialization + epssm_coeff_u = 0.0_RKIND + epssm_coeff_w = 0.0_RKIND + height_u_levels = 0.0_RKIND + height_w_levels= 0.0_RKIND + +!!! transition_point = 60000.0_RKIND ! Changed the value for WACCM simulations + transition_point = 10500.0_RKIND ! (Model top 21K) ! Height of center of transition zone +!!! transition_width =20000.0_RKIND ! Changed the value for WACCM simulations + transition_width = 5000.0_RKIND ! (Model top 21K) ! Width of transition zone + min_coeff = 0.1_RKIND ! Value of epssm below transition zone + max_coeff = 0.5_RKIND ! Value of epssm above transition zone +!--------------------- +!--------------------- +!SK_epssm +! calculate height at levels + call calculate_computational_height(n, rdzw, height_u_levels, height_w_levels) +! call the ramping function + call calculate_coef_height(n,height_u_levels, min_coeff,& + & max_coeff, transition_width, transition_point,epssm_coeff_u) + call calculate_coef_height(n+1,height_w_levels, min_coeff,& + & max_coeff, transition_width, transition_point,epssm_coeff_w) +!-------------------------------- +!SK_epssm: commented the lines below so that we can have epssm varying with height +! etp(:) = 0.5*(1.0+epssm_b) +! etm(:) = 0.5*(1.0-epssm_b) +! ewp(:) = 0.5*(1.0+epssm_b) +! ewm(:) = 0.5*(1.0-epssm_b) +! SK_epssm: Not sure about the epssm_coeff value I am using for ewp and ewm + do ilevel =1, n + etp(ilevel) = 0.5*(1.0 + epssm_coeff_u(ilevel)) + etm(ilevel) = 0.5*(1.0 - epssm_coeff_u(ilevel)) + ewp(ilevel) = 0.5*(1.0 + epssm_coeff_w(ilevel)) + ewm(ilevel) = 0.5*(1.0 - epssm_coeff_w(ilevel)) + + call mpas_log_write('ilevel, height_u_levels, etp, etm, height_w_levels, ewp, ewm $i $r $r $r $r $r $r', intARgs=(/ilevel/), & + realArgs=(/height_u_levels(ilevel),etp(ilevel),etm(ilevel), & + height_w_levels(ilevel),ewp(ilevel),ewm(ilevel)/)) + enddo + ewp(n+1) = 0.5*(1.0 + epssm_coeff_w(n+1)) + ewm(n+1) = 0.5*(1.0 - epssm_coeff_w(n+1)) + + end subroutine set_epssm + +!---------------- +!SK_epssm + subroutine calculate_computational_height(nVertLevels,rdzw,height_u_levels, height_w_levels) + integer, intent(in) :: nVertLevels + real (kind=RKIND), dimension( nVertLevels ), intent(in) :: rdzw + real (kind=RKIND), dimension( nVertLevels ), intent(inout) :: height_u_levels + real (kind=RKIND), dimension( nVertLevels+1 ), intent(inout) :: height_w_levels + +!define variables + integer k +!initialization + height_w_levels(:) = 0.0_RKIND + +!calculate the computational height at u and w levels +!DIR$ IVDEP + 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)) +! call mpas_log_write('test $i $r', intARgs=(/k/),realArgs=(/rdzw(k)/)) + enddo + +! do k =1, nVertLevels +! call mpas_log_write('test1 $i $r $r $r', intARgs=(/k/),realArgs=(/height_u_levels(k),& +! height_w_levels(k),height_w_levels(k+1)/)) +! enddo + + end subroutine calculate_computational_height + +! ramp a value base on height + subroutine calculate_coef_height(nVertLevels,height_levels, min_coeff,& + & max_coeff, transition_width, transition_point,damp_coeff) + integer, intent(in)::nVertLevels + real (kind=RKIND) :: max_coeff, min_coeff, transition_width, transition_point + real (kind=RKIND), dimension(nVertLevels), intent (in) :: height_levels + real (kind=RKIND), dimension(nVertLevels), intent (out) :: damp_coeff +! local parameters + integer :: k, counter_transition + real (kind=RKIND) :: dx,x, transition_lower_bound, transition_upper_bound +! initialization + damp_coeff = 0.0_RKIND + counter_transition = 0 + dx = 0.0_RKIND + x = -pii/2.0_RKIND + transition_lower_bound = transition_point - transition_width*0.5_RKIND + transition_upper_bound = transition_point + transition_width*0.5_RKIND + +! calculate how many levels are in the transition region + do k =1, nVertLevels + if(((height_levels(k)-transition_lower_bound)>= 0.0_RKIND).and.((height_levels(k)-transition_upper_bound)<=0.0_RKIND)) then + counter_transition = counter_transition +1 + endif + enddo + dx = pii/REAL(counter_transition-1) +! call mpas_log_write('check $i $r $r $r $r', intARgs=(/counter_transition/),realArgs=(/dx,x,transition_lower_bound,transition_upper_bound/)) +! start calculating the coefficient + do k =1, nVertLevels + if((height_levels(k)-transition_lower_bound) < 0.0_RKIND ) then + damp_coeff(k)= min_coeff + elseif((height_levels(k)-transition_upper_bound) > 0.0_RKIND ) then + damp_coeff(k) =max_coeff + else +!!!JK_epssm +! alternative calculation of damp_coeff in transition zone +! damp_coeff(k) = min_coeff+ (sin(x)+1)*0.5_RKIND*(max_coeff-min_coeff) +! x = x + dx +! + x = (height_levels(k)-transition_lower_bound)/transition_width + damp_coeff(k) = min_coeff + sin(0.5_RKIND*pii*x)**2*(max_coeff-min_coeff) +!!!JK_epssm end + endif +! call mpas_log_write('damping factor $i $r $r', intARgs=(/k/),realArgs=(/damp_coeff(k), height_levels(k)/)) + enddo + + end subroutine calculate_coef_height + end module atm_time_integration diff --git a/src/core_atmosphere/mpas_atm_core.F b/src/core_atmosphere/mpas_atm_core.F index 065e74c547..0a59a02fe5 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_min_coeff', min_coeff) + call mpas_pool_get_config(configs, 'config_max_coeff', max_coeff) + call mpas_pool_get_config(configs, 'config_transition_lower_bound', transition_lower_bound) + call mpas_pool_get_config(configs, 'config_transition_upper_bound', 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,53 @@ subroutine atm_compute_damping_coefs(mesh, configs) end do end do +!!!JBK-vepssm - code to compute variable epssm profiles at model startup + +! min_coeff = 0.1_RKIND +! max_coeff = 0.5_RKIND +! transition_lower_bound = 8000._RKIND +! transition_upper_bound = 13000._RKIND + 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)) + +! call mpas_log_write('test $i $r', intARgs=(/k/),realArgs=(/rdzw(k)/)) + + enddo + +! Height dependent values of epssm + + 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/mpas_init_atm_cases.F b/src/core_init_atmosphere/mpas_init_atm_cases.F index e3e1ba56ea..e76f3c81bc 100644 --- a/src/core_init_atmosphere/mpas_init_atm_cases.F +++ b/src/core_init_atmosphere/mpas_init_atm_cases.F @@ -2380,6 +2380,19 @@ 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 removed reference state +! + 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 +! +!!!JBK +! end do call mpas_log_write(' ***** base state sounding ***** ') From 191470911185a2254002bc727b80b64e9994551f Mon Sep 17 00:00:00 2001 From: Joseph Klemp Date: Wed, 6 Dec 2023 11:20:21 -0700 Subject: [PATCH 2/3] Modifications to clean version 8 of MPAS coe: 1) Code added to ../src/core_atmosphere/mpas_atm_core.F to compute vertical profiles of epssm that vary with height at the start of the time integration. These profiles are determined by the parameters config_min_coeff, config_max_coeff, config_transition_lower_bound, and config_transition_upper_bound, and are specified in namelist.atmosphere. 2) Code added to allow thermodynamic variables to be represented either as perurbations from a specified reference state or as the full variables without a reference state. This option is controlled by the logical config_reference_sounding as specified in namelist.init_atmosphere. When configured using the full variables, the horizontal and vertical pressure gradients computed on the large time steps are written in terms of log(p). --- src/core_atmosphere/Registry.xml | 13 ++-- .../dynamics/mpas_atm_time_integration.F | 60 +++++++++---------- src/core_atmosphere/mpas_atm_core.F | 2 +- src/core_init_atmosphere/Registry.xml | 9 +++ .../mpas_init_atm_cases.F | 36 +++++++---- 5 files changed, 70 insertions(+), 50 deletions(-) diff --git a/src/core_atmosphere/Registry.xml b/src/core_atmosphere/Registry.xml index 009893c627..5057db80f1 100644 --- a/src/core_atmosphere/Registry.xml +++ b/src/core_atmosphere/Registry.xml @@ -197,11 +197,6 @@ description="Whether to advect scalar fields" possible_values=".true. or .false."/> - - - + + + #ifdef MPAS_CAM_DYCORE @@ -656,6 +653,7 @@ + #ifdef MPAS_CAM_DYCORE @@ -1445,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 e01e4e8555..1bde629720 100644 --- a/src/core_atmosphere/dynamics/mpas_atm_time_integration.F +++ b/src/core_atmosphere/dynamics/mpas_atm_time_integration.F @@ -1635,16 +1635,14 @@ subroutine atm_compute_vert_imp_coefs(state, mesh, diag, configs, nVertLevels, d !!!JBK-vepssm and reference sounding real (kind=RKIND), dimension(:), pointer :: etp, etm, ewp, ewm - logical, pointer :: reference_sounding + 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_config(configs, 'config_reference_sounding', reference_sounding) -!!!JBK - + 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) @@ -1674,13 +1672,12 @@ 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, & + etp, etm, ewp, ewm, reference_sounding_switch, & !!!JBK cellSolveStart, cellSolveEnd, edgeSolveStart, edgeSolveEnd) @@ -1693,7 +1690,7 @@ subroutine atm_compute_vert_imp_coefs_work(nCells, moist_start, moist_end, dts, a_tri, alpha_tri, gamma_tri, cofrz, rdzw, fzm, fzp, rdzu, scalars, & cellStart, cellEnd, edgeStart, edgeEnd, & !!!JBK - etp, etm, ewp, ewm, reference_sounding, & + etp, etm, ewp, ewm, reference_sounding_switch, & !!!JBK cellSolveStart, cellSolveEnd, edgeSolveStart, edgeSolveEnd) @@ -1737,7 +1734,7 @@ subroutine atm_compute_vert_imp_coefs_work(nCells, moist_start, moist_end, dts, 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 - logical :: reference_sounding + integer :: reference_sounding_switch !!!JBK @@ -1793,8 +1790,8 @@ subroutine atm_compute_vert_imp_coefs_work(nCells, moist_start, moist_end, dts, coftz(nVertLevels+1,iCell) = 0.0 !DIR$ IVDEP -!!!JBK Mod to use full variables and log(p) - if(reference_sounding) then +!!!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 @@ -4087,14 +4084,13 @@ 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 -!!!JBK - Mod for reference sounding switch - logical, pointer :: reference_sounding + 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) call mpas_pool_get_config(configs, 'config_mix_full', config_mix_full) @@ -4114,9 +4110,7 @@ subroutine atm_compute_dyn_tend(tend, tend_physics, state, diag, mesh, configs, call mpas_pool_get_config(configs, 'config_rayleigh_damp_u_timescale_days', config_rayleigh_damp_u_timescale_days) call mpas_pool_get_config(configs, 'config_number_rayleigh_damp_u_levels', config_number_rayleigh_damp_u_levels) call mpas_pool_get_config(configs, 'config_number_cam_damping_levels', config_number_cam_damping_levels) -!!!JBK Mod foe reference sounding switch - call mpas_pool_get_config(configs, 'config_reference_sounding', reference_sounding) -!!!JBK end + call mpas_pool_get_array(state, 'rho_zz', rho_zz, 2) call mpas_pool_get_array(state, 'u', u, 2) call mpas_pool_get_array(state, 'w', w, 2) @@ -4146,7 +4140,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) @@ -4239,8 +4235,8 @@ 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 switch - reference_sounding, & +!!!JBK - Mod for reference sounding + reference_sounding_switch, & !!!JBK end rthdynten, & cellStart, cellEnd, vertexStart, vertexEnd, edgeStart, edgeEnd, & @@ -4267,8 +4263,8 @@ subroutine atm_compute_dyn_tend_work(nCells, nEdges, nVertices, nVertLevels_dumm 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 switch - reference_sounding, & +!!!JBK - Mod for reference sounding + reference_sounding_switch, & !!!JBK end rthdynten, & cellStart, cellEnd, vertexStart, vertexEnd, edgeStart, edgeEnd, & @@ -4428,11 +4424,11 @@ subroutine atm_compute_dyn_tend_work(nCells, nEdges, nVertices, nVertLevels_dumm real (kind=RKIND) :: q_im2, q_im1, q_i, q_ip1, ua, coef3 !!!JBK - mod for reference sounding switch - logical :: reference_sounding -! Declaration for plog + integer :: reference_sounding_switch +! Declaration for plog real (kind=RKIND) :: p0 real (kind=RKIND), dimension(nVertLevels,nCells+1) :: plog -!!!JBK end +!!!JBK end flux4(q_im2, q_im1, q_i, q_ip1, ua) = & ua*( 7.*(q_i + q_im1) - (q_ip1 + q_im2) )/12.0 @@ -4449,10 +4445,10 @@ 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... + p0 = 1.0e+05 ! this should come from somewhere else... !!!JBK - Temporary mod to print reference_sounding - call mpas_log_write('reference_sounding = $l ', logicArgs=(/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 @@ -4556,7 +4552,7 @@ subroutine atm_compute_dyn_tend_work(nCells, nEdges, nVertices, nVertLevels_dumm 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) end do - if(reference_sounding) then + 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 @@ -4591,7 +4587,7 @@ subroutine atm_compute_dyn_tend_work(nCells, nEdges, nVertices, nVertLevels_dumm ! 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) then + 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))) & @@ -5017,7 +5013,7 @@ subroutine atm_compute_dyn_tend_work(nCells, nEdges, nVertices, nVertLevels_dumm ! - (fzm(k)*dpdz(k,iCell) + fzp(k)*dpdz(k-1,iCell)) ) ! dpdz is the buoyancy term here. ! end do - if(reference_sounding) then + 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)) & diff --git a/src/core_atmosphere/mpas_atm_core.F b/src/core_atmosphere/mpas_atm_core.F index 0a59a02fe5..2a0129f200 100644 --- a/src/core_atmosphere/mpas_atm_core.F +++ b/src/core_atmosphere/mpas_atm_core.F @@ -1284,7 +1284,7 @@ subroutine atm_compute_damping_coefs(mesh, configs, nVertLevels) enddo -! Height dependent values of epssm +! 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 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 e76f3c81bc..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. @@ -2381,15 +2393,17 @@ subroutine init_atm_case_mtn_wave(mesh, nCells, nVertLevels, state, diag, config cqw(k,i) = 1. end do ! -!!!JBK Mods to removed reference state +!!!JBK Mods to remove reference state ! - 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 + 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 ! From 24e589dd09a2bc5296fa8b4355690a451c8d3979 Mon Sep 17 00:00:00 2001 From: Joseph Klemp Date: Thu, 28 Dec 2023 12:46:42 -0700 Subject: [PATCH 3/3] Revisions to the pull request that contain additions to the clean MPAS version 8 to include: 1) an option to use full thermodynamic variables instead of perturbations from a reference profile and 2) a capability for specifying an epssm profile that is variable with height. These revisions include the changes suggested by Bill Skamarock in his review of the pull request. --- src/core_atmosphere/Registry.xml | 20 +- .../dynamics/mpas_atm_time_integration.F | 245 ------------------ src/core_atmosphere/mpas_atm_core.F | 17 +- 3 files changed, 14 insertions(+), 268 deletions(-) diff --git a/src/core_atmosphere/Registry.xml b/src/core_atmosphere/Registry.xml index 5057db80f1..a8f08f9c96 100644 --- a/src/core_atmosphere/Registry.xml +++ b/src/core_atmosphere/Registry.xml @@ -265,23 +265,23 @@ description="Maximum w-damping coefficient at model top" possible_values="0 $\leq$ config_xnutr $\leq$ 1"/> - - - - @@ -1385,16 +1385,16 @@ description="Reciprocal dzw"/> + description="etp = (1+epssm(z))/2 at theta point"/> + description="etm = (1-epssm(z))/2 at theta point"/> + description="ewp = (1+epssm(z))/2 at w point"/> + description="ewm = (1-epssm(z))/2 at w point"/> diff --git a/src/core_atmosphere/dynamics/mpas_atm_time_integration.F b/src/core_atmosphere/dynamics/mpas_atm_time_integration.F index 1bde629720..337489f9c1 100644 --- a/src/core_atmosphere/dynamics/mpas_atm_time_integration.F +++ b/src/core_atmosphere/dynamics/mpas_atm_time_integration.F @@ -94,9 +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 - !WCM_epssm -! real (kind=RKIND), dimension(:), allocatable :: ewp, ewm, etp, etm ! variable_epssm - contains @@ -422,16 +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 -!---------------------- -!SK_epssm: setting the 'mesh' pool pointer - call mpas_pool_get_subpool(domain % blocklist % structs, 'mesh', mesh) -!---------------------------- - -!------------------------ -!SK_epssm: added this to retrieve the pointer for rdzw - call mpas_pool_get_array(mesh, 'rdzw', rdzw) -!------------------ - ! ! Retrieve dimensions ! Note: nCellsSolve and nVerticesSolve are not currently used in this function @@ -503,14 +490,6 @@ subroutine atm_srk3(domain, dt, itimestep, exchange_halo_group) call mpas_pool_get_array(tend_physics, 'tend_ru_physics', tend_ru_physics) tend_ru_physics(:,nEdges+1) = 0.0_RKIND - ! WCS_epssm - ! allocate storage for epssm column arrays and set values -! allocate(ewp(nVertLevels+1)) -! allocate(ewm(nVertLevels+1)) -! allocate(etp(nVertLevels)) -! allocate(etm(nVertLevels)) -! call set_epssm(ewp,ewm,etp,etm,nVertLevels,rdzw) - ! ! Initialize RK weights ! @@ -1281,12 +1260,6 @@ subroutine atm_srk3(domain, dt, itimestep, exchange_halo_group) end if ! regional_MPAS addition - ! WCS_epssm -! deallocate(ewp) -! deallocate(ewm) -! deallocate(etp) -! deallocate(etm) - call summarize_timestep(domain) end subroutine atm_srk3 @@ -1745,16 +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 -! temporary print to check etm and ewm - - integer :: iprint=0 - if(iprint.eq.0) then - do k=1,nVertLevels - call mpas_log_write(' k, etp(k), ewp(k) = $i $r $r', intArgs=(/k/), realArgs=(/etp(k),ewp(k)/)) - end do - iprint = 1 - end if - ! set coefficients dtseps = .5*dts*(1.+epssm) rcv = rgas/(cp-rgas) @@ -1771,18 +1734,11 @@ subroutine atm_compute_vert_imp_coefs_work(nCells, moist_start, moist_end, dts, !DIR$ IVDEP do k=2,nVertLevels - - ! WCS_epssm - ! 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 - ! WCS_epssm - ! 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)) @@ -1801,9 +1757,6 @@ subroutine atm_compute_vert_imp_coefs_work(nCells, moist_start, moist_end, dts, ! end do qtotal = qtot(k,iCell) - ! WCS_epssm - ! 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)) @@ -1829,40 +1782,24 @@ 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) - - ! WCS_epssm a_tri(k,iCell) = a_tri(k,iCell)*etp(k-1)*ewp(k-1) - ! WiCS_epssm - ! 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)) 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)) - ! WCS_epssm 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 ) - - ! WCS_epssm c_tri(k) = c_tri(k)*etp(k)*ewp(k+1) end do - ! WCS_epssm c_tri(nVertLevels) = 0.0 !MGD VECTOR DEPENDENCE do k=2,nVertLevels - ! WCS_epssm - ! 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 @@ -2285,16 +2222,6 @@ subroutine atm_advance_acoustic_step_work(nCells, nEdges, nCellsSolve, cellStart real (kind=RKIND) :: pgrad, flux, resm, rdts -! temporary print to check etm and ewm - - integer :: iprint=0 - if(iprint.eq.0) then - do k=1,nVertLevels - call mpas_log_write(' k, etm(k), ewm(k) = $i $r $r', intArgs=(/k/), realArgs=(/etm(k),ewm(k)/)) - end do - iprint = 1 - end if - rcv = rgas / (cp - rgas) c2 = cp * rcv resm = (1.0 - epssm) / (1.0 + epssm) @@ -2407,15 +2334,8 @@ subroutine atm_advance_acoustic_step_work(nCells, nEdges, nCellsSolve, cellStart !DIR$ IVDEP do k=1, nVertLevels - ! WCS_epssm - ! 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)) 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)) - ! WCS_epssm - ! 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)) 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)) @@ -2423,23 +2343,11 @@ subroutine atm_advance_acoustic_step_work(nCells, nEdges, nCellsSolve, cellStart !DIR$ IVDEP do k=2, nVertLevels - ! WCS_epssm - ! 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 - ! WCS_epssm - ! 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)) & @@ -2457,8 +2365,6 @@ subroutine atm_advance_acoustic_step_work(nCells, nEdges, nCellsSolve, cellStart !MGD VECTOR DEPENDENCE do k=2,nVertLevels - ! WCS_epssm - ! 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 @@ -2482,8 +2388,6 @@ 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 - ! WCS_epssm - ! 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 @@ -2491,10 +2395,6 @@ subroutine atm_advance_acoustic_step_work(nCells, nEdges, nCellsSolve, cellStart !DIR$ IVDEP do k=1,nVertLevels - ! WCS_epssm - ! 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) & @@ -2507,8 +2407,6 @@ 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) - ! WCS_epssm - ! 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 @@ -6838,147 +6736,4 @@ subroutine summarize_timestep(domain) end subroutine summarize_timestep -!---------------------- -!----------------------- - subroutine set_epssm(ewp,ewm,etp,etm,n,rdzw) - implicit none - integer :: n - real (kind=RKIND), dimension(n+1) :: ewp, ewm - real (kind=RKIND), dimension(n ) :: etp, etm -!----------------------- -!SK_epssm: added rdzw to the function call - real (kind=RKIND), dimension(n) :: rdzw - -!---------------------------------- -!SK_epssm: created local variables - real(kind=RKIND), dimension(n):: height_u_levels - real(kind=RKIND), dimension(n+1):: height_w_levels - real (kind=RKIND), dimension(n) :: epssm_coeff_u - real (kind=RKIND), dimension(n+1):: epssm_coeff_w - real (kind=RKIND) :: max_coeff, min_coeff, transition_width, transition_point - integer :: ilevel -!-------------------- -! SK_epssm: commented this out since we are transitiong to height based epssm -! real (kind=RKIND), parameter :: epssm_b = 0.1 - -!------------------ -!SK_epssm: intialization - epssm_coeff_u = 0.0_RKIND - epssm_coeff_w = 0.0_RKIND - height_u_levels = 0.0_RKIND - height_w_levels= 0.0_RKIND - -!!! transition_point = 60000.0_RKIND ! Changed the value for WACCM simulations - transition_point = 10500.0_RKIND ! (Model top 21K) ! Height of center of transition zone -!!! transition_width =20000.0_RKIND ! Changed the value for WACCM simulations - transition_width = 5000.0_RKIND ! (Model top 21K) ! Width of transition zone - min_coeff = 0.1_RKIND ! Value of epssm below transition zone - max_coeff = 0.5_RKIND ! Value of epssm above transition zone -!--------------------- -!--------------------- -!SK_epssm -! calculate height at levels - call calculate_computational_height(n, rdzw, height_u_levels, height_w_levels) -! call the ramping function - call calculate_coef_height(n,height_u_levels, min_coeff,& - & max_coeff, transition_width, transition_point,epssm_coeff_u) - call calculate_coef_height(n+1,height_w_levels, min_coeff,& - & max_coeff, transition_width, transition_point,epssm_coeff_w) -!-------------------------------- -!SK_epssm: commented the lines below so that we can have epssm varying with height -! etp(:) = 0.5*(1.0+epssm_b) -! etm(:) = 0.5*(1.0-epssm_b) -! ewp(:) = 0.5*(1.0+epssm_b) -! ewm(:) = 0.5*(1.0-epssm_b) -! SK_epssm: Not sure about the epssm_coeff value I am using for ewp and ewm - do ilevel =1, n - etp(ilevel) = 0.5*(1.0 + epssm_coeff_u(ilevel)) - etm(ilevel) = 0.5*(1.0 - epssm_coeff_u(ilevel)) - ewp(ilevel) = 0.5*(1.0 + epssm_coeff_w(ilevel)) - ewm(ilevel) = 0.5*(1.0 - epssm_coeff_w(ilevel)) - - call mpas_log_write('ilevel, height_u_levels, etp, etm, height_w_levels, ewp, ewm $i $r $r $r $r $r $r', intARgs=(/ilevel/), & - realArgs=(/height_u_levels(ilevel),etp(ilevel),etm(ilevel), & - height_w_levels(ilevel),ewp(ilevel),ewm(ilevel)/)) - enddo - ewp(n+1) = 0.5*(1.0 + epssm_coeff_w(n+1)) - ewm(n+1) = 0.5*(1.0 - epssm_coeff_w(n+1)) - - end subroutine set_epssm - -!---------------- -!SK_epssm - subroutine calculate_computational_height(nVertLevels,rdzw,height_u_levels, height_w_levels) - integer, intent(in) :: nVertLevels - real (kind=RKIND), dimension( nVertLevels ), intent(in) :: rdzw - real (kind=RKIND), dimension( nVertLevels ), intent(inout) :: height_u_levels - real (kind=RKIND), dimension( nVertLevels+1 ), intent(inout) :: height_w_levels - -!define variables - integer k -!initialization - height_w_levels(:) = 0.0_RKIND - -!calculate the computational height at u and w levels -!DIR$ IVDEP - 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)) -! call mpas_log_write('test $i $r', intARgs=(/k/),realArgs=(/rdzw(k)/)) - enddo - -! do k =1, nVertLevels -! call mpas_log_write('test1 $i $r $r $r', intARgs=(/k/),realArgs=(/height_u_levels(k),& -! height_w_levels(k),height_w_levels(k+1)/)) -! enddo - - end subroutine calculate_computational_height - -! ramp a value base on height - subroutine calculate_coef_height(nVertLevels,height_levels, min_coeff,& - & max_coeff, transition_width, transition_point,damp_coeff) - integer, intent(in)::nVertLevels - real (kind=RKIND) :: max_coeff, min_coeff, transition_width, transition_point - real (kind=RKIND), dimension(nVertLevels), intent (in) :: height_levels - real (kind=RKIND), dimension(nVertLevels), intent (out) :: damp_coeff -! local parameters - integer :: k, counter_transition - real (kind=RKIND) :: dx,x, transition_lower_bound, transition_upper_bound -! initialization - damp_coeff = 0.0_RKIND - counter_transition = 0 - dx = 0.0_RKIND - x = -pii/2.0_RKIND - transition_lower_bound = transition_point - transition_width*0.5_RKIND - transition_upper_bound = transition_point + transition_width*0.5_RKIND - -! calculate how many levels are in the transition region - do k =1, nVertLevels - if(((height_levels(k)-transition_lower_bound)>= 0.0_RKIND).and.((height_levels(k)-transition_upper_bound)<=0.0_RKIND)) then - counter_transition = counter_transition +1 - endif - enddo - dx = pii/REAL(counter_transition-1) -! call mpas_log_write('check $i $r $r $r $r', intARgs=(/counter_transition/),realArgs=(/dx,x,transition_lower_bound,transition_upper_bound/)) -! start calculating the coefficient - do k =1, nVertLevels - if((height_levels(k)-transition_lower_bound) < 0.0_RKIND ) then - damp_coeff(k)= min_coeff - elseif((height_levels(k)-transition_upper_bound) > 0.0_RKIND ) then - damp_coeff(k) =max_coeff - else -!!!JK_epssm -! alternative calculation of damp_coeff in transition zone -! damp_coeff(k) = min_coeff+ (sin(x)+1)*0.5_RKIND*(max_coeff-min_coeff) -! x = x + dx -! - x = (height_levels(k)-transition_lower_bound)/transition_width - damp_coeff(k) = min_coeff + sin(0.5_RKIND*pii*x)**2*(max_coeff-min_coeff) -!!!JK_epssm end - endif -! call mpas_log_write('damping factor $i $r $r', intARgs=(/k/),realArgs=(/damp_coeff(k), height_levels(k)/)) - enddo - - end subroutine calculate_coef_height - end module atm_time_integration diff --git a/src/core_atmosphere/mpas_atm_core.F b/src/core_atmosphere/mpas_atm_core.F index 2a0129f200..542e9e71b1 100644 --- a/src/core_atmosphere/mpas_atm_core.F +++ b/src/core_atmosphere/mpas_atm_core.F @@ -1244,10 +1244,10 @@ subroutine atm_compute_damping_coefs(mesh, configs, nVertLevels) 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_min_coeff', min_coeff) - call mpas_pool_get_config(configs, 'config_max_coeff', max_coeff) - call mpas_pool_get_config(configs, 'config_transition_lower_bound', transition_lower_bound) - call mpas_pool_get_config(configs, 'config_transition_upper_bound', transition_upper_bound) + 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) @@ -1265,12 +1265,6 @@ subroutine atm_compute_damping_coefs(mesh, configs, nVertLevels) end do end do -!!!JBK-vepssm - code to compute variable epssm profiles at model startup - -! min_coeff = 0.1_RKIND -! max_coeff = 0.5_RKIND -! transition_lower_bound = 8000._RKIND -! transition_upper_bound = 13000._RKIND transition_width = transition_upper_bound - transition_lower_bound ! initialization for heights of coordinate at u and w levels @@ -1279,9 +1273,6 @@ subroutine atm_compute_damping_coefs(mesh, configs, nVertLevels) 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)) - -! call mpas_log_write('test $i $r', intARgs=(/k/),realArgs=(/rdzw(k)/)) - enddo ! Height dependent values of epssm; profiles stored in etp, etm, ewp, and ewm,