From 19219b028563e51e8f5599d79789743b8375c904 Mon Sep 17 00:00:00 2001 From: Hing Ong Date: Thu, 21 Nov 2024 21:57:31 -0700 Subject: [PATCH 1/6] Add a namelist option for deep atmosphere --- src/core_atmosphere/Registry.xml | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/core_atmosphere/Registry.xml b/src/core_atmosphere/Registry.xml index 6378596797..ddcc94d159 100644 --- a/src/core_atmosphere/Registry.xml +++ b/src/core_atmosphere/Registry.xml @@ -242,6 +242,11 @@ description="Scale eddy viscosities with mesh-density function for horizontal diffusion" possible_values=".true. or .false."/> + + Date: Fri, 22 Nov 2024 05:19:41 -0700 Subject: [PATCH 2/6] Add nondimensional radii for deep atmosphere --- src/core_atmosphere/Registry.xml | 21 +++++++ src/core_atmosphere/mpas_atm_core.F | 88 +++++++++++++++++++++++++++++ 2 files changed, 109 insertions(+) diff --git a/src/core_atmosphere/Registry.xml b/src/core_atmosphere/Registry.xml index ddcc94d159..e067a1201f 100644 --- a/src/core_atmosphere/Registry.xml +++ b/src/core_atmosphere/Registry.xml @@ -489,6 +489,10 @@ + + + + @@ -899,6 +903,10 @@ + + + + @@ -1437,6 +1445,19 @@ + + + + + + + + + diff --git a/src/core_atmosphere/mpas_atm_core.F b/src/core_atmosphere/mpas_atm_core.F index 997d7ca8ba..577b9dd6fb 100644 --- a/src/core_atmosphere/mpas_atm_core.F +++ b/src/core_atmosphere/mpas_atm_core.F @@ -568,6 +568,8 @@ subroutine atm_mpas_init_block(dminfo, stream_manager, block, mesh, dt) call atm_compute_mesh_scaling(mesh, block % configs) + call atm_compute_nondimensional_radii(mesh, block % configs) + call atm_compute_damping_coefs(mesh, block % configs) ! @@ -1224,6 +1226,92 @@ subroutine atm_compute_signs(mesh) end subroutine atm_compute_signs + subroutine atm_compute_nondimensional_radii(mesh, configs) + + implicit none + + type (mpas_pool_type), intent(inout) :: mesh + type (mpas_pool_type), intent(in) :: configs + + logical, pointer :: config_deep_atmosphere, on_a_sphere + integer :: iCell, iCell1, iCell2, iCell3, i, k, nz + integer, pointer :: nCells, nEdges, nVertices, nz1 + integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex + real (kind=RKIND) :: r_k, r_kp1 + real (kind=RKIND), pointer :: sphere_radius + real (kind=RKIND), dimension(:), pointer :: areaTriangle, dcEdge + real (kind=RKIND), dimension(:,:), pointer :: rTildeCell, rTildeLayer, rTildeEdge, rTildeVertex + real (kind=RKIND), dimension(:,:), pointer :: zxu, zgrid, kiteAreasOnVertex + + call mpas_pool_get_config(configs, 'config_deep_atmosphere', config_deep_atmosphere) + call mpas_pool_get_config(mesh, 'on_a_sphere', on_a_sphere) + call mpas_pool_get_config(mesh, 'sphere_radius', sphere_radius) + + call mpas_pool_get_dimension(mesh, 'nCells', nCells) + call mpas_pool_get_dimension(mesh, 'nEdges', nEdges) + call mpas_pool_get_dimension(mesh, 'nVertices', nVertices) + call mpas_pool_get_dimension(mesh, 'nVertLevels', nz1) + nz = nz1 + 1 + + call mpas_pool_get_array(mesh, 'rTildeCell', rTildeCell) + call mpas_pool_get_array(mesh, 'rTildeLayer', rTildeLayer) + call mpas_pool_get_array(mesh, 'rTildeEdge', rTildeEdge) + call mpas_pool_get_array(mesh, 'rTildeVertex', rTildeVertex) + call mpas_pool_get_array(mesh, 'zxu', zxu) + call mpas_pool_get_array(mesh, 'zgrid', zgrid) + call mpas_pool_get_array(mesh, 'kiteAreasOnVertex', kiteAreasOnVertex) + call mpas_pool_get_array(mesh, 'areaTriangle', areaTriangle) + call mpas_pool_get_array(mesh, 'dcEdge', dcEdge) + call mpas_pool_get_array(mesh, 'cellsOnEdge', cellsOnEdge) + call mpas_pool_get_array(mesh, 'cellsOnVertex', cellsOnVertex) + + if (config_deep_atmosphere.and.on_a_sphere) then + do iCell=1,nCells + do k=1,nz1 + r_k = sphere_radius + zgrid(k,iCell) + r_kp1 = sphere_radius + zgrid(k+1,iCell) + rTildeCell(k,iCell) = sqrt( (r_k**2 + r_k*r_kp1 + r_kp1**2) / (3.*sphere_radius**2) ) + end do + do k=1,nz + rTildeLayer(k,iCell) = 1.0 + (zgrid(k,iCell)/sphere_radius) + end do + end do + do i=1, nEdges + iCell1 = cellsOnEdge(1,i) + iCell2 = cellsOnEdge(2,i) + do k=1,nz1 + rTildeEdge(k,i) = 1.0 + 0.25*(zgrid(k,iCell1)+zgrid(k,iCell2)+zgrid(k+1,iCell1)+zgrid(k+1,iCell2))/sphere_radius + end do + end do + do i=1, nVertices + iCell1 = cellsOnVertex(1,i) + iCell2 = cellsOnVertex(2,i) + iCell3 = cellsOnVertex(3,i) + do k=1,nz1 + rTildeVertex(k,i) = 0.5 * ( (kiteAreasOnVertex(1,i)+kiteAreasOnVertex(2,i)) * rTildeCell(k,iCell3) & + + (kiteAreasOnVertex(2,i)+kiteAreasOnVertex(3,i)) * rTildeCell(k,iCell1) & + + (kiteAreasOnVertex(3,i)+kiteAreasOnVertex(1,i)) * rTildeCell(k,iCell2) ) & + / areaTriangle(i) + end do + end do + else ! shallow or on a plane + rTildeCell(:,:) = 1.0_RKIND + rTildeEdge(:,:) = 1.0_RKIND + rTildeVertex(:,:) = 1.0_RKIND + rTildeLayer(:,:) = 1.0_RKIND + end if + + do i=1, nEdges + iCell1 = cellsOnEdge(1,i) + iCell2 = cellsOnEdge(2,i) + do k=1,nz1 + zxu (k,i) = 0.5 * (zgrid(k,iCell2)-zgrid(k,iCell1) + zgrid(k+1,iCell2)-zgrid(k+1,iCell1)) / (dcEdge(i)*rTildeEdge(k,i)) + end do + end do + + end subroutine atm_compute_nondimensional_radii + + subroutine atm_compute_damping_coefs(mesh, configs) implicit none From 63d41a199457d11e1f8e3ce2abab4fda5736ec58 Mon Sep 17 00:00:00 2001 From: Hing Ong Date: Sat, 23 Nov 2024 08:14:43 -0700 Subject: [PATCH 3/6] Declare the nondimensional radii as needed and use them in mpas_atm_core and mpas_pv_diagnostics --- .../diagnostics/mpas_pv_diagnostics.F | 5 +- .../dynamics/mpas_atm_time_integration.F | 67 ++++++++++++++++++- src/core_atmosphere/mpas_atm_core.F | 5 +- 3 files changed, 71 insertions(+), 6 deletions(-) diff --git a/src/core_atmosphere/diagnostics/mpas_pv_diagnostics.F b/src/core_atmosphere/diagnostics/mpas_pv_diagnostics.F index d21061b0fb..f1c3b98d51 100644 --- a/src/core_atmosphere/diagnostics/mpas_pv_diagnostics.F +++ b/src/core_atmosphere/diagnostics/mpas_pv_diagnostics.F @@ -1268,7 +1268,7 @@ subroutine atm_compute_pv_diagnostics(state, time_lev, diag, mesh) integer :: iCell, k integer, pointer :: nCells, nVertLevels, index_qv real (kind=RKIND) :: pvuVal, missingVal, stratoPV - real (kind=RKIND), dimension(:,:), pointer :: theta, rho, theta_m, rho_zz, zz, dtheta_dt_mix, tend_theta_euler + real (kind=RKIND), dimension(:,:), pointer :: theta, rho, theta_m, rho_zz, zz, dtheta_dt_mix, tend_theta_euler, rTildeCell type (field2DReal), pointer :: theta_f, uReconstructX_f, uReconstructY_f, uReconstructZ_f, w_f, epv_f real (kind=RKIND), dimension(:,:,:), pointer :: scalars @@ -1284,11 +1284,12 @@ subroutine atm_compute_pv_diagnostics(state, time_lev, diag, mesh) call mpas_pool_get_array(diag, 'rho', rho) call mpas_pool_get_array(mesh, 'zz', zz) + call mpas_pool_get_array(mesh, 'rTildeCell', rTildeCell) do iCell=1,nCells do k=1,nVertLevels theta(k,iCell) = theta_m(k,iCell) / (1._RKIND + rvord * scalars(index_qv,k,iCell)) - rho(k,iCell) = rho_zz(k,iCell) * zz(k,iCell) + rho(k,iCell) = rho_zz(k,iCell) * zz(k,iCell) / rTildeCell(k,iCell)**2 ! deep end do end do diff --git a/src/core_atmosphere/dynamics/mpas_atm_time_integration.F b/src/core_atmosphere/dynamics/mpas_atm_time_integration.F index e2bafe8752..5188d598fd 100644 --- a/src/core_atmosphere/dynamics/mpas_atm_time_integration.F +++ b/src/core_atmosphere/dynamics/mpas_atm_time_integration.F @@ -1765,6 +1765,7 @@ subroutine atm_compute_vert_imp_coefs(state, mesh, diag, configs, nVertLevels, d real (kind=RKIND), dimension(:,:), pointer :: zz, cqw, p, t, rb, rtb, pb, rt real (kind=RKIND), dimension(:,:), pointer :: cofwr, cofwz, coftz, cofwt, a_tri, alpha_tri, gamma_tri + real (kind=RKIND), dimension(:,:), pointer :: rTildeCell, rTildeLayer real (kind=RKIND), dimension(:), pointer :: cofrz, rdzw, fzm, fzp, rdzu real (kind=RKIND), dimension(:,:,:), pointer :: scalars @@ -1781,6 +1782,9 @@ subroutine atm_compute_vert_imp_coefs(state, mesh, diag, configs, nVertLevels, d call mpas_pool_get_array(mesh, 'fzp', fzp) call mpas_pool_get_array(mesh, 'zz', zz) + call mpas_pool_get_array(mesh, 'rTildeCell', rTildeCell) + call mpas_pool_get_array(mesh, 'rTildeLayer', rTildeLayer) + call mpas_pool_get_array(diag, 'cqw', cqw) call mpas_pool_get_array(diag, 'exner', p) call mpas_pool_get_array(diag, 'exner_base', pb) @@ -1807,6 +1811,7 @@ subroutine atm_compute_vert_imp_coefs(state, mesh, diag, configs, nVertLevels, d call atm_compute_vert_imp_coefs_work(nCells, moist_start, moist_end, dts, epssm, & zz, cqw, p, t, rb, rtb, pb, rt, cofwr, cofwz, coftz, cofwt, & + rTildeCell, rTildeLayer, & a_tri, alpha_tri, gamma_tri, cofrz, rdzw, fzm, fzp, rdzu, scalars, & cellStart, cellEnd, edgeStart, edgeEnd, & cellSolveStart, cellSolveEnd, edgeSolveStart, edgeSolveEnd) @@ -1817,6 +1822,7 @@ end subroutine atm_compute_vert_imp_coefs subroutine atm_compute_vert_imp_coefs_work(nCells, moist_start, moist_end, dts, epssm, & zz, cqw, p, t, rb, rtb, pb, rt, cofwr, cofwz, coftz, cofwt, & + rTildeCell, rTildeLayer, & a_tri, alpha_tri, gamma_tri, cofrz, rdzw, fzm, fzp, rdzu, scalars, & cellStart, cellEnd, edgeStart, edgeEnd, & cellSolveStart, cellSolveEnd, edgeSolveStart, edgeSolveEnd) @@ -1848,6 +1854,8 @@ subroutine atm_compute_vert_imp_coefs_work(nCells, moist_start, moist_end, dts, real (kind=RKIND), dimension(nVertLevels,nCells+1) :: a_tri real (kind=RKIND), dimension(nVertLevels,nCells+1) :: alpha_tri real (kind=RKIND), dimension(nVertLevels,nCells+1) :: gamma_tri + real (kind=RKIND), dimension(nVertLevels,nCells+1) :: rTildeCell + real (kind=RKIND), dimension(nVertLevels+1,nCells+1) :: rTildeLayer real (kind=RKIND), dimension(nVertLevels) :: cofrz real (kind=RKIND), dimension(nVertLevels) :: rdzw real (kind=RKIND), dimension(nVertLevels) :: fzm @@ -2157,6 +2165,8 @@ subroutine atm_advance_acoustic_step( state, diag, tend, mesh, configs, nCells, real (kind=RKIND), dimension(:,:), pointer :: rw real (kind=RKIND), dimension(:,:), pointer :: rw_save + real (kind=RKIND), dimension(:,:), pointer :: rTildeCell, rTildeEdge + real (kind=RKIND), dimension(:), pointer :: fzm, fzp, rdzw, dcEdge, invDcEdge, invAreaCell, cofrz, dvEdge integer, dimension(:), pointer :: nEdgesOnCell @@ -2227,6 +2237,9 @@ subroutine atm_advance_acoustic_step( state, diag, tend, mesh, configs, nCells, call mpas_pool_get_array(mesh, 'cf2', cf2) call mpas_pool_get_array(mesh, 'cf3', cf3) + call mpas_pool_get_array(mesh, 'rTildeCell', rTildeCell) + call mpas_pool_get_array(mesh, 'rTildeEdge', rTildeEdge) + ! redefine ru_p to be perturbation from time t, change 3b ! temporary call mpas_pool_get_array(diag, 'ru', ru) call mpas_pool_get_array(diag, 'ru_save', ru_save) @@ -2244,6 +2257,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, & + rTildeCell, rTildeEdge, & specZoneMaskEdge, specZoneMaskCell & ) @@ -2257,6 +2271,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, & + rTildeCell, rTildeEdge, & specZoneMaskEdge, specZoneMaskCell & ) @@ -2307,6 +2322,9 @@ subroutine atm_advance_acoustic_step_work(nCells, nEdges, nCellsSolve, cellStart real (kind=RKIND), dimension(nVertLevels+1,nCells+1) :: rw real (kind=RKIND), dimension(nVertLevels+1,nCells+1) :: rw_save + real (kind=RKIND), dimension(nVertLevels,nCells+1) :: rTildeCell + real (kind=RKIND), dimension(nVertLevels,nEdges+1) :: rTildeEdge + real (kind=RKIND), dimension(nVertLevels) :: fzm real (kind=RKIND), dimension(nVertLevels) :: fzp real (kind=RKIND), dimension(nVertLevels) :: rdzw @@ -2543,6 +2561,7 @@ subroutine atm_divergence_damping_3d( state, diag, mesh, configs, dts, edgeStart integer, intent(in) :: edgeStart, edgeEnd real (kind=RKIND), dimension(:,:), pointer :: theta_m, ru_p, rtheta_pp, rtheta_pp_old + real (kind=RKIND), dimension(:,:), pointer :: rTildeCell, rTildeEdge ! real (kind=RKIND), dimension(:), pointer :: dcEdge real (kind=RKIND), pointer :: smdiv, config_len_disp real (kind=RKIND), dimension(:), pointer :: specZoneMaskEdge @@ -2556,6 +2575,8 @@ subroutine atm_divergence_damping_3d( state, diag, mesh, configs, dts, edgeStart call mpas_pool_get_array(mesh, 'cellsOnEdge', cellsOnEdge) call mpas_pool_get_array(mesh, 'specZoneMaskEdge', specZoneMaskEdge) + call mpas_pool_get_array(mesh, 'rTildeCell', rTildeCell) + call mpas_pool_get_array(mesh, 'rTildeEdge', rTildeEdge) ! call mpas_pool_get_array(mesh, 'dcEdge', dcEdge) call mpas_pool_get_array(state, 'theta_m', theta_m, 1) call mpas_pool_get_array(diag, 'rtheta_pp', rtheta_pp) @@ -2628,6 +2649,9 @@ subroutine atm_recover_large_step_variables( state, diag, tend, mesh, configs, d rho_pp, rho_zz, rho_base, ruAvg, ru_save, ru_p, u, ru, & exner, exner_base, rtheta_base, pressure_p, & zz, theta_m, pressure_b + + real (kind=RKIND), dimension(:,:), pointer :: rTildeCell, rTildeEdge + real (kind=RKIND), dimension(:,:,:), pointer :: scalars real (kind=RKIND), dimension(:), pointer :: fzm, fzp real (kind=RKIND), dimension(:,:,:), pointer :: zb, zb3, zb_cell, zb3_cell @@ -2698,6 +2722,8 @@ subroutine atm_recover_large_step_variables( state, diag, tend, mesh, configs, d call mpas_pool_get_array(mesh, 'cf2', cf2) call mpas_pool_get_array(mesh, 'cf3', cf3) + call mpas_pool_get_array(mesh, 'rTildeCell', rTildeCell) + call mpas_pool_get_array(mesh, 'rTildeEdge', rTildeEdge) call atm_recover_large_step_variables_work(nCells, nEdges, nCellsSolve, nEdgesSolve, dt, ns, rk_step, & wwAvg, rw_save, w, rw, rw_p, rtheta_p, rtheta_pp, rtheta_p_save, rt_diabatic_tend, rho_p, & @@ -2705,6 +2731,7 @@ subroutine atm_recover_large_step_variables( state, diag, tend, mesh, configs, d rtheta_base, pressure_p, zz, theta_m, pressure_b, scalars, fzm, fzp, & zb, zb3, zb_cell, zb3_cell, edgesOnCell_sign, cellsOnEdge, edgesOnCell, nEdgesOnCell, & cf1, cf2, cf3, & + rTildeCell, rTildeEdge, & bdyMaskCell, & ! addition for regional_MPAS cellStart, cellEnd, vertexStart, vertexEnd, edgeStart, edgeEnd, & cellSolveStart, cellSolveEnd, vertexSolveStart, vertexSolveEnd, edgeSolveStart, edgeSolveEnd) @@ -2718,6 +2745,7 @@ subroutine atm_recover_large_step_variables_work(nCells, nEdges, nCellsSolve, nE rtheta_base, pressure_p, zz, theta_m, pressure_b, scalars, fzm, fzp, & zb, zb3, zb_cell, zb3_cell, edgesOnCell_sign, cellsOnEdge, edgesOnCell, nEdgesOnCell, & cf1, cf2, cf3, & + rTildeCell, rTildeEdge, & bdyMaskCell, & ! addition for regional_MPAS cellStart, cellEnd, vertexStart, vertexEnd, edgeStart, edgeEnd, & cellSolveStart, cellSolveEnd, vertexSolveStart, vertexSolveEnd, edgeSolveStart, edgeSolveEnd) @@ -2762,6 +2790,8 @@ subroutine atm_recover_large_step_variables_work(nCells, nEdges, nCellsSolve, nE real (kind=RKIND), dimension(nVertLevels,nCells+1) :: zz real (kind=RKIND), dimension(nVertLevels,nCells+1) :: theta_m real (kind=RKIND), dimension(nVertLevels,nCells+1) :: pressure_b + real (kind=RKIND), dimension(nVertLevels,nCells+1) :: rTildeCell + real (kind=RKIND), dimension(nVertLevels,nEdges+1) :: rTildeEdge real (kind=RKIND), dimension(num_scalars,nVertLevels,nCells+1) :: scalars real (kind=RKIND), dimension(nVertLevels) :: fzm real (kind=RKIND), dimension(nVertLevels) :: fzp @@ -4290,6 +4320,8 @@ subroutine atm_compute_dyn_tend(tend, tend_physics, state, diag, mesh, configs, real(kind=RKIND), dimension(:,:), pointer :: tend_w_pgf, tend_w_buoy + real(kind=RKIND), dimension(:,:), pointer :: rTildeCell, rTildeLayer, rTildeEdge + real (kind=RKIND), pointer :: coef_3rd_order, c_s logical, pointer :: config_mix_full character (len=StrKIND), pointer :: config_horiz_mixing @@ -4303,6 +4335,8 @@ subroutine atm_compute_dyn_tend(tend, tend_physics, state, diag, mesh, configs, real (kind=RKIND), pointer :: config_mpas_cam_coef logical, pointer :: config_rayleigh_damp_u + logical, pointer :: config_deep_atmosphere + logical, pointer :: on_a_sphere real (kind=RKIND), pointer :: config_rayleigh_damp_u_timescale_days integer, pointer :: config_number_rayleigh_damp_u_levels, config_number_cam_damping_levels @@ -4326,6 +4360,8 @@ 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) + call mpas_pool_get_config(configs, 'config_deep_atmosphere', config_deep_atmosphere) + call mpas_pool_get_config(mesh, 'on_a_sphere', on_a_sphere) call mpas_pool_get_array(state, 'rho_zz', rho_zz, 2) call mpas_pool_get_array(state, 'u', u, 2) @@ -4432,6 +4468,10 @@ subroutine atm_compute_dyn_tend(tend, tend_physics, state, diag, mesh, configs, call mpas_pool_get_array(mesh, 'cf2', cf2) call mpas_pool_get_array(mesh, 'cf3', cf3) + call mpas_pool_get_array(mesh, 'rTildeCell', rTildeCell) + call mpas_pool_get_array(mesh, 'rTildeLayer', rTildeLayer) + call mpas_pool_get_array(mesh, 'rTildeEdge', rTildeEdge) + call atm_compute_dyn_tend_work(nCells, nEdges, nVertices, nVertLevels, & nCellsSolve, nEdgesSolve, vertexDegree, maxEdges, maxEdges2, num_scalars, moist_start, moist_end, & fEdge, dvEdge, dcEdge, invDcEdge, invDvEdge, invAreaCell, invAreaTriangle, meshScalingDel2, meshScalingDel4, & @@ -4443,12 +4483,14 @@ subroutine atm_compute_dyn_tend(tend, tend_physics, state, diag, mesh, configs, cellsOnEdge, verticesOnEdge, edgesOnCell, edgesOnEdge, cellsOnCell, edgesOnVertex, nEdgesOnCell, nEdgesOnEdge, & latCell, latEdge, angleEdge, u_init, v_init, advCellsForEdge, nAdvCellsForEdge, adv_coefs, adv_coefs_3rd, & rdzu, rdzw, fzm, fzp, qv_init, t_init, cf1, cf2, cf3, r_earth, ur_cell, vr_cell, defc_a, defc_b, & + rTildeCell, rTildeLayer, rTildeEdge, & tend_w_pgf, tend_w_buoy, coef_3rd_order, c_s, config_mix_full, config_horiz_mixing, config_del4u_div_factor, & config_h_mom_eddy_visc2, config_v_mom_eddy_visc2, config_h_theta_eddy_visc2, config_v_theta_eddy_visc2, & config_h_theta_eddy_visc4, config_h_mom_eddy_visc4, config_visc4_2dsmag, config_len_disp, rk_step, dt, & 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, & + config_deep_atmosphere, on_a_sphere, & rthdynten, & cellStart, cellEnd, vertexStart, vertexEnd, edgeStart, edgeEnd, & cellSolveStart, cellSolveEnd, vertexSolveStart, vertexSolveEnd, edgeSolveStart, edgeSolveEnd) @@ -4467,12 +4509,14 @@ subroutine atm_compute_dyn_tend_work(nCells, nEdges, nVertices, nVertLevels_dumm cellsOnEdge, verticesOnEdge, edgesOnCell, edgesOnEdge, cellsOnCell, edgesOnVertex, nEdgesOnCell, nEdgesOnEdge, & latCell, latEdge, angleEdge, u_init, v_init, advCellsForEdge, nAdvCellsForEdge, adv_coefs, adv_coefs_3rd, & rdzu, rdzw, fzm, fzp, qv_init, t_init, cf1, cf2, cf3, r_earth, ur_cell, vr_cell, defc_a, defc_b, & + rTildeCell, rTildeLayer, rTildeEdge, & tend_w_pgf, tend_w_buoy, coef_3rd_order, c_s, config_mix_full, config_horiz_mixing, config_del4u_div_factor, & config_h_mom_eddy_visc2, config_v_mom_eddy_visc2, config_h_theta_eddy_visc2, config_v_theta_eddy_visc2, & config_h_theta_eddy_visc4, config_h_mom_eddy_visc4, config_visc4_2dsmag, config_len_disp, rk_step, dt, & 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, & + config_deep_atmosphere, on_a_sphere, & rthdynten, & cellStart, cellEnd, vertexStart, vertexEnd, edgeStart, edgeEnd, & cellSolveStart, cellSolveEnd, vertexSolveStart, vertexSolveEnd, edgeSolveStart, edgeSolveEnd) @@ -4579,6 +4623,10 @@ subroutine atm_compute_dyn_tend_work(nCells, nEdges, nVertices, nVertLevels_dumm real (kind=RKIND), dimension(nVertLevels+1,nCells+1) :: tend_w_pgf real (kind=RKIND), dimension(nVertLevels+1,nCells+1) :: tend_w_buoy + real (kind=RKIND), dimension(nVertLevels,nCells+1) :: rTildeCell + real (kind=RKIND), dimension(nVertLevels+1,nCells+1) :: rTildeLayer + real (kind=RKIND), dimension(nVertLevels,nEdges+1) :: rTildeEdge + real (kind=RKIND) :: coef_3rd_order, c_s logical :: config_mix_full character (len=StrKIND) :: config_horiz_mixing @@ -4595,6 +4643,8 @@ subroutine atm_compute_dyn_tend_work(nCells, nEdges, nVertices, nVertLevels_dumm real (kind=RKIND) :: config_mpas_cam_coef logical, intent(in) :: config_rayleigh_damp_u + logical, intent(in) :: config_deep_atmosphere + logical, intent(in) :: on_a_sphere real (kind=RKIND), intent(in) :: config_rayleigh_damp_u_timescale_days integer, intent(in) :: config_number_rayleigh_damp_u_levels, config_number_cam_damping_levels @@ -4615,7 +4665,7 @@ subroutine atm_compute_dyn_tend_work(nCells, nEdges, nVertices, nVertLevels_dumm real (kind=RKIND), dimension( nVertLevels ) :: ru_edge_w, q, u_mix real (kind=RKIND) :: theta_turb_flux, w_turb_flux, r real (kind=RKIND) :: scalar_weight - real (kind=RKIND) :: inv_r_earth + real (kind=RKIND) :: inv_r_earth, r_c real (kind=RKIND) :: invDt, flux, workpv real (kind=RKIND) :: edge_sign, pr_scale, r_dc, r_dv, u_mix_scale @@ -5423,6 +5473,7 @@ subroutine atm_compute_solve_diagnostics(dt, state, time_lev, diag, mesh, config real (kind=RKIND), dimension(:,:), pointer :: weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, & vorticity, ke, pv_edge, pv_vertex, pv_cell, gradPVn, gradPVt, & divergence + real (kind=RKIND), dimension(:,:), pointer :: rTildeCell, rTildeEdge, rTildeVertex integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex, & kiteForCell, verticesOnCell real (kind=RKIND), dimension(:,:), pointer :: edgesOnVertex_sign, edgesOnCell_sign @@ -5470,6 +5521,10 @@ subroutine atm_compute_solve_diagnostics(dt, state, time_lev, diag, mesh, config call mpas_pool_get_array(mesh, 'fVertex', fVertex) call mpas_pool_get_array(mesh, 'fEdge', fEdge) + call mpas_pool_get_array(mesh, 'rTildeCell', rTildeCell) + call mpas_pool_get_array(mesh, 'rTildeEdge', rTildeEdge) + call mpas_pool_get_array(mesh, 'rTildeVertex', rTildeVertex) + call mpas_pool_get_dimension(mesh, 'nCells', nCells) call mpas_pool_get_dimension(mesh, 'nEdges', nEdges) call mpas_pool_get_dimension(mesh, 'nVertices', nVertices) @@ -5479,6 +5534,7 @@ subroutine atm_compute_solve_diagnostics(dt, state, time_lev, diag, mesh, config call atm_compute_solve_diagnostics_work(nCells, nEdges, nVertices, & vertexDegree, dt, config_apvm_upwinding, & fVertex, fEdge, invAreaTriangle, invAreaCell, dvEdge, dcEdge, invDvEdge, invDcEdge, & + rTildeCell, rTildeEdge, rTildeVertex, & weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, vorticity, ke, pv_edge, pv_vertex, pv_cell, & gradPVn, gradPVt, divergence, cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, edgesOnEdge, & edgesOnVertex, kiteForCell, verticesOnCell, edgesOnVertex_sign, edgesOnCell_sign, nEdgesOnCell, nEdgesOnEdge, & @@ -5491,6 +5547,7 @@ end subroutine atm_compute_solve_diagnostics subroutine atm_compute_solve_diagnostics_work(nCells, nEdges, nVertices, & vertexDegree, dt, config_apvm_upwinding, & fVertex, fEdge, invAreaTriangle, invAreaCell, dvEdge, dcEdge, invDvEdge, invDcEdge, & + rTildeCell, rTildeEdge, rTildeVertex, & weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, vorticity, ke, pv_edge, pv_vertex, pv_cell, & gradPVn, gradPVt, divergence, cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, edgesOnEdge, & edgesOnVertex, kiteForCell, verticesOnCell, edgesOnVertex_sign, edgesOnCell_sign, nEdgesOnCell, nEdgesOnEdge, & @@ -5528,6 +5585,9 @@ subroutine atm_compute_solve_diagnostics_work(nCells, nEdges, nVertices, & real (kind=RKIND), dimension(nVertLevels,nEdges+1) :: gradPVn real (kind=RKIND), dimension(nVertLevels,nEdges+1) :: gradPVt real (kind=RKIND), dimension(nVertLevels,nCells+1) :: divergence + real (kind=RKIND), dimension(nVertLevels,nCells+1) :: rTildeCell + real (kind=RKIND), dimension(nVertLevels,nEdges+1) :: rTildeEdge + real (kind=RKIND), dimension(nVertLevels,nVertices+1) :: rTildeVertex integer, dimension(2,nEdges+1) :: cellsOnEdge integer, dimension(3,nVertices+1) :: cellsOnVertex integer, dimension(2,nEdges+1) :: verticesOnEdge @@ -5840,6 +5900,8 @@ subroutine atm_init_coupled_diagnostics(state, time_lev, diag, mesh, configs, & real (kind=RKIND), dimension(:,:), pointer :: pressure_base real (kind=RKIND), dimension(:,:), pointer :: exner real (kind=RKIND), dimension(:,:), pointer :: exner_base + real (kind=RKIND), dimension(:,:), pointer :: rTildeCell + real (kind=RKIND), dimension(:,:), pointer :: rTildeEdge real (kind=RKIND), dimension(:), pointer :: fzm, fzp real (kind=RKIND), dimension(:,:,:), pointer :: zb, zb3, zb_cell, zb3_cell @@ -5879,7 +5941,8 @@ subroutine atm_init_coupled_diagnostics(state, time_lev, diag, mesh, configs, & call mpas_pool_get_array(mesh, 'zb3', zb3) call mpas_pool_get_array(mesh, 'zb_cell', zb_cell) call mpas_pool_get_array(mesh, 'zb3_cell', zb3_cell) - + call mpas_pool_get_array(mesh, 'rTildeCell', rTildeCell) + call mpas_pool_get_array(mesh, 'rTildeEdge', rTildeEdge) rcv = rgas / (cp-rgas) p0 = 1.e5 ! this should come from somewhere else... diff --git a/src/core_atmosphere/mpas_atm_core.F b/src/core_atmosphere/mpas_atm_core.F index 577b9dd6fb..ebf6eebdee 100644 --- a/src/core_atmosphere/mpas_atm_core.F +++ b/src/core_atmosphere/mpas_atm_core.F @@ -907,7 +907,7 @@ subroutine atm_compute_output_diagnostics(state, time_lev, diag, mesh) integer :: iCell, k integer, pointer :: nCells, nVertLevels, index_qv - real (kind=RKIND), dimension(:,:), pointer :: theta, rho, theta_m, rho_zz, zz + real (kind=RKIND), dimension(:,:), pointer :: theta, rho, theta_m, rho_zz, zz, rTildeCell real (kind=RKIND), dimension(:,:), pointer :: pressure_base, pressure_p, pressure real (kind=RKIND), dimension(:,:,:), pointer :: scalars @@ -926,11 +926,12 @@ subroutine atm_compute_output_diagnostics(state, time_lev, diag, mesh) call mpas_pool_get_array(diag, 'pressure', pressure) call mpas_pool_get_array(mesh, 'zz', zz) + call mpas_pool_get_array(mesh, 'rTildeCell', rTildeCell) do iCell=1,nCells do k=1,nVertLevels theta(k,iCell) = theta_m(k,iCell) / (1._RKIND + rvord * scalars(index_qv,k,iCell)) - rho(k,iCell) = rho_zz(k,iCell) * zz(k,iCell) + rho(k,iCell) = rho_zz(k,iCell) * zz(k,iCell) / rTildeCell(k,iCell)**2 ! deep pressure(k,iCell) = pressure_base(k,iCell) + pressure_p(k,iCell) end do end do From 8d4bfcd23bb7631b5a35cb137ecde3937dae7876 Mon Sep 17 00:00:00 2001 From: Hing Ong Date: Sat, 23 Nov 2024 09:39:20 -0700 Subject: [PATCH 4/6] Move call atm_compute_nondimensional_radii before initialization --- src/core_atmosphere/mpas_atm_core.F | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/core_atmosphere/mpas_atm_core.F b/src/core_atmosphere/mpas_atm_core.F index ebf6eebdee..65484de0d8 100644 --- a/src/core_atmosphere/mpas_atm_core.F +++ b/src/core_atmosphere/mpas_atm_core.F @@ -469,7 +469,9 @@ subroutine atm_mpas_init_block(dminfo, stream_manager, block, mesh, dt) call atm_adv_coef_compression(mesh) call atm_couple_coef_3rd_order(mesh, block % configs) - + + call atm_compute_nondimensional_radii(mesh, block % configs) + call mpas_pool_get_dimension(state, 'nVertices', nVertices) call mpas_pool_get_dimension(state, 'nVertLevels', nVertLevels) @@ -568,8 +570,6 @@ subroutine atm_mpas_init_block(dminfo, stream_manager, block, mesh, dt) call atm_compute_mesh_scaling(mesh, block % configs) - call atm_compute_nondimensional_radii(mesh, block % configs) - call atm_compute_damping_coefs(mesh, block % configs) ! From 0dac860bb1a5a74f7e531ee6926e2f4b73df0587 Mon Sep 17 00:00:00 2001 From: Hing Ong Date: Sat, 23 Nov 2024 10:38:01 -0700 Subject: [PATCH 5/6] Implement the deep atmosphere equations --- .../dynamics/mpas_atm_time_integration.F | 130 ++++++++++-------- 1 file changed, 70 insertions(+), 60 deletions(-) diff --git a/src/core_atmosphere/dynamics/mpas_atm_time_integration.F b/src/core_atmosphere/dynamics/mpas_atm_time_integration.F index 5188d598fd..1fa5ff7755 100644 --- a/src/core_atmosphere/dynamics/mpas_atm_time_integration.F +++ b/src/core_atmosphere/dynamics/mpas_atm_time_integration.F @@ -1895,7 +1895,7 @@ subroutine atm_compute_vert_imp_coefs_work(nCells, moist_start, moist_end, dts, !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)) + *rdzu(k)*cqw(k,iCell)*(fzm(k)*p (k,iCell)+fzp(k)*p (k-1,iCell))*rTildeLayer(k,iCell)**2 coftz(k,iCell) = dtseps* (fzm(k)*t (k,iCell)+fzp(k)*t (k-1,iCell)) end do coftz(nVertLevels+1,iCell) = 0.0 @@ -1908,7 +1908,7 @@ subroutine atm_compute_vert_imp_coefs_work(nCells, moist_start, moist_end, dts, ! end do qtotal = qtot(k,iCell) - cofwt(k,iCell) = .5*dtseps*rcv*zz(k,iCell)*gravity*rb(k,iCell)/(1.+qtotal) & + cofwt(k,iCell) = .5*dtseps*rcv*zz(k,iCell)*gravity/rTildeCell(k,iCell)**2*rb(k,iCell)/(1.+qtotal) & *p(k,iCell)/((rtb(k,iCell)+rt(k,iCell))*pb(k,iCell)) ! cofwt(k,iCell) = 0. end do @@ -2383,9 +2383,11 @@ subroutine atm_advance_acoustic_step_work(nCells, nEdges, nCellsSolve, cellStart !DIR$ IVDEP do k=1,nVertLevels - pgrad = ((rtheta_pp(k,cell2)-rtheta_pp(k,cell1))*invDcEdge(iEdge) )/(.5*(zz(k,cell2)+zz(k,cell1))) + pgrad = ((rtheta_pp(k,cell2)/rTildeCell(k,cell2)**2-rtheta_pp(k,cell1)/rTildeCell(k,cell1)**2)*invDcEdge(iEdge))& + /(.5*(zz(k,cell2)+zz(k,cell1))) pgrad = cqu(k,iEdge)*0.5*c2*(exner(k,cell1)+exner(k,cell2))*pgrad - pgrad = pgrad + 0.5*zxu(k,iEdge)*gravity*(rho_pp(k,cell1)+rho_pp(k,cell2)) + pgrad = pgrad + 0.5*zxu(k,iEdge)*gravity/rTildeEdge(k,iEdge)**2 & + *(rho_pp(k,cell1)/rTildeCell(k,cell1)**2+rho_pp(k,cell2)/rTildeCell(k,cell2)**2) ru_p(k,iEdge) = ru_p(k,iEdge) + dts*(tend_ru(k,iEdge) - (1.0_RKIND - specZoneMaskEdge(iEdge))*pgrad) end do @@ -2483,12 +2485,10 @@ subroutine atm_advance_acoustic_step_work(nCells, nEdges, nCellsSolve, cellStart !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))) & + - cofwz(k,iCell)*( zz(k ,iCell)*(ts(k )+resm*rtheta_pp(k, iCell))/rTildeCell(k ,iCell)**2 & + -zz(k-1,iCell)*(ts(k-1)+resm*rtheta_pp(k-1,iCell))/rTildeCell(k-1,iCell)**2 ) & + - cofwr(k,iCell)*((rs(k)/rTildeCell(k ,iCell)**2+rs(k-1)/rTildeCell(k-1,iCell)**2) & + +resm*(rho_pp(k,iCell)/rTildeCell(k ,iCell)**2+rho_pp(k-1,iCell)/rTildeCell(k-1,iCell)**2)) & + cofwt(k ,iCell)*(ts(k )+resm*rtheta_pp(k ,iCell)) & + cofwt(k-1,iCell)*(ts(k-1)+resm*rtheta_pp(k-1,iCell)) end do @@ -2610,10 +2610,10 @@ subroutine atm_divergence_damping_3d( state, diag, mesh, configs, dts, edgeStart !! /(theta_m(k,cell1)+theta_m(k,cell2)) !! scaled 3d divergence damping - divCell1 = -(rtheta_pp(k,cell1)-rtheta_pp_old(k,cell1)) - divCell2 = -(rtheta_pp(k,cell2)-rtheta_pp_old(k,cell2)) + divCell1 = -(rtheta_pp(k,cell1)-rtheta_pp_old(k,cell1))/rTildeCell(k,cell1)**2 + divCell2 = -(rtheta_pp(k,cell2)-rtheta_pp_old(k,cell2))/rTildeCell(k,cell2)**2 ru_p(k,iEdge) = ru_p(k,iEdge) + coef_divdamp*(divCell2-divCell1)*(1.0_RKIND - specZoneMaskEdge(iEdge)) & - /(theta_m(k,cell1)+theta_m(k,cell2)) + /(theta_m(k,cell1)+theta_m(k,cell2)) * rTildeEdge(k,iEdge) end do end if ! edges for block-owned cells @@ -2860,9 +2860,9 @@ subroutine atm_recover_large_step_variables_work(nCells, nEdges, nCellsSolve, nE rtheta_p(k,iCell) = rtheta_p_save(k,iCell) + rtheta_pp(k,iCell) & - dt * rho_zz(k,iCell) * rt_diabatic_tend(k,iCell) theta_m(k,iCell) = (rtheta_p(k,iCell) + rtheta_base(k,iCell))/rho_zz(k,iCell) - exner(k,iCell) = (zz(k,iCell)*(rgas/p0)*(rtheta_p(k,iCell)+rtheta_base(k,iCell)))**rcv + exner(k,iCell) = (zz(k,iCell)/rTildeCell(k,iCell)**2*(rgas/p0)*(rtheta_p(k,iCell)+rtheta_base(k,iCell)))**rcv ! pressure_p is perturbation pressure - pressure_p(k,iCell) = zz(k,iCell) * rgas * (exner(k,iCell)*rtheta_p(k,iCell)+rtheta_base(k,iCell) & + pressure_p(k,iCell) = zz(k,iCell) / rTildeCell(k,iCell)**2 * rgas * (exner(k,iCell)*rtheta_p(k,iCell)+rtheta_base(k,iCell) & * (exner(k,iCell)-exner_base(k,iCell))) end do else @@ -2890,7 +2890,8 @@ subroutine atm_recover_large_step_variables_work(nCells, nEdges, nCellsSolve, nE do k = 1, nVertLevels ruAvg(k,iEdge) = ru_save(k,iEdge) + (ruAvg(k,iEdge) * invNs) ru(k,iEdge) = ru_save(k,iEdge) + ru_p(k,iEdge) - u(k,iEdge) = 2.*ru(k,iEdge)/(rho_zz(k,cell1)+rho_zz(k,cell2)) + u(k,iEdge) = 2.*ru(k,iEdge)/(rho_zz(k,cell1)/rTildeCell(k,cell1)**2+rho_zz(k,cell2)/rTildeCell(k,cell2)**2) & + /rTildeEdge(k,iEdge) end do end do @@ -4611,7 +4612,7 @@ subroutine atm_compute_dyn_tend_work(nCells, nEdges, nVertices, nVertLevels_dumm real (kind=RKIND), dimension(nVertLevels,nCells+1) :: t_init real (kind=RKIND) :: cf1, cf2, cf3 - real (kind=RKIND) :: prandtl_inv, r_areaCell, rgas_cprcv + real (kind=RKIND) :: prandtl_inv, r_areaCell!, rgas_cprcv ! not used real (kind=RKIND) :: r_earth real (kind=RKIND), dimension(nVertLevels,nCells+1) :: ur_cell @@ -4689,7 +4690,11 @@ subroutine atm_compute_dyn_tend_work(nCells, nEdges, nVertices, nVertLevels_dumm prandtl_inv = 1.0_RKIND / prandtl invDt = 1.0_RKIND / dt - inv_r_earth = 1.0_RKIND / r_earth + if (on_a_sphere) then + inv_r_earth = 1.0_RKIND / r_earth + else ! on a plane, r_earth -> infinity + inv_r_earth = 0.0_RKIND + endif v_mom_eddy_visc2 = config_v_mom_eddy_visc2 v_theta_eddy_visc2 = config_v_theta_eddy_visc2 @@ -4783,13 +4788,13 @@ subroutine atm_compute_dyn_tend_work(nCells, nEdges, nVertices, nVertLevels_dumm ! if(rk_step == 1) then - rgas_cprcv = rgas*cp/cv + !rgas_cprcv = rgas*cp/cv ! not used do iCell = cellStart,cellEnd !DIR$ IVDEP 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))) + dpdz(k,iCell) = -gravity/rTildeCell(k,iCell)**2*(rb(k,iCell)*(qtot(k,iCell)) + rr_save(k,iCell)*(1.+qtot(k,iCell))) end do end do end if @@ -4811,7 +4816,7 @@ subroutine atm_compute_dyn_tend_work(nCells, nEdges, nVertices, nVertLevels_dumm !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)) ) + -0.5*zxu(k,iEdge)*(dpdz(k,cell1)+dpdz(k,cell2))/rTildeEdge(k,iEdge)**2 ) end do end if @@ -4832,7 +4837,7 @@ subroutine atm_compute_dyn_tend_work(nCells, nEdges, nVertices, nVertLevels_dumm !DIR$ IVDEP do k=1,nVertLevels - tend_u(k,iEdge) = - rdzw(k)*(wduz(k+1)-wduz(k)) ! first use of tend_u + tend_u(k,iEdge) = - rdzw(k)*(wduz(k+1)-wduz(k))/rTildeEdge(k,iEdge) ! first use of tend_u end do ! Next, nonlinear Coriolis term (q) following Ringler et al JCP 2009 @@ -4853,19 +4858,23 @@ subroutine atm_compute_dyn_tend_work(nCells, nEdges, nVertices, nVertLevels_dumm ! horizontal ke gradient and vorticity terms in the vector invariant formulation ! of the horizontal momentum equation - tend_u(k,iEdge) = tend_u(k,iEdge) + rho_edge(k,iEdge)* (q(k) - (ke(k,cell2) - ke(k,cell1)) & - * invDcEdge(iEdge)) & - - u(k,iEdge)*0.5*(h_divergence(k,cell1)+h_divergence(k,cell2)) -#ifdef CURVATURE - ! curvature terms for the sphere - tend_u(k,iEdge) = tend_u(k,iEdge) & - - 2.*omega*cos(angleEdge(iEdge))*cos(latEdge(iEdge)) & - *rho_edge(k,iEdge)*.25*(w(k,cell1)+w(k+1,cell1)+w(k,cell2)+w(k+1,cell2)) & - - u(k,iEdge)*.25*(w(k+1,cell1)+w(k,cell1)+w(k,cell2)+w(k+1,cell2)) & - *rho_edge(k,iEdge) * inv_r_earth -#endif + tend_u(k,iEdge) = tend_u(k,iEdge) + 0.5 * (rho_zz(k,cell1)/rTildeCell(k,cell1)**2 & + +rho_zz(k,cell2)/rTildeCell(k,cell2)**2) & + * (q(k) - (ke(k,cell2)-ke(k,cell1))*invDcEdge(iEdge)) & + - u(k,iEdge)*0.5*(h_divergence(k,cell1)+h_divergence(k,cell2))/rTildeEdge(k,iEdge) end do + if (config_deep_atmosphere) then + do k=1,nVertLevels + tend_u(k,iEdge) = tend_u(k,iEdge) + ( & + - 2.*omega*cos(angleEdge(iEdge))*cos(latEdge(iEdge)) & + *.25*(w(k,cell1)+w(k+1,cell1)+w(k,cell2)+w(k+1,cell2)) & + - u(k,iEdge)*.25*(w(k+1,cell1)+w(k,cell1)+w(k,cell2)+w(k+1,cell2)) & + * inv_r_earth / rTildeEdge(k,iEdge) ) & + * rho_edge(k,iEdge) / rTildeEdge(k,iEdge) + end do + endif + end do @@ -5102,22 +5111,6 @@ subroutine atm_compute_dyn_tend_work(nCells, nEdges, nVertices, nVertLevels_dumm end do end do -#ifdef CURVATURE - do iCell = cellSolveStart, cellSolveEnd -!DIR$ IVDEP - do k=2,nVertLevels - tend_w(k,iCell) = tend_w(k,iCell) + (rho_zz(k,iCell)*fzm(k)+rho_zz(k-1,iCell)*fzp(k))* & - ( (fzm(k)*ur_cell(k,iCell)+fzp(k)*ur_cell(k-1,iCell))**2. & - +(fzm(k)*vr_cell(k,iCell)+fzp(k)*vr_cell(k-1,iCell))**2. )/r_earth & - + 2.*omega*cos(latCell(iCell)) & - *(fzm(k)*ur_cell(k,iCell)+fzp(k)*ur_cell(k-1,iCell)) & - *(rho_zz(k,iCell)*fzm(k)+rho_zz(k-1,iCell)*fzp(k)) - - end do - end do -#endif - - ! ! horizontal mixing for w - we could combine this with advection directly (i.e. as a turbulent flux), ! but here we can also code in hyperdiffusion if we wish (2nd order at present) @@ -5218,13 +5211,29 @@ subroutine atm_compute_dyn_tend_work(nCells, nEdges, nVertices, nVertLevels_dumm !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)) & + rdzu(k)*(pp(k,iCell)-pp(k-1,iCell))*rTildeLayer(k,iCell)**2 & - (fzm(k)*dpdz(k,iCell) + fzp(k)*dpdz(k-1,iCell)) ) ! dpdz is the buoyancy term here. end do end if end do + if (config_deep_atmosphere) then + do iCell = cellSolveStart, cellSolveEnd +!DIR$ IVDEP + do k=2,nVertLevels + r_c = rTildeLayer(k,iCell) + tend_w(k,iCell) = tend_w(k,iCell) + (r_c**2)*( & + fzm(k)*( (ur_cell(k ,iCell)**2 + vr_cell(k ,iCell)**2)*inv_r_earth/r_c & + + 2.*omega*cos(latCell(iCell))*ur_cell(k ,iCell)) & + *rho_zz(k ,iCell)/rTildeCell(k ,iCell)**2 & + +fzp(k)*( (ur_cell(k-1,iCell)**2 + vr_cell(k-1,iCell)**2)*inv_r_earth/r_c & + + 2.*omega*cos(latCell(iCell))*ur_cell(k-1,iCell)) & + *rho_zz(k-1,iCell)/rTildeCell(k-1,iCell)**2 ) + end do + end do + endif + if (rk_step == 1) then if ( v_mom_eddy_visc2 > 0.0 ) then @@ -5648,12 +5657,12 @@ subroutine atm_compute_solve_diagnostics_work(nCells, nEdges, nVertices, & s = edgesOnVertex_sign(i,iVertex) * dcEdge(iEdge) !DIR$ IVDEP do k=1,nVertLevels - vorticity(k,iVertex) = vorticity(k,iVertex) + s * u(k,iEdge) + vorticity(k,iVertex) = vorticity(k,iVertex) + s * u(k,iEdge) * rTildeEdge(k,iEdge) end do end do !DIR$ IVDEP do k=1,nVertLevels - vorticity(k,iVertex) = vorticity(k,iVertex) * invAreaTriangle(iVertex) + vorticity(k,iVertex) = vorticity(k,iVertex) * invAreaTriangle(iVertex) / rTildeVertex(k,iVertex)**2 end do end do @@ -5668,12 +5677,12 @@ subroutine atm_compute_solve_diagnostics_work(nCells, nEdges, nVertices, & s = edgesOnCell_sign(i,iCell) * dvEdge(iEdge) !DIR$ IVDEP do k=1,nVertLevels - divergence(k,iCell) = divergence(k,iCell) + s * u(k,iEdge) + divergence(k,iCell) = divergence(k,iCell) + s * u(k,iEdge) * rTildeEdge(k,iEdge) end do end do r = invAreaCell(iCell) do k = 1,nVertLevels - divergence(k,iCell) = divergence(k,iCell) * r + divergence(k,iCell) = divergence(k,iCell) * r / rTildeCell(k,iCell)**2 end do end do @@ -5950,7 +5959,7 @@ subroutine atm_init_coupled_diagnostics(state, time_lev, diag, mesh, configs, & do iCell=cellStart,cellEnd do k=1,nVertLevels theta_m(k,iCell) = theta(k,iCell) * (1._RKIND + rvord * scalars(index_qv,k,iCell)) - rho_zz(k,iCell) = rho(k,iCell) / zz(k,iCell) + rho_zz(k,iCell) = rho(k,iCell) / zz(k,iCell) * rTildeCell(k,iCell)**2 end do end do @@ -5960,7 +5969,8 @@ subroutine atm_init_coupled_diagnostics(state, time_lev, diag, mesh, configs, & cell1 = cellsOnEdge(1,iEdge) cell2 = cellsOnEdge(2,iEdge) do k=1,nVertLevels - ru(k,iEdge) = 0.5 * u(k,iEdge) * (rho_zz(k,cell1) + rho_zz(k,cell2)) + ru(k,iEdge) = 0.5 * u(k,iEdge) * rTildeEdge(k,iEdge) * (rho_zz(k,cell1) / rTildeCell(k,cell1)**2 & + + rho_zz(k,cell2) / rTildeCell(k,cell2)**2) end do end do @@ -6013,18 +6023,18 @@ subroutine atm_init_coupled_diagnostics(state, time_lev, diag, mesh, configs, & do iCell=cellStart,cellEnd do k=1,nVertLevels - exner(k,iCell) = (zz(k,iCell) * (rgas/p0) * (rtheta_p(k,iCell) + rtheta_base(k,iCell)))**rcv - exner_base(k,iCell) = (zz(k,iCell) * (rgas/p0) * (rtheta_base(k,iCell)))**rcv ! WCS addition 20180403 + exner(k,iCell) = (zz(k,iCell) / rTildeCell(k,iCell)**2 * (rgas/p0) * (rtheta_p(k,iCell) + rtheta_base(k,iCell)))**rcv + exner_base(k,iCell) = (zz(k,iCell) / rTildeCell(k,iCell)**2 * (rgas/p0) * (rtheta_base(k,iCell)))**rcv ! WCS addition 20180403 end do end do do iCell=cellStart,cellEnd do k=1,nVertLevels - pressure_p(k,iCell) = zz(k,iCell) * rgas & + pressure_p(k,iCell) = zz(k,iCell) / rTildeCell(k,iCell)**2 * rgas & * ( exner(k,iCell) * rtheta_p(k,iCell) & + rtheta_base(k,iCell) * (exner(k,iCell) - exner_base(k,iCell)) & ) - pressure_base(k,iCell) = zz(k,iCell) * rgas * exner_base(k,iCell) * rtheta_base(k,iCell) ! WCS addition 20180403 + pressure_base(k,iCell) = zz(k,iCell) / rTildeCell(k,iCell)**2 * rgas * exner_base(k,iCell) * rtheta_base(k,iCell) ! WCS addition 20180403 end do end do From 1847051ad1606b711b620c4177ac4b91ae1ab563 Mon Sep 17 00:00:00 2001 From: Hing Ong Date: Sat, 23 Nov 2024 12:09:06 -0700 Subject: [PATCH 6/6] Correct physics momentum tendency with rTildeEdge --- src/core_atmosphere/dynamics/mpas_atm_time_integration.F | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/core_atmosphere/dynamics/mpas_atm_time_integration.F b/src/core_atmosphere/dynamics/mpas_atm_time_integration.F index 1fa5ff7755..339105df62 100644 --- a/src/core_atmosphere/dynamics/mpas_atm_time_integration.F +++ b/src/core_atmosphere/dynamics/mpas_atm_time_integration.F @@ -5069,7 +5069,7 @@ subroutine atm_compute_dyn_tend_work(nCells, nEdges, nVertices, nVertLevels_dumm !DIR$ IVDEP do k=1,nVertLevels ! tend_u(k,iEdge) = tend_u(k,iEdge) + tend_u_euler(k,iEdge) - tend_u(k,iEdge) = tend_u(k,iEdge) + tend_u_euler(k,iEdge) + tend_ru_physics(k,iEdge) + tend_u(k,iEdge) = tend_u(k,iEdge) + tend_u_euler(k,iEdge) + tend_ru_physics(k,iEdge)*rTildeEdge(k,iEdge) end do end do