diff --git a/CONTRIBUTORS.md b/CONTRIBUTORS.md index da42a74b5..7d96743f2 100644 --- a/CONTRIBUTORS.md +++ b/CONTRIBUTORS.md @@ -5,3 +5,4 @@ | james-bruten-mo | James Bruten | Met Office | 2025-12-09 | | jennyhickson | Jenny Hickson | Met Office | 2025-12-10 | | mo-marqh | mark Hedley | Met Office | 2025-12-11 | +| tommbendall | Thomas Bendall | Met Office | 2026-01-28 | \ No newline at end of file diff --git a/interfaces/physics_schemes_interface/source/algorithm/bl_exp_alg_mod.x90 b/interfaces/physics_schemes_interface/source/algorithm/bl_exp_alg_mod.x90 index a514beb49..cf9381f37 100644 --- a/interfaces/physics_schemes_interface/source/algorithm/bl_exp_alg_mod.x90 +++ b/interfaces/physics_schemes_interface/source/algorithm/bl_exp_alg_mod.x90 @@ -19,7 +19,9 @@ module bl_exp_alg_mod use fs_continuity_mod, only: W2, W3, Wtheta use sci_geometric_constants_mod, only: get_height_fv, get_da_at_w2, & get_delta_at_wtheta, & - get_dz_at_wtheta + get_dz_at_wtheta, & + get_face_selector_ew, & + get_face_selector_ns use physics_constants_mod, only: get_max_diff, get_rdz_fd1, & get_rdz_w3, get_dtrdz_wth use map_fd_to_prognostics_alg_mod, only: set_wind @@ -208,6 +210,9 @@ contains type( field_type ), pointer :: w2_rmultiplicity type( field_type ), pointer :: dz_wth + type( integer_field_type ), pointer :: face_selector_ew + type( integer_field_type ), pointer :: face_selector_ns + ! local variables type(mesh_type), pointer :: mesh type( field_type ) :: dtl_mphys, dmt_mphys, ngstress_w2, & @@ -321,6 +326,8 @@ contains height_wth => get_height_fv( Wtheta, mesh%get_id() ) height_w3 => get_height_fv( W3, mesh%get_id() ) + face_selector_ew => get_face_selector_ew( mesh%get_id() ) + face_selector_ns => get_face_selector_ns( mesh%get_id() ) ! delta is calculated in sci_calc_delta_at_wtheta_kernel_mod.F90 ! in lfric core as the mininum of the cell horizontal edge lengths delta => get_delta_at_wtheta( mesh%get_id() ) @@ -428,7 +435,8 @@ contains call invoke(bl_exp_du_kernel_type(4_i_def, tau_w2, tau_land_w2, tau_ssi_w2,& rhokm_w2, rdz_fd1, u_physics, & surf_interp_w2, ngstress_w2, fd_tau_w2, & - sea_current_w2_ptr) ) + sea_current_w2_ptr, & + face_selector_ew, face_selector_ns) ) if (l_calc_tau_at_p) then call dz_wth%copy_field_properties(rdz_wth) @@ -446,11 +454,13 @@ contains bl_exp_du_kernel_type(1_i_def, taux, taux_land, taux_ssi, & rhokm_bl, rdz_wth, u_in_w3, & surf_interp, ngstress_bl, fd_taux, & - sea_u_current_ptr), & + sea_u_current_ptr, & + face_selector_ew, face_selector_ns), & bl_exp_du_kernel_type(1_i_def, tauy, tauy_land, tauy_ssi, & rhokm_bl, rdz_wth, v_in_w3, & surf_interp, ngstress_bl, fd_tauy, & - sea_v_current_ptr) ) + sea_v_current_ptr, & + face_selector_ew, face_selector_ns) ) end if if ( subroutine_timers ) call timer('bl_exp_alg') diff --git a/interfaces/physics_schemes_interface/source/algorithm/bl_imp_alg_mod.x90 b/interfaces/physics_schemes_interface/source/algorithm/bl_imp_alg_mod.x90 index db83e0d3d..63c921545 100644 --- a/interfaces/physics_schemes_interface/source/algorithm/bl_imp_alg_mod.x90 +++ b/interfaces/physics_schemes_interface/source/algorithm/bl_imp_alg_mod.x90 @@ -13,7 +13,9 @@ module bl_imp_alg_mod use integer_field_mod, only: integer_field_type use field_collection_mod, only: field_collection_type use fs_continuity_mod, only: W1, W2, W3, Wtheta - use sci_geometric_constants_mod, only: get_height_fv, get_da_at_w2 + use sci_geometric_constants_mod, only: get_height_fv, get_da_at_w2, & + get_face_selector_ew, & + get_face_selector_ns use physics_constants_mod, only: get_rdz_fd1, get_dtrdz_fd2, & get_rdz_w3 use mr_indices_mod, only: nummr, imr_v, imr_cl, imr_ci, & @@ -198,6 +200,9 @@ contains type( field_type ), pointer :: dtrdz => null() type( field_type ), pointer :: rdz_w3 => null() + type( integer_field_type ), pointer :: face_selector_ew + type( integer_field_type ), pointer :: face_selector_ns + ! Local variables type( field_type ) :: u_physics_star type( field_type ) :: du_conv_w2, dw_conv @@ -312,6 +317,9 @@ contains mesh => zh%get_mesh() + face_selector_ew => get_face_selector_ew(mesh%get_id()) + face_selector_ns => get_face_selector_ns(mesh%get_id()) + ! Set up the BL increments call theta%copy_field_properties(dtheta_bl) @@ -367,7 +375,8 @@ contains rhokm_w2, rdz, dtrdz, wetrho_in_w2, & u_physics, u_physics_star, & surf_interp_w2, du_conv_w2, dw_bl, dA, & - height_w1, height_w2), & + height_w1, height_w2, & + face_selector_ew, face_selector_ns), & ! Map the wind back into dynamics quantities inc_X_times_Y(du_bl, dA), & inc_X_times_Y(dissip, dA) ) diff --git a/interfaces/physics_schemes_interface/source/kernel/bl_exp_du_kernel_mod.F90 b/interfaces/physics_schemes_interface/source/kernel/bl_exp_du_kernel_mod.F90 index cd98c4507..8d1b43373 100644 --- a/interfaces/physics_schemes_interface/source/kernel/bl_exp_du_kernel_mod.F90 +++ b/interfaces/physics_schemes_interface/source/kernel/bl_exp_du_kernel_mod.F90 @@ -13,12 +13,14 @@ module bl_exp_du_kernel_mod ANY_SPACE_1, ANY_SPACE_2, & ANY_SPACE_3, ANY_SPACE_4, & ANY_SPACE_5, ANY_SPACE_6, & + ANY_DISCONTINUOUS_SPACE_3, & GH_REAL, GH_WRITE, GH_SCALAR, GH_INTEGER use constants_mod, only: r_def, i_def, r_bl use fs_continuity_mod, only: W1, W2 use kernel_mod, only: kernel_type use nlsizes_namelist_mod, only: bl_levels use jules_surface_config_mod, only: formdrag, formdrag_dist_drag + use reference_element_mod, only: N implicit none @@ -30,19 +32,21 @@ module bl_exp_du_kernel_mod !> Kernel metadata type. type, public, extends(kernel_type) :: bl_exp_du_kernel_type private - type(arg_type) :: meta_args(11) = (/ & - arg_type(GH_SCALAR, GH_INTEGER, GH_READ), &! nfaces - arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_SPACE_1), &! tau w3/w2 - arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_SPACE_2), &! tau_land 2d w3/w2 - arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_SPACE_2), &! tau_ssi 2d w3/w2 - arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_SPACE_3), &! rhokm wth/w2 - arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_SPACE_4), &! rdz wth/fd1 - arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_SPACE_1), &! u_physics w3/w2 - arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_SPACE_5), &! surf_interp mult w3/w2 - arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_SPACE_3), &! ngstress wth/w2 - arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_SPACE_1), &! fd_tau w3/w2 - arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_SPACE_6) &! sea_current mix w3/w2 - /) + type(arg_type) :: meta_args(13) = (/ & + arg_type(GH_SCALAR, GH_INTEGER, GH_READ), &! nfaces + arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_SPACE_1), &! tau w3/w2 + arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_SPACE_2), &! tau_land 2d w3/w2 + arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_SPACE_2), &! tau_ssi 2d w3/w2 + arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_SPACE_3), &! rhokm wth/w2 + arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_SPACE_4), &! rdz wth/fd1 + arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_SPACE_1), &! u_physics w3/w2 + arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_SPACE_5), &! surf_interp mult w3/w2 + arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_SPACE_3), &! ngstress wth/w2 + arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_SPACE_1), &! fd_tau w3/w2 + arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_SPACE_6), &! sea_current mix w3/w2 + arg_type(GH_FIELD, GH_INTEGER, GH_READ, ANY_DISCONTINUOUS_SPACE_3), &! face_selector_ew + arg_type(GH_FIELD, GH_INTEGER, GH_READ, ANY_DISCONTINUOUS_SPACE_3) &! face_selector_ns + /) integer :: operates_on = CELL_COLUMN contains procedure, nopass :: bl_exp_du_code @@ -69,6 +73,10 @@ module bl_exp_du_kernel_mod !> @param[in] ngstress NG stress function !> @param[in] fd_tau Stress from turbulent form-drag !> @param[in] sea_current Ocean surface current + !> @param[in] face_selector_ew 2D field indicating which W/E faces + !! to loop over in this column + !> @param[in] face_selector_ns 2D field indicating which N/S faces + !! to loop over in this column !> @param[in] ndf_half Number of DOFs per cell for half levels !> @param[in] undf_half Number of unique DOFs for half levels !> @param[in] map_half Dofmap for the cell at the base of the column @@ -87,36 +95,44 @@ module bl_exp_du_kernel_mod !> @param[in] ndf_curr Number of DOFs per cell for current space !> @param[in] undf_curr Number of unique DOFs for current space !> @param[in] map_curr Dofmap for the cell at the base of the column - subroutine bl_exp_du_code(nlayers, & - nfaces, & - tau, & - tau_land, & - tau_ssi, & - rhokm, & - rdz, & - u_physics, & - surf_interp, & - ngstress, & - fd_tau, & - sea_current, & - ndf_half, & - undf_half, & - map_half, & - ndf_2d, & - undf_2d, & - map_2d, & - ndf_full, & - undf_full, & - map_full, & - ndf_rdz, & - undf_rdz, & - map_rdz, & - ndf_surf, & - undf_surf, & - map_surf, & - ndf_curr, & - undf_curr, & - map_curr) + !> @param[in] ndf_w3_2d Num of DoFs for 2D W3 per cell + !> @param[in] undf_w3_2d Num of DoFs for this partition for 2D W3 + !> @param[in] map_w3_2d Map for 2D W3 + subroutine bl_exp_du_code(nlayers, & + nfaces, & + tau, & + tau_land, & + tau_ssi, & + rhokm, & + rdz, & + u_physics, & + surf_interp, & + ngstress, & + fd_tau, & + sea_current, & + face_selector_ew, & + face_selector_ns, & + ndf_half, & + undf_half, & + map_half, & + ndf_2d, & + undf_2d, & + map_2d, & + ndf_full, & + undf_full, & + map_full, & + ndf_rdz, & + undf_rdz, & + map_rdz, & + ndf_surf, & + undf_surf, & + map_surf, & + ndf_curr, & + undf_curr, & + map_curr, & + ndf_w3_2d, & + undf_w3_2d, & + map_w3_2d ) !--------------------------------------- ! UM modules containing switches or global constants @@ -141,6 +157,8 @@ subroutine bl_exp_du_code(nlayers, & integer(kind=i_def), intent(in) :: map_full(ndf_full) integer(kind=i_def), intent(in) :: ndf_surf, undf_surf integer(kind=i_def), intent(in) :: map_surf(ndf_surf) + integer(kind=i_def), intent(in) :: ndf_w3_2d, undf_w3_2d + integer(kind=i_def), intent(in) :: map_w3_2d(ndf_w3_2d) real(kind=r_def), dimension(undf_half), intent(inout) :: tau real(kind=r_def), dimension(undf_2d), intent(inout) :: tau_land, tau_ssi @@ -152,70 +170,75 @@ subroutine bl_exp_du_code(nlayers, & real(kind=r_def), dimension(undf_surf), intent(in) :: surf_interp real(kind=r_def), dimension(undf_curr), intent(in) :: sea_current + integer(kind=i_def), dimension(undf_w3_2d), intent(in) :: face_selector_ew + integer(kind=i_def), dimension(undf_w3_2d), intent(in) :: face_selector_ns + ! Internal variables - integer(kind=i_def) :: df, k + integer(kind=i_def) :: df, k, j, total_faces real(kind=r_bl) :: rhokm_land, rhokm_ssi, fland, flandfac, fseafac real(kind=r_bl) :: zh_nonloc(1) real(kind=r_bl), dimension(0:bl_levels-1) :: tau_grad, tau_non_grad, u_sp, & rdz_sp, rhokm_sp, ngstress_sp, tau_sp, fd_tau_sp + total_faces = MIN( & + nfaces, face_selector_ew(map_w3_2d(1)) + face_selector_ns(map_w3_2d(1)) & + ) + ! loop over all faces of the cell - do df = 1,nfaces - - ! Only calculate face if it's not already been done - if (tau(map_half(df)) == 0.0_r_def) then - - !================================================================ - ! In the UM this happens in bdy_expl3 - !================================================================ - fland = surf_interp(map_surf(df) + 0) - rhokm_land = surf_interp(map_surf(df) + 1) - rhokm_ssi = surf_interp(map_surf(df) + 2) - flandfac = surf_interp(map_surf(df) + 3) - fseafac = surf_interp(map_surf(df) + 4) - zh_nonloc = surf_interp(map_surf(df) + 5) - - do k = 0, bl_levels-1 - u_sp(k) = u_physics(map_half(df)+k) - rdz_sp(k) = rdz(map_rdz(df)+k) - rhokm_sp(k) = rhokm(map_full(df)+k) - ngstress_sp(k) = ngstress(map_full(df)+k) - fd_tau_sp(k) = fd_tau(map_half(df)+k) - tau_sp(k) = tau(map_half(df)+k) - end do - - tau_land(map_2d(df)) = rhokm_land * u_sp(0) * flandfac - - ! If sea surface current (sea_current) has not been obtained - ! from coupling fields or from an ancillary file then it will simply - ! contain uniform zeros. - tau_ssi(map_2d(df)) = rhokm_ssi * & - (u_sp(0) - sea_current(map_curr(df))) * fseafac - tau_sp(0) = fland * tau_land(map_2d(df)) & - + (1.0_r_bl - fland) * tau_ssi(map_2d(df)) - - if (formdrag == formdrag_dist_drag) then - if (fland > 0.0_r_bl) then - tau_land(map_2d(df)) = tau_land(map_2d(df)) & - + fd_tau(map_half(df)) / fland - end if + do j = 1, total_faces + df = j + if (j == 3 .and. face_selector_ns(map_w3_2d(1)) == 2 & + .and. face_selector_ew(map_w3_2d(1)) == 1) df = N + + !================================================================ + ! In the UM this happens in bdy_expl3 + !================================================================ + fland = surf_interp(map_surf(df) + 0) + rhokm_land = surf_interp(map_surf(df) + 1) + rhokm_ssi = surf_interp(map_surf(df) + 2) + flandfac = surf_interp(map_surf(df) + 3) + fseafac = surf_interp(map_surf(df) + 4) + zh_nonloc = surf_interp(map_surf(df) + 5) + + do k = 0, bl_levels-1 + u_sp(k) = u_physics(map_half(df)+k) + rdz_sp(k) = rdz(map_rdz(df)+k) + rhokm_sp(k) = rhokm(map_full(df)+k) + ngstress_sp(k) = ngstress(map_full(df)+k) + fd_tau_sp(k) = fd_tau(map_half(df)+k) + tau_sp(k) = tau(map_half(df)+k) + end do + + tau_land(map_2d(df)) = rhokm_land * u_sp(0) * flandfac + + ! If sea surface current (sea_current) has not been obtained + ! from coupling fields or from an ancillary file then it will simply + ! contain uniform zeros. + tau_ssi(map_2d(df)) = rhokm_ssi * & + (u_sp(0) - sea_current(map_curr(df))) * fseafac + tau_sp(0) = fland * tau_land(map_2d(df)) & + + (1.0_r_bl - fland) * tau_ssi(map_2d(df)) + + if (formdrag == formdrag_dist_drag) then + if (fland > 0.0_r_bl) then + tau_land(map_2d(df)) = tau_land(map_2d(df)) & + + fd_tau(map_half(df)) / fland end if - - call ex_flux_uv(pdims, pdims, pdims, bl_levels, & - u_sp(0:bl_levels-1), & - zh_nonloc, & - rdz_sp(1:bl_levels-1), & - rhokm_sp(0:bl_levels-1), & - ngstress_sp(1:bl_levels-1), & - fd_tau_sp(0:bl_levels-1), & - tau_sp(0:bl_levels-1), & - tau_grad(0:bl_levels-1), tau_non_grad(0:bl_levels-1)) - - do k = 0, bl_levels-1 - tau(map_half(df)+k) = tau_sp(k) - end do - - end if ! this face needs calculating + end if + + call ex_flux_uv(pdims, pdims, pdims, bl_levels, & + u_sp(0:bl_levels-1), & + zh_nonloc, & + rdz_sp(1:bl_levels-1), & + rhokm_sp(0:bl_levels-1), & + ngstress_sp(1:bl_levels-1), & + fd_tau_sp(0:bl_levels-1), & + tau_sp(0:bl_levels-1), & + tau_grad(0:bl_levels-1), tau_non_grad(0:bl_levels-1)) + + do k = 0, bl_levels-1 + tau(map_half(df)+k) = tau_sp(k) + end do end do ! loop over df diff --git a/interfaces/physics_schemes_interface/source/kernel/bl_imp_du_kernel_mod.F90 b/interfaces/physics_schemes_interface/source/kernel/bl_imp_du_kernel_mod.F90 index 490e9b597..ba7e9d081 100644 --- a/interfaces/physics_schemes_interface/source/kernel/bl_imp_du_kernel_mod.F90 +++ b/interfaces/physics_schemes_interface/source/kernel/bl_imp_du_kernel_mod.F90 @@ -8,9 +8,11 @@ module bl_imp_du_kernel_mod use kernel_mod, only: kernel_type - use argument_mod, only: arg_type, func_type, & - GH_FIELD, GH_READ, CELL_COLUMN, & - ANY_SPACE_1, ANY_SPACE_2, & + use argument_mod, only: arg_type, func_type, & + GH_FIELD, GH_READ, CELL_COLUMN, & + ANY_SPACE_1, ANY_SPACE_2, & + ANY_DISCONTINUOUS_SPACE_3, & + ANY_DISCONTINUOUS_SPACE_2, & GH_INTEGER, GH_REAL, GH_SCALAR, GH_WRITE use constants_mod, only: r_def, i_def, r_bl use extrusion_config_mod, only: planet_radius @@ -19,6 +21,7 @@ module bl_imp_du_kernel_mod use nlsizes_namelist_mod, only: bl_levels use timestepping_config_mod, only: outer_iterations use blayer_config_mod, only: fric_heating, bl_mix_w + use reference_element_mod, only: N implicit none @@ -32,29 +35,31 @@ module bl_imp_du_kernel_mod !> Kernel metadata type. type, public, extends(kernel_type) :: bl_imp_du_kernel_type private - type(arg_type) :: meta_args(21) = (/ & - arg_type(GH_SCALAR, GH_INTEGER, GH_READ), &! outer - arg_type(GH_FIELD, GH_REAL, GH_WRITE, W2), &! du_bl - arg_type(GH_FIELD, GH_REAL, GH_WRITE, W2), &! dissip - arg_type(GH_FIELD, GH_REAL, GH_WRITE, W2), &! tau - arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_SPACE_1), &! wind10m - arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_SPACE_1), &! wind10m_neut - arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_SPACE_1), &! tau_land - arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_SPACE_1), &! tau_ssi - arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_SPACE_1), &! pseudotau - arg_type(GH_FIELD, GH_REAL, GH_READ, W2), &! rhokm - arg_type(GH_FIELD, GH_REAL, GH_READ, W1), &! rdz - arg_type(GH_FIELD, GH_REAL, GH_READ, W2), &! dtrdz - arg_type(GH_FIELD, GH_REAL, GH_READ, W2), &! wetrho - arg_type(GH_FIELD, GH_REAL, GH_READ, W2), &! u_physics - arg_type(GH_FIELD, GH_REAL, GH_READ, W2), &! u_physics_star - arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_SPACE_2), &! surf_interp - arg_type(GH_FIELD, GH_REAL, GH_READ, W2), &! du_conv - arg_type(GH_FIELD, GH_REAL, GH_READ, WTheta), &! dw_bl - arg_type(GH_FIELD, GH_REAL, GH_READ, W2), &! dA - arg_type(GH_FIELD, GH_REAL, GH_READ, W1), &! height_w1 - arg_type(GH_FIELD, GH_REAL, GH_READ, W2) &! height_w2 - /) + type(arg_type) :: meta_args(23) = (/ & + arg_type(GH_SCALAR, GH_INTEGER, GH_READ), &! outer + arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_2), &! du_bl + arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_2), &! dissip + arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_2), &! tau + arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_SPACE_1), &! wind10m + arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_SPACE_1), &! wind10m_neut + arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_SPACE_1), &! tau_land + arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_SPACE_1), &! tau_ssi + arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_SPACE_1), &! pseudotau + arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_2), &! rhokm + arg_type(GH_FIELD, GH_REAL, GH_READ, W1), &! rdz + arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_2), &! dtrdz + arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_2), &! wetrho + arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_2), &! u_physics + arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_2), &! u_physics_star + arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_SPACE_2), &! surf_interp + arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_2), &! du_conv + arg_type(GH_FIELD, GH_REAL, GH_READ, WTheta), &! dw_bl + arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_2), &! dA + arg_type(GH_FIELD, GH_REAL, GH_READ, W1), &! height_w1 + arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_2), &! height_w2 + arg_type(GH_FIELD, GH_INTEGER, GH_READ, ANY_DISCONTINUOUS_SPACE_3), &! face_selector_ew + arg_type(GH_FIELD, GH_INTEGER, GH_READ, ANY_DISCONTINUOUS_SPACE_3) &! face_selector_ns + /) integer :: operates_on = CELL_COLUMN contains procedure, nopass :: bl_imp_du_code @@ -91,6 +96,10 @@ module bl_imp_du_kernel_mod !> @param[in] dA Area of faces !> @param[in] height_w1 Height of cell top/bottom above surface !> @param[in] height_w2 Height of cell centre above surface + !> @param[in] face_selector_ew 2D field indicating which W/E faces + !! to loop over in this column + !> @param[in] face_selector_ns 2D field indicating which N/S faces + !! to loop over in this column !> @param[in] ndf_w2 Number of DOFs per cell for w2 space !> @param[in] undf_w2 Number of unique DOFs for w2 space !> @param[in] map_w2 Dofmap for the cell at the base of the column @@ -106,43 +115,51 @@ module bl_imp_du_kernel_mod !> @param[in] ndf_wth No of DOFs per cell for wtheta space !> @param[in] undf_wth No of unique DOFs for wtheta space !> @param[in] map_wth DOFmap for cell at base of wtheta column - subroutine bl_imp_du_code(nlayers, & - outer, & - du_bl, & - dissip, & - tau, & - wind10m, & - wind10m_neut, & - tau_land, & - tau_ssi, & - pseudotau, & - rhokm, & - rdz, & - dtrdz, & - wetrho, & - u_physics, & - u_physics_star,& - surf_interp, & - du_conv, & - dw_bl, & - dA, & - height_w1, & - height_w2, & - ndf_w2, & - undf_w2, & - map_w2, & - ndf_w2_2d, & - undf_w2_2d, & - map_w2_2d, & - ndf_w1, & - undf_w1, & - map_w1, & - ndf_w2_surf, & - undf_w2_surf, & - map_w2_surf, & - ndf_wth, & - undf_wth, & - map_wth) + !> @param[in] ndf_w3_2d Num of DoFs for 2D W3 per cell + !> @param[in] undf_w3_2d Num of DoFs for this partition for 2D W3 + !> @param[in] map_w3_2d Map for 2D W3 + subroutine bl_imp_du_code(nlayers, & + outer, & + du_bl, & + dissip, & + tau, & + wind10m, & + wind10m_neut, & + tau_land, & + tau_ssi, & + pseudotau, & + rhokm, & + rdz, & + dtrdz, & + wetrho, & + u_physics, & + u_physics_star, & + surf_interp, & + du_conv, & + dw_bl, & + dA, & + height_w1, & + height_w2, & + face_selector_ew, & + face_selector_ns, & + ndf_w2, & + undf_w2, & + map_w2, & + ndf_w2_2d, & + undf_w2_2d, & + map_w2_2d, & + ndf_w1, & + undf_w1, & + map_w1, & + ndf_w2_surf, & + undf_w2_surf, & + map_w2_surf, & + ndf_wth, & + undf_wth, & + map_wth, & + ndf_w3_2d, & + undf_w3_2d, & + map_w3_2d) !--------------------------------------- ! UM modules containing switches or global constants @@ -164,6 +181,8 @@ subroutine bl_imp_du_code(nlayers, & integer(kind=i_def), intent(in) :: map_w1(ndf_w1) integer(kind=i_def), intent(in) :: ndf_wth, undf_wth integer(kind=i_def), intent(in) :: map_wth(ndf_wth) + integer(kind=i_def), intent(in) :: ndf_w3_2d, undf_w3_2d + integer(kind=i_def), intent(in) :: map_w3_2d(ndf_w3_2d) real(kind=r_def), dimension(undf_w2), intent(inout) :: du_bl, dissip, tau real(kind=r_def), dimension(undf_w2_2d), intent(inout) :: wind10m, & @@ -176,8 +195,11 @@ subroutine bl_imp_du_code(nlayers, & real(kind=r_def), dimension(undf_w1), intent(in) :: height_w1, rdz real(kind=r_def), dimension(undf_wth), intent(in) :: dw_bl + integer(kind=i_def), dimension(undf_w3_2d), intent(in) :: face_selector_ew + integer(kind=i_def), dimension(undf_w3_2d), intent(in) :: face_selector_ns + ! Internal variables - integer(kind=i_def) :: k, df, k_blend + integer(kind=i_def) :: k, j, df, k_blend real(kind=r_bl) :: pnonl, i1, e1, e2, gamma1, gamma2, cdr10m, cdr10m_neut,& du_1, cq_cm_1, tau_land_star, tau_ssi_star, fb_surf, fland, & tau_land_loc, tau_ssi_loc, cd10m_neut @@ -188,211 +210,208 @@ subroutine bl_imp_du_code(nlayers, & ! In the UM this happens in imp_solver - predictor section !================================================================ ! loop over all faces of the cell - do df = 1,4 - - ! Only calculate face if it's not already been done - if (du_bl(map_w2(df)) == 0.0_r_def .and. & - du_bl(map_w2(df)+1) == 0.0_r_def) then - - fland = surf_interp(map_w2_surf(df) + 0) - fb_surf = surf_interp(map_w2_surf(df) + 6) - k_blend = int(surf_interp(map_w2_surf(df) + 7), i_def) - cdr10m = surf_interp(map_w2_surf(df) + 8) - cd10m_neut = surf_interp(map_w2_surf(df) + 9) - cdr10m_neut = surf_interp(map_w2_surf(df) + 10) - - if (fb_surf > 0.0_r_bl) then - pnonl = puns - else - pnonl = pstb - end if - - ! Set implicit weights for predictor section - i1 = (1.0_r_bl + 1.0_r_bl/sqrt2)*(1.0_r_bl + pnonl) - e1 = (1.0_r_bl + 1.0_r_bl/sqrt2)*(pnonl + 1.0_r_bl/sqrt2 + & - sqrt(pnonl*(sqrt2 - 1.0_r_bl) + 0.5_r_bl) ) - e2 = (1.0_r_bl + 1.0_r_bl/sqrt2)*(pnonl + 1.0_r_bl/sqrt2 - & - sqrt(pnonl*(sqrt2 - 1.0_r_bl) + 0.5_r_bl) ) - gamma1 = i1 - gamma2 = i1 - e1 - - ! Take local copies so we don't over-write input values - ! on each outer iteration - tau_loc(0:bl_levels-1) = tau(map_w2(df):map_w2(df)+bl_levels-1) - tau_ssi_loc = tau_ssi(map_w2_2d(df)) - tau_land_loc = tau_land(map_w2_2d(df)) - - !================================================================ - ! In the UM this happens in bdy_impl3 - predictor section - !================================================================ - - do k = 0, bl_levels-1 - du_nt(k) = u_physics_star(map_w2(df) + k) - & - u_physics(map_w2(df) + k) + & - du_conv(map_w2(df) + k) - dtr_rhodz(k) = dtrdz(map_w2(df) + k) / wetrho(map_w2(df) + k) - rhokm_sp(k) = rhokm(map_w2(df) + k) - du_bl_sp(k) = du_bl(map_w2(df) + k) - rdz_sp(k) = rdz(map_w1(df) + k) - end do - - call matrix_sweep(du_bl_sp, du_1, cq_cm, cq_cm_1, df, gamma1, gamma2, & - height_w1, dtr_rhodz, tau_loc, du_nt,rhokm_sp,rdz_sp,& - k_blend, ndf_w1, undf_w1, map_w1) - - !================================================================ - ! In Jules this happens in im_sf_pt2 - predictor section - !================================================================ - if (fland > 0.0_r_bl) then - tau_land_star = ( gamma2 * tau_land_loc + & - gamma1 * rhokm_sp(0) * du_1 ) / & - (1.0_r_bl + gamma1 * rhokm_sp(0) * cq_cm_1) - else - tau_land_star = 0.0_r_bl - end if - if (fland < 1.0_r_bl) then - tau_ssi_star = ( gamma2 * tau_ssi_loc + & - gamma1 * rhokm_sp(0) * du_1 ) / & - (1.0_r_bl + gamma1 * rhokm_sp(0) * cq_cm_1 ) - else - tau_ssi_star = 0.0_r_bl - end if - tau_star(0) = fland * tau_land_star + (1.0_r_bl - fland) * tau_ssi_star - - !================================================================ - ! In the UM this happens in bdy_impl4 - predictor section - !================================================================ - du_bl_sp(0) = du_bl_sp(0) - cq_cm(0) * tau_star(0) + do j = 1, face_selector_ew(map_w3_2d(1)) + face_selector_ns(map_w3_2d(1)) + df = j + if (j == 3 .and. face_selector_ns(map_w3_2d(1)) == 2 & + .and. face_selector_ew(map_w3_2d(1)) == 1) df = N + + fland = surf_interp(map_w2_surf(df) + 0) + fb_surf = surf_interp(map_w2_surf(df) + 6) + k_blend = int(surf_interp(map_w2_surf(df) + 7), i_def) + cdr10m = surf_interp(map_w2_surf(df) + 8) + cd10m_neut = surf_interp(map_w2_surf(df) + 9) + cdr10m_neut = surf_interp(map_w2_surf(df) + 10) + + if (fb_surf > 0.0_r_bl) then + pnonl = puns + else + pnonl = pstb + end if + + ! Set implicit weights for predictor section + i1 = (1.0_r_bl + 1.0_r_bl/sqrt2)*(1.0_r_bl + pnonl) + e1 = (1.0_r_bl + 1.0_r_bl/sqrt2)*(pnonl + 1.0_r_bl/sqrt2 + & + sqrt(pnonl*(sqrt2 - 1.0_r_bl) + 0.5_r_bl) ) + e2 = (1.0_r_bl + 1.0_r_bl/sqrt2)*(pnonl + 1.0_r_bl/sqrt2 - & + sqrt(pnonl*(sqrt2 - 1.0_r_bl) + 0.5_r_bl) ) + gamma1 = i1 + gamma2 = i1 - e1 + + ! Take local copies so we don't over-write input values + ! on each outer iteration + tau_loc(0:bl_levels-1) = tau(map_w2(df):map_w2(df)+bl_levels-1) + tau_ssi_loc = tau_ssi(map_w2_2d(df)) + tau_land_loc = tau_land(map_w2_2d(df)) + + !================================================================ + ! In the UM this happens in bdy_impl3 - predictor section + !================================================================ + + do k = 0, bl_levels-1 + du_nt(k) = u_physics_star(map_w2(df) + k) - & + u_physics(map_w2(df) + k) + & + du_conv(map_w2(df) + k) + dtr_rhodz(k) = dtrdz(map_w2(df) + k) / wetrho(map_w2(df) + k) + rhokm_sp(k) = rhokm(map_w2(df) + k) + du_bl_sp(k) = du_bl(map_w2(df) + k) + rdz_sp(k) = rdz(map_w1(df) + k) + end do - do k = 1, bl_levels-1 - du_bl_sp(k) = du_bl_sp(k) - & - cq_cm(k) * du_bl_sp(k-1) - tau_star(k) = gamma2 * tau_loc(k) + & - gamma1 * rhokm_sp(k) * rdz_sp(k) * & - ( du_bl_sp(k) - du_bl_sp(k-1) ) - end do + call matrix_sweep(du_bl_sp, du_1, cq_cm, cq_cm_1, df, gamma1, gamma2, & + height_w1, dtr_rhodz, tau_loc, du_nt,rhokm_sp,rdz_sp, & + k_blend, ndf_w1, undf_w1, map_w1) + + !================================================================ + ! In Jules this happens in im_sf_pt2 - predictor section + !================================================================ + if (fland > 0.0_r_bl) then + tau_land_star = ( gamma2 * tau_land_loc + & + gamma1 * rhokm_sp(0) * du_1 ) / & + (1.0_r_bl + gamma1 * rhokm_sp(0) * cq_cm_1) + else + tau_land_star = 0.0_r_bl + end if + if (fland < 1.0_r_bl) then + tau_ssi_star = ( gamma2 * tau_ssi_loc + & + gamma1 * rhokm_sp(0) * du_1 ) / & + (1.0_r_bl + gamma1 * rhokm_sp(0) * cq_cm_1 ) + else + tau_ssi_star = 0.0_r_bl + end if + tau_star(0) = fland * tau_land_star + (1.0_r_bl - fland) * tau_ssi_star + + !================================================================ + ! In the UM this happens in bdy_impl4 - predictor section + !================================================================ + du_bl_sp(0) = du_bl_sp(0) - cq_cm(0) * tau_star(0) + + do k = 1, bl_levels-1 + du_bl_sp(k) = du_bl_sp(k) - & + cq_cm(k) * du_bl_sp(k-1) + tau_star(k) = gamma2 * tau_loc(k) + & + gamma1 * rhokm_sp(k) * rdz_sp(k) * & + ( du_bl_sp(k) - du_bl_sp(k-1) ) + end do - do k = 0, bl_levels-1 - du_star(k) = du_bl_sp(k) - end do + do k = 0, bl_levels-1 + du_star(k) = du_bl_sp(k) + end do - !================================================================ - ! In the UM this happens in bdy_impl3 - corrector section - !================================================================ + !================================================================ + ! In the UM this happens in bdy_impl3 - corrector section + !================================================================ - gamma2 = i1 - e2 + gamma2 = i1 - e2 - ! Update explicit fluxes using predictor X* value as needed by the - ! 2nd stage of the scheme. - do k = 1, bl_levels-1 - tau_loc(k) = tau_loc(k) + rhokm_sp(k) * & - ( du_bl_sp(k) - du_bl_sp(k-1) ) & - * rdz_sp(k) - end do + ! Update explicit fluxes using predictor X* value as needed by the + ! 2nd stage of the scheme. + do k = 1, bl_levels-1 + tau_loc(k) = tau_loc(k) + rhokm_sp(k) * & + ( du_bl_sp(k) - du_bl_sp(k-1) ) & + * rdz_sp(k) + end do - call matrix_sweep(du_bl_sp, du_1, cq_cm, cq_cm_1, df, gamma1, gamma2, & - height_w1, dtr_rhodz, tau_loc, du_nt,rhokm_sp,rdz_sp,& - k_blend, ndf_w1, undf_w1, map_w1) - - !================================================================ - ! In Jules this happens in im_sf_pt2 - corrector section - !================================================================ - if (fland > 0.0_r_bl) then - tau_land_loc = ( gamma2 * ( tau_land_loc + & - rhokm_sp(0) * du_star(0) ) & - + gamma1 * rhokm_sp(0) * du_1 ) / & - (1.0_r_bl + gamma1 * rhokm_sp(0) * cq_cm_1 ) - else - tau_land_loc = 0.0_r_bl - end if - if (fland < 1.0_r_bl) then - tau_ssi_loc = ( gamma2 * ( tau_ssi_loc + & - rhokm_sp(0) * du_star(0) ) & - + gamma1 * rhokm_sp(0) * du_1 ) / & + call matrix_sweep(du_bl_sp, du_1, cq_cm, cq_cm_1, df, gamma1, gamma2, & + height_w1, dtr_rhodz, tau_loc, du_nt,rhokm_sp,rdz_sp,& + k_blend, ndf_w1, undf_w1, map_w1) + + !================================================================ + ! In Jules this happens in im_sf_pt2 - corrector section + !================================================================ + if (fland > 0.0_r_bl) then + tau_land_loc = ( gamma2 * ( tau_land_loc + & + rhokm_sp(0) * du_star(0) ) & + + gamma1 * rhokm_sp(0) * du_1 ) / & (1.0_r_bl + gamma1 * rhokm_sp(0) * cq_cm_1 ) - else - tau_ssi_loc = 0.0_r_bl - end if - tau_loc(0) = fland * tau_land_loc + (1.0_r_bl - fland) * tau_ssi_loc - - !================================================================ - ! In the UM this happens in bdy_impl4 - corrector section - !================================================================ - du_bl_sp(0) = du_bl_sp(0) - cq_cm(0) * tau_loc(0) - tau_loc(0) = tau_loc(0) + tau_star(0) + else + tau_land_loc = 0.0_r_bl + end if + if (fland < 1.0_r_bl) then + tau_ssi_loc = ( gamma2 * ( tau_ssi_loc + & + rhokm_sp(0) * du_star(0) ) & + + gamma1 * rhokm_sp(0) * du_1 ) / & + (1.0_r_bl + gamma1 * rhokm_sp(0) * cq_cm_1 ) + else + tau_ssi_loc = 0.0_r_bl + end if + tau_loc(0) = fland * tau_land_loc + (1.0_r_bl - fland) * tau_ssi_loc + + !================================================================ + ! In the UM this happens in bdy_impl4 - corrector section + !================================================================ + du_bl_sp(0) = du_bl_sp(0) - cq_cm(0) * tau_loc(0) + tau_loc(0) = tau_loc(0) + tau_star(0) + + do k = 1, bl_levels-1 + du_bl_sp(k) = du_bl_sp(k) - & + cq_cm(k) * du_bl_sp(k-1) + end do - do k = 1, bl_levels-1 - du_bl_sp(k) = du_bl_sp(k) - & - cq_cm(k) * du_bl_sp(k-1) - end do + do k = 0, bl_levels-1 + du_bl_sp(k) = du_bl_sp(k) + du_star(k) + end do - do k = 0, bl_levels-1 - du_bl_sp(k) = du_bl_sp(k) + du_star(k) - end do + ! Calculate the molecular dissipation rate for frictional heating - ! Calculate the molecular dissipation rate for frictional heating + ! first calculate tau + do k = 1, bl_levels-1 + tau_loc(k) = tau_star(k) + gamma2 * tau_loc(k) + & + gamma1 * rhokm_sp(k) * rdz_sp(k) * & + ( du_bl_sp(k) - du_bl_sp(k-1) ) + end do - ! first calculate tau - do k = 1, bl_levels-1 - tau_loc(k) = tau_star(k) + gamma2 * tau_loc(k) + & - gamma1 * rhokm_sp(k) * rdz_sp(k) * & - ( du_bl_sp(k) - du_bl_sp(k-1) ) + !================================================================ + ! In the UM this happens in imp_solver - end + !================================================================ + if (fric_heating) then + do k = 0, bl_levels-2 + dissip(map_w2(df) + k) = tau_loc(k+1) * rdz_sp(k+1) * & + ( u_physics(map_w2(df) + k+1) + du_bl_sp(k+1) - & + u_physics(map_w2(df) + k) - du_bl_sp(k) ) end do - !================================================================ - ! In the UM this happens in imp_solver - end - !================================================================ - if (fric_heating) then - do k = 0, bl_levels-2 - dissip(map_w2(df) + k) = tau_loc(k+1) * rdz_sp(k+1) * & - ( u_physics(map_w2(df) + k+1) + du_bl_sp(k+1) - & - u_physics(map_w2(df) + k) - du_bl_sp(k) ) - end do - - ! Add on dissipation in the surface layer - dissip(map_w2(df)) = ( tau_loc(0) * & - (u_physics(map_w2(df)) + du_bl_sp(0) ) + & - dissip(map_w2(df)) * & - (height_w2(map_w2(df)+1) - height_w2(map_w2(df))) )/ & - (height_w2(map_w2(df)+1) - height_w1(map_w1(df))) - - ! And nothing on the top level - dissip(map_w2(df)+bl_levels-1) = 0.0_r_bl - end if - - ! Diagnostic calculations at final iteration only - if (outer == outer_iterations) then - tau_ssi(map_w2_2d(df)) = (tau_ssi_loc + tau_ssi_star) & - * dA(map_w2(df)) - tau_land_loc = (tau_land_loc + tau_land_star) & - * dA(map_w2(df)) - wind10m(map_w2_2d(df)) = (u_physics(map_w2(df)) + du_bl_sp(0)) & - * cdr10m * dA(map_w2(df)) - wind10m_neut(map_w2_2d(df)) = & - (u_physics(map_w2(df)) + du_bl_sp(0)) & - * cdr10m_neut * dA(map_w2(df)) - ! tau in bottom level is grid-bix mean surface stress for diagnostics - tau(map_w2(df)) = fland * tau_land_loc + & - (1.0_r_bl - fland) * tau_ssi(map_w2_2d(df)) - ! save end-of-timestep stress profile for diagnostics - do k = 1, bl_levels-1 - tau(map_w2(df) + k) = tau_loc(k)*dA(map_w2(df)) - end do - pseudotau(map_w2_2d(df)) = tau(map_w2(df)) / cd10m_neut - end if - - ! final increment calculation - do k = 0, bl_levels-1 - du_bl(map_w2(df) + k) = du_bl_sp(k) - & - ( u_physics_star(map_w2(df) + k) - u_physics(map_w2(df) + k) ) + ! Add on dissipation in the surface layer + dissip(map_w2(df)) = ( tau_loc(0) * & + (u_physics(map_w2(df)) + du_bl_sp(0) ) + & + dissip(map_w2(df)) * & + (height_w2(map_w2(df)+1) - height_w2(map_w2(df))) )/ & + (height_w2(map_w2(df)+1) - height_w1(map_w1(df))) + + ! And nothing on the top level + dissip(map_w2(df)+bl_levels-1) = 0.0_r_bl + end if + + ! Diagnostic calculations at final iteration only + if (outer == outer_iterations) then + tau_ssi(map_w2_2d(df)) = (tau_ssi_loc + tau_ssi_star) & + * dA(map_w2(df)) + tau_land_loc = (tau_land_loc + tau_land_star) & + * dA(map_w2(df)) + wind10m(map_w2_2d(df)) = (u_physics(map_w2(df)) + du_bl_sp(0)) & + * cdr10m * dA(map_w2(df)) + wind10m_neut(map_w2_2d(df)) = & + (u_physics(map_w2(df)) + du_bl_sp(0)) & + * cdr10m_neut * dA(map_w2(df)) + ! tau in bottom level is grid-bix mean surface stress for diagnostics + tau(map_w2(df)) = fland * tau_land_loc + & + (1.0_r_bl - fland) * tau_ssi(map_w2_2d(df)) + ! save end-of-timestep stress profile for diagnostics + do k = 1, bl_levels-1 + tau(map_w2(df) + k) = tau_loc(k)*dA(map_w2(df)) end do + pseudotau(map_w2_2d(df)) = tau(map_w2(df)) / cd10m_neut + end if - ! Above BL levels it's just the convection increment - do k = bl_levels, nlayers-1 - du_bl(map_w2(df) + k) = du_conv(map_w2(df) + k) - end do + ! final increment calculation + do k = 0, bl_levels-1 + du_bl(map_w2(df) + k) = du_bl_sp(k) - & + ( u_physics_star(map_w2(df) + k) - u_physics(map_w2(df) + k) ) + end do - end if ! this face needs calculating + ! Above BL levels it's just the convection increment + do k = bl_levels, nlayers-1 + du_bl(map_w2(df) + k) = du_conv(map_w2(df) + k) + end do end do ! loop over df diff --git a/rose-stem/site/meto/kgos/gungho_model/azspice/checksum_gungho_model_deep-hot-jupiter-C24_MG_azspice_gnu_fast-debug-64bit.txt b/rose-stem/site/meto/kgos/gungho_model/azspice/checksum_gungho_model_deep-hot-jupiter-C24_MG_azspice_gnu_fast-debug-64bit.txt index 862dd9a4e..0084c94ea 100644 --- a/rose-stem/site/meto/kgos/gungho_model/azspice/checksum_gungho_model_deep-hot-jupiter-C24_MG_azspice_gnu_fast-debug-64bit.txt +++ b/rose-stem/site/meto/kgos/gungho_model/azspice/checksum_gungho_model_deep-hot-jupiter-C24_MG_azspice_gnu_fast-debug-64bit.txt @@ -1,3 +1,3 @@ -Inner product checksum rho = 40EEEF2D423E478C -Inner product checksum theta = 42E4D05A8BCEE836 -Inner product checksum u = 475B46A9C47FB04E +Inner product checksum rho = 40EEEF2D423E47AA +Inner product checksum theta = 42E4D05A8BCEE122 +Inner product checksum u = 475B46A9C47FC554 diff --git a/rose-stem/site/meto/kgos/gungho_model/ex1a/checksum_gungho_model_deep-hot-jupiter-C24_MG_ex1a_gnu_fast-debug-64bit.txt b/rose-stem/site/meto/kgos/gungho_model/ex1a/checksum_gungho_model_deep-hot-jupiter-C24_MG_ex1a_gnu_fast-debug-64bit.txt index abd067848..e22e60d87 100644 --- a/rose-stem/site/meto/kgos/gungho_model/ex1a/checksum_gungho_model_deep-hot-jupiter-C24_MG_ex1a_gnu_fast-debug-64bit.txt +++ b/rose-stem/site/meto/kgos/gungho_model/ex1a/checksum_gungho_model_deep-hot-jupiter-C24_MG_ex1a_gnu_fast-debug-64bit.txt @@ -1,3 +1,3 @@ Inner product checksum rho = 40EEEF2D423E47A9 -Inner product checksum theta = 42E4D05A8BCEE302 -Inner product checksum u = 475B46A9C47FC6D6 +Inner product checksum theta = 42E4D05A8BCEE3DA +Inner product checksum u = 475B46A9C47FBB25 diff --git a/rose-stem/site/meto/kgos/lfric_atm/azspice/checksum_lfric_atm_ral3-seuk_MG_azspice_gnu_fast-debug-32bit.txt b/rose-stem/site/meto/kgos/lfric_atm/azspice/checksum_lfric_atm_ral3-seuk_MG_azspice_gnu_fast-debug-32bit.txt index f22aad7e4..d71a48f21 100644 --- a/rose-stem/site/meto/kgos/lfric_atm/azspice/checksum_lfric_atm_ral3-seuk_MG_azspice_gnu_fast-debug-32bit.txt +++ b/rose-stem/site/meto/kgos/lfric_atm/azspice/checksum_lfric_atm_ral3-seuk_MG_azspice_gnu_fast-debug-32bit.txt @@ -1,9 +1,9 @@ -Inner product checksum rho = 483997CA -Inner product checksum theta = 513548FD -Inner product checksum u = 612F5B34 -Inner product checksum mr1 = 4090F669 -Inner product checksum mr2 = 3589856C -Inner product checksum mr3 = 2FA53E3D -Inner product checksum mr4 = 33F4FD6A -Inner product checksum mr5 = BE7D04E +Inner product checksum rho = 483997C5 +Inner product checksum theta = 513548FF +Inner product checksum u = 612F5B3E +Inner product checksum mr1 = 4090F6B7 +Inner product checksum mr2 = 359C01DF +Inner product checksum mr3 = 2FA585D1 +Inner product checksum mr4 = 33F4FD71 +Inner product checksum mr5 = BE7F5A1 Inner product checksum mr6 = 0 diff --git a/rose-stem/site/meto/kgos/lfric_atm/azspice/checksum_lfric_atm_ral3_ens-seuk_MG_azspice_gnu_fast-debug-32bit.txt b/rose-stem/site/meto/kgos/lfric_atm/azspice/checksum_lfric_atm_ral3_ens-seuk_MG_azspice_gnu_fast-debug-32bit.txt index 42525bc73..c242261da 100644 --- a/rose-stem/site/meto/kgos/lfric_atm/azspice/checksum_lfric_atm_ral3_ens-seuk_MG_azspice_gnu_fast-debug-32bit.txt +++ b/rose-stem/site/meto/kgos/lfric_atm/azspice/checksum_lfric_atm_ral3_ens-seuk_MG_azspice_gnu_fast-debug-32bit.txt @@ -1,9 +1,9 @@ -Inner product checksum rho = 483999D4 -Inner product checksum theta = 51354858 -Inner product checksum u = 612F6179 -Inner product checksum mr1 = 40916617 -Inner product checksum mr2 = 364FE1DF -Inner product checksum mr3 = 2FABFD16 -Inner product checksum mr4 = 33F4FDB5 -Inner product checksum mr5 = BF2197F +Inner product checksum rho = 483999C2 +Inner product checksum theta = 5135485B +Inner product checksum u = 612F617E +Inner product checksum mr1 = 409165ED +Inner product checksum mr2 = 367067AC +Inner product checksum mr3 = 2FAC1DDF +Inner product checksum mr4 = 33F4FDB4 +Inner product checksum mr5 = BF21841 Inner product checksum mr6 = 0 diff --git a/rose-stem/site/meto/kgos/lfric_atm/azspice/checksum_lfric_atm_ral3_mixmol-seuk_MG_azspice_gnu_fast-debug-32bit.txt b/rose-stem/site/meto/kgos/lfric_atm/azspice/checksum_lfric_atm_ral3_mixmol-seuk_MG_azspice_gnu_fast-debug-32bit.txt index cb46157c6..d02f537fc 100644 --- a/rose-stem/site/meto/kgos/lfric_atm/azspice/checksum_lfric_atm_ral3_mixmol-seuk_MG_azspice_gnu_fast-debug-32bit.txt +++ b/rose-stem/site/meto/kgos/lfric_atm/azspice/checksum_lfric_atm_ral3_mixmol-seuk_MG_azspice_gnu_fast-debug-32bit.txt @@ -1,9 +1,9 @@ -Inner product checksum rho = 48398C72 +Inner product checksum rho = 48398C6D Inner product checksum theta = 513542FA -Inner product checksum u = 612EE1BF -Inner product checksum mr1 = 409102F8 -Inner product checksum mr2 = 358C6492 -Inner product checksum mr3 = 300104A8 -Inner product checksum mr4 = 3404ABE7 -Inner product checksum mr5 = C14576D +Inner product checksum u = 612EE1B8 +Inner product checksum mr1 = 40910332 +Inner product checksum mr2 = 358F2478 +Inner product checksum mr3 = 3001071A +Inner product checksum mr4 = 3404ABD5 +Inner product checksum mr5 = C145771 Inner product checksum mr6 = 0 diff --git a/rose-stem/site/meto/kgos/lfric_atm/azspice/checksum_lfric_atm_rce-BiP64x64-1500x1500_MG_azspice_gnu_fast-debug-32bit.txt b/rose-stem/site/meto/kgos/lfric_atm/azspice/checksum_lfric_atm_rce-BiP64x64-1500x1500_MG_azspice_gnu_fast-debug-32bit.txt index c8d190573..265d2a968 100644 --- a/rose-stem/site/meto/kgos/lfric_atm/azspice/checksum_lfric_atm_rce-BiP64x64-1500x1500_MG_azspice_gnu_fast-debug-32bit.txt +++ b/rose-stem/site/meto/kgos/lfric_atm/azspice/checksum_lfric_atm_rce-BiP64x64-1500x1500_MG_azspice_gnu_fast-debug-32bit.txt @@ -1,9 +1,9 @@ -Inner product checksum rho = 4830B5EE -Inner product checksum theta = 513FD8B1 -Inner product checksum u = 5F74A875 -Inner product checksum mr1 = 4155D3E1 -Inner product checksum mr2 = 3A38BFC2 -Inner product checksum mr3 = 309DE635 +Inner product checksum rho = 4830B606 +Inner product checksum theta = 513FD8B9 +Inner product checksum u = 5F74ACFB +Inner product checksum mr1 = 4155C865 +Inner product checksum mr2 = 3A3CCB3E +Inner product checksum mr3 = 3098107A Inner product checksum mr4 = 0 Inner product checksum mr5 = 0 Inner product checksum mr6 = 0 diff --git a/rose-stem/site/meto/kgos/lfric_atm/ex1a/checksum_lfric_atm_comp_tran_ref_3d_l120-BiP64x64-1500x1500_MG_ex1a_cce_fast-debug-32bit.txt b/rose-stem/site/meto/kgos/lfric_atm/ex1a/checksum_lfric_atm_comp_tran_ref_3d_l120-BiP64x64-1500x1500_MG_ex1a_cce_fast-debug-32bit.txt index 91bb01bc5..d79ff4df1 100644 --- a/rose-stem/site/meto/kgos/lfric_atm/ex1a/checksum_lfric_atm_comp_tran_ref_3d_l120-BiP64x64-1500x1500_MG_ex1a_cce_fast-debug-32bit.txt +++ b/rose-stem/site/meto/kgos/lfric_atm/ex1a/checksum_lfric_atm_comp_tran_ref_3d_l120-BiP64x64-1500x1500_MG_ex1a_cce_fast-debug-32bit.txt @@ -1,9 +1,9 @@ -Inner product checksum rho = 489A8317 -Inner product checksum theta = 51BA5074 -Inner product checksum u = 5F1192E8 -Inner product checksum mr1 = 4187BC72 -Inner product checksum mr2 = 3AF9734F -Inner product checksum mr3 = 34BC1253 +Inner product checksum rho = 489A840E +Inner product checksum theta = 51BA5026 +Inner product checksum u = 5F1175D9 +Inner product checksum mr1 = 41877CAC +Inner product checksum mr2 = 3B0397AF +Inner product checksum mr3 = 34F1BEC7 Inner product checksum mr4 = 0 Inner product checksum mr5 = 0 Inner product checksum mr6 = 0 diff --git a/rose-stem/site/meto/kgos/lfric_atm/ex1a/checksum_lfric_atm_ral3-seuk_MG_ex1a_cce_fast-debug-32bit.txt b/rose-stem/site/meto/kgos/lfric_atm/ex1a/checksum_lfric_atm_ral3-seuk_MG_ex1a_cce_fast-debug-32bit.txt index 1e9885931..c38119564 100644 --- a/rose-stem/site/meto/kgos/lfric_atm/ex1a/checksum_lfric_atm_ral3-seuk_MG_ex1a_cce_fast-debug-32bit.txt +++ b/rose-stem/site/meto/kgos/lfric_atm/ex1a/checksum_lfric_atm_ral3-seuk_MG_ex1a_cce_fast-debug-32bit.txt @@ -1,9 +1,9 @@ -Inner product checksum rho = 483997DA -Inner product checksum theta = 5135490A -Inner product checksum u = 612F5E0C -Inner product checksum mr1 = 4090F5EC -Inner product checksum mr2 = 3599E934 -Inner product checksum mr3 = 2FA5315A -Inner product checksum mr4 = 33F4FD4E -Inner product checksum mr5 = BE85FCC +Inner product checksum rho = 483997DF +Inner product checksum theta = 51354909 +Inner product checksum u = 612F5E10 +Inner product checksum mr1 = 4090F616 +Inner product checksum mr2 = 3597B63A +Inner product checksum mr3 = 2FA4E0DE +Inner product checksum mr4 = 33F4FD55 +Inner product checksum mr5 = BE85FEF Inner product checksum mr6 = 0 diff --git a/rose-stem/site/meto/kgos/lfric_atm/ex1a/checksum_lfric_atm_ral3_ens-seuk_MG_ex1a_cce_fast-debug-32bit.txt b/rose-stem/site/meto/kgos/lfric_atm/ex1a/checksum_lfric_atm_ral3_ens-seuk_MG_ex1a_cce_fast-debug-32bit.txt index 9f0b2b052..68fdb02b6 100644 --- a/rose-stem/site/meto/kgos/lfric_atm/ex1a/checksum_lfric_atm_ral3_ens-seuk_MG_ex1a_cce_fast-debug-32bit.txt +++ b/rose-stem/site/meto/kgos/lfric_atm/ex1a/checksum_lfric_atm_ral3_ens-seuk_MG_ex1a_cce_fast-debug-32bit.txt @@ -1,9 +1,9 @@ -Inner product checksum rho = 483997EC -Inner product checksum theta = 513548FD -Inner product checksum u = 612F5EB9 -Inner product checksum mr1 = 40910DEE -Inner product checksum mr2 = 35DE4A2C -Inner product checksum mr3 = 2FA75F15 -Inner product checksum mr4 = 33F4FFD3 -Inner product checksum mr5 = BE04B38 +Inner product checksum rho = 483997EA +Inner product checksum theta = 513548FE +Inner product checksum u = 612F5EAC +Inner product checksum mr1 = 40910E73 +Inner product checksum mr2 = 35E0BD53 +Inner product checksum mr3 = 2FA7843A +Inner product checksum mr4 = 33F4FFDE +Inner product checksum mr5 = BE04B43 Inner product checksum mr6 = 0 diff --git a/rose-stem/site/meto/kgos/lfric_atm/ex1a/checksum_lfric_atm_ral3_mixmol-seuk_MG_ex1a_cce_fast-debug-32bit.txt b/rose-stem/site/meto/kgos/lfric_atm/ex1a/checksum_lfric_atm_ral3_mixmol-seuk_MG_ex1a_cce_fast-debug-32bit.txt index b3bf80080..aca81e36e 100644 --- a/rose-stem/site/meto/kgos/lfric_atm/ex1a/checksum_lfric_atm_ral3_mixmol-seuk_MG_ex1a_cce_fast-debug-32bit.txt +++ b/rose-stem/site/meto/kgos/lfric_atm/ex1a/checksum_lfric_atm_ral3_mixmol-seuk_MG_ex1a_cce_fast-debug-32bit.txt @@ -1,9 +1,9 @@ Inner product checksum rho = 48398C6A Inner product checksum theta = 51354302 -Inner product checksum u = 612EDCBA -Inner product checksum mr1 = 4091034B -Inner product checksum mr2 = 3592F0C2 -Inner product checksum mr3 = 300105E3 -Inner product checksum mr4 = 3404AC1D -Inner product checksum mr5 = C14A22D +Inner product checksum u = 612EDCB9 +Inner product checksum mr1 = 4091035B +Inner product checksum mr2 = 35894423 +Inner product checksum mr3 = 30010BA1 +Inner product checksum mr4 = 3404AC29 +Inner product checksum mr5 = C14A22F Inner product checksum mr6 = 0 diff --git a/rose-stem/site/meto/kgos/lfric_atm/ex1a/checksum_lfric_atm_rce-BiP64x64-1500x1500_MG_ex1a_cce_fast-debug-32bit.txt b/rose-stem/site/meto/kgos/lfric_atm/ex1a/checksum_lfric_atm_rce-BiP64x64-1500x1500_MG_ex1a_cce_fast-debug-32bit.txt index 9006b2022..367dad278 100644 --- a/rose-stem/site/meto/kgos/lfric_atm/ex1a/checksum_lfric_atm_rce-BiP64x64-1500x1500_MG_ex1a_cce_fast-debug-32bit.txt +++ b/rose-stem/site/meto/kgos/lfric_atm/ex1a/checksum_lfric_atm_rce-BiP64x64-1500x1500_MG_ex1a_cce_fast-debug-32bit.txt @@ -1,9 +1,9 @@ -Inner product checksum rho = 4830B646 -Inner product checksum theta = 513FD832 -Inner product checksum u = 5F74BAEC -Inner product checksum mr1 = 4155B86C -Inner product checksum mr2 = 3A48BCB8 -Inner product checksum mr3 = 30A6CEC1 +Inner product checksum rho = 4830B62F +Inner product checksum theta = 513FD82E +Inner product checksum u = 5F74B33C +Inner product checksum mr1 = 4155B6EF +Inner product checksum mr2 = 3A44F0C2 +Inner product checksum mr3 = 30953C70 Inner product checksum mr4 = 0 Inner product checksum mr5 = 0 Inner product checksum mr6 = 0 diff --git a/science/gungho/source/algorithm/diffusion/leonard_term_alg_mod.x90 b/science/gungho/source/algorithm/diffusion/leonard_term_alg_mod.x90 index c9174a68e..b2109a54a 100644 --- a/science/gungho/source/algorithm/diffusion/leonard_term_alg_mod.x90 +++ b/science/gungho/source/algorithm/diffusion/leonard_term_alg_mod.x90 @@ -10,13 +10,16 @@ module leonard_term_alg_mod use constants_mod, only: i_def, r_def use field_mod, only: field_type use field_collection_mod, only: field_collection_type + use integer_field_mod, only: integer_field_type use mr_indices_mod, only: nummr, imr_v, imr_cl use formulation_config_mod, only: moisture_formulation, & moisture_formulation_dry use fs_continuity_mod, only: W1, W2, W3, Wtheta use sci_geometric_constants_mod, & - only: get_height_fv, get_panel_id + only: get_height_fv, get_panel_id, & + get_face_selector_ew, & + get_face_selector_ns use physics_constants_mod, only: get_dtrdz_fd2 use mixing_config_mod, only: leonard_kl @@ -91,6 +94,9 @@ subroutine leonard_term_alg(mt_inc_leonard, thetal_inc_leonard, & type( field_type ), pointer :: height_wth => null() type( field_type ), pointer :: panel_id => null() + type( integer_field_type ), pointer :: face_selector_ew + type( integer_field_type ), pointer :: face_selector_ns + ! local variables integer(kind=i_def) :: mesh_id integer(kind=i_def), parameter :: stencil_depth = 1 @@ -123,6 +129,8 @@ subroutine leonard_term_alg(mt_inc_leonard, thetal_inc_leonard, & height_w2 => get_height_fv(W2, mesh_id) height_w1 => get_height_fv(W1, mesh_id) panel_id => get_panel_id(mesh_id) + face_selector_ew => get_face_selector_ew(mesh_id) + face_selector_ns => get_face_selector_ns(mesh_id) ! Set-up arrays for local fields call rho%copy_field_properties(kl) @@ -173,6 +181,8 @@ subroutine leonard_term_alg(mt_inc_leonard, thetal_inc_leonard, & wetrho_in_w2, & panel_id, & stencil_depth, & + face_selector_ew, & + face_selector_ns, & planet_radius, & leonard_kl, & dt, bl_levels ), & diff --git a/science/gungho/source/algorithm/diffusion/mixing_alg_mod.x90 b/science/gungho/source/algorithm/diffusion/mixing_alg_mod.x90 index e24f3487d..f14b46d7c 100644 --- a/science/gungho/source/algorithm/diffusion/mixing_alg_mod.x90 +++ b/science/gungho/source/algorithm/diffusion/mixing_alg_mod.x90 @@ -29,8 +29,11 @@ contains use log_mod, only: log_event, & LOG_LEVEL_INFO use field_mod, only: field_type + use integer_field_mod, only: integer_field_type use field_collection_mod, only: field_collection_type - use sci_geometric_constants_mod, only: get_dx_at_w2 + use sci_geometric_constants_mod, only: get_dx_at_w2, & + get_face_selector_ew, & + get_face_selector_ns use mesh_mod, only: mesh_type use mr_indices_mod, only: nummr use mixing_config_mod, only: viscosity, smagorinsky, viscosity_mu @@ -55,9 +58,10 @@ contains ! Increment fields type( field_type ) :: du, dtheta - type( field_type ), pointer :: dx_at_w2 => null() - - type( mesh_type ), pointer :: mesh => null() + type( field_type ), pointer :: dx_at_w2 + type( integer_field_type ), pointer :: face_selector_ew + type( integer_field_type ), pointer :: face_selector_ns + type( mesh_type ), pointer :: mesh !-------------------------------------------------------------------- ! Apply horizontal Smagorinsky subgrid mixing @@ -87,10 +91,15 @@ contains call du%initialise( u%get_function_space() ) dx_at_w2 => get_dx_at_w2(mesh%get_id()) + face_selector_ew => get_face_selector_ew(mesh%get_id()) + face_selector_ns => get_face_selector_ns(mesh%get_id()) + call invoke( setval_c(du, 0.0_r_def), & tracer_viscosity_kernel_type( dtheta, theta, 1, dx_at_w2, & viscosity_mu ), & momentum_viscosity_kernel_type( du, u, 1, dx_at_w2, 1, & + face_selector_ew, & + face_selector_ns, & viscosity_mu ), & inc_X_plus_bY( theta, dt, dtheta ), & inc_X_plus_bY( u, dt, du ) ) diff --git a/science/gungho/source/algorithm/diffusion/smagorinsky_alg_mod.x90 b/science/gungho/source/algorithm/diffusion/smagorinsky_alg_mod.x90 index d40e14ab7..78979e46e 100644 --- a/science/gungho/source/algorithm/diffusion/smagorinsky_alg_mod.x90 +++ b/science/gungho/source/algorithm/diffusion/smagorinsky_alg_mod.x90 @@ -10,6 +10,7 @@ module smagorinsky_alg_mod use constants_mod, only: i_def, r_def use field_mod, only: field_type use field_collection_mod, only: field_collection_type + use integer_field_mod, only: integer_field_type use mr_indices_mod, only: nummr, imr_v, imr_cl use formulation_config_mod, only: moisture_formulation, & @@ -33,7 +34,8 @@ module smagorinsky_alg_mod use sci_geometric_constants_mod, & only: get_height_fv, get_panel_id, & get_delta_at_wtheta, get_dx_at_w2, & - get_dz_at_wtheta + get_dz_at_wtheta, get_face_selector_ew, & + get_face_selector_ns use fs_continuity_mod, only: W1, W2 implicit none @@ -84,6 +86,9 @@ subroutine smagorinsky_alg(dtheta_io, du_io, mr, theta, u, & type( field_type ), pointer :: panel_id => null() type( field_type ) :: dmr_v, dmr_cl + type( integer_field_type ), pointer :: face_selector_ew + type( integer_field_type ), pointer :: face_selector_ns + type( mesh_type ), pointer :: mesh => null() ! Stencil depth for the Smagorinsky diffusion kernel @@ -112,6 +117,9 @@ subroutine smagorinsky_alg(dtheta_io, du_io, mr, theta, u, & dx_at_w2 => get_dx_at_w2(mesh%get_id()) panel_id => get_panel_id(mesh%get_id()) + face_selector_ew => get_face_selector_ew(mesh%get_id()) + face_selector_ns => get_face_selector_ns(mesh%get_id()) + if ( boundary_layer == boundary_layer_um ) then ! Use stability-dependent diffusion coefficient from UM BL scheme @@ -153,7 +161,9 @@ subroutine smagorinsky_alg(dtheta_io, du_io, mr, theta, u, & visc_m, & smag_stencil_depth, & panel_id, & - smag_stencil_depth) ) + smag_stencil_depth, & + face_selector_ew, & + face_selector_ns ) ) else ! Without UM BL scheme @@ -213,7 +223,9 @@ subroutine smagorinsky_alg(dtheta_io, du_io, mr, theta, u, & visc_m, & smag_stencil_depth, & panel_id, & - smag_stencil_depth) ) + smag_stencil_depth, & + face_selector_ew, & + face_selector_ns ) ) end if ! With or without BL scheme diff --git a/science/gungho/source/algorithm/physics/external_forcing_alg_mod.X90 b/science/gungho/source/algorithm/physics/external_forcing_alg_mod.X90 index 1a0eb304d..f25997fe6 100644 --- a/science/gungho/source/algorithm/physics/external_forcing_alg_mod.X90 +++ b/science/gungho/source/algorithm/physics/external_forcing_alg_mod.X90 @@ -49,14 +49,16 @@ module external_forcing_alg_mod use constants_mod, only: i_def, r_def, l_def use field_mod, only: field_type use field_collection_mod, only: field_collection_type - + use integer_field_mod, only: integer_field_type use fs_continuity_mod, only: W1, W2, Wtheta use sci_fem_constants_mod, only: get_rmultiplicity_fv - use sci_geometric_constants_mod, only: get_coordinates, & - get_da_at_w2, & - get_panel_id, & - get_height_fv, & - get_dx_at_w2 + use sci_geometric_constants_mod, only: get_coordinates, & + get_da_at_w2, & + get_panel_id, & + get_height_fv, & + get_dx_at_w2, & + get_face_selector_ew, & + get_face_selector_ns use mesh_mod, only: mesh_type use mr_indices_mod, only: nummr use io_config_mod, only: subroutine_timers, write_diag, & @@ -146,6 +148,9 @@ contains type( field_type ), pointer :: height_w2 => null() type( field_type ), pointer :: dx_at_w2 => null() + type( integer_field_type ), pointer :: face_selector_ew + type( integer_field_type ), pointer :: face_selector_ns + type( field_type ) :: u_smoothed, dx, visc_m integer( kind=i_def ) :: order, nlayers @@ -310,6 +315,8 @@ contains height_w2 => get_height_fv(W2, mesh%get_id()) dx_at_w2 => get_dx_at_w2(mesh%get_id()) + face_selector_ew => get_face_selector_ew(mesh%get_id()) + face_selector_ns => get_face_selector_ns(mesh%get_id()) call visc_m%initialise( vector_space = theta%get_function_space()) call dx%initialise( vector_space = dx_at_w2%get_function_space() ) @@ -326,17 +333,19 @@ contains call invoke(setval_X(u_smoothed, u)) do order = 1, diffusion_order - call invoke(momentum_smagorinsky_kernel_type(du_forcing, & - u_smoothed, & - 1, & - dx_at_w2, & - 1, & - height_w2, & - height_w1, & - visc_m, & - 1, & - panel_id, & - 1)) + call invoke(momentum_smagorinsky_kernel_type(du_forcing, & + u_smoothed, & + 1, & + dx_at_w2, & + 1, & + height_w2, & + height_w1, & + visc_m, & + 1, & + panel_id, & + 1, & + face_selector_ew, & + face_selector_ns) ) call invoke(set_any_dof_kernel_type(du_forcing, T, 0.0_r_def)) call invoke(inc_X_plus_Y(u_smoothed, du_forcing), & setval_c(du_forcing, 0.0_r_def)) diff --git a/science/gungho/source/algorithm/timestepping/rk_alg_timestep_mod.x90 b/science/gungho/source/algorithm/timestepping/rk_alg_timestep_mod.x90 index 45605e512..69c3a2fca 100644 --- a/science/gungho/source/algorithm/timestepping/rk_alg_timestep_mod.x90 +++ b/science/gungho/source/algorithm/timestepping/rk_alg_timestep_mod.x90 @@ -29,9 +29,11 @@ module rk_alg_timestep_mod use derived_config_mod, only: bundle_size use sci_fem_constants_mod, only: get_qr_fe, & get_inverse_mass_matrix_fe - use sci_geometric_constants_mod, only: get_coordinates, & - get_panel_id, & - get_dx_at_w2 + use sci_geometric_constants_mod, only: get_coordinates, & + get_panel_id, & + get_dx_at_w2, & + get_face_selector_ew, & + get_face_selector_ns use transport_config_mod, only: operators, & operators_fv use mixing_config_mod, only: viscosity, viscosity_mu @@ -62,6 +64,7 @@ module rk_alg_timestep_mod use field_array_mod, only: field_array_type use field_mod, only: field_type use field_collection_mod, only: field_collection_type + use integer_field_mod, only: integer_field_type use function_space_collection_mod, only: function_space_collection use driver_modeldb_mod, only: modeldb_type use model_clock_mod, only: model_clock_type @@ -322,6 +325,8 @@ contains type(field_type), pointer :: dx_at_w2 type(operator_type), pointer :: m3_inv type(operator_type), pointer :: coriolis + type(integer_field_type), pointer :: face_selector_ew + type(integer_field_type), pointer :: face_selector_ns type(quadrature_rule_gaussian_type) :: gaussian_quadrature @@ -444,11 +449,15 @@ contains if ( viscosity ) then call log_event( 'Applying Viscosity', LOG_LEVEL_INFO ) dx_at_w2 => get_dx_at_w2(mesh_id) + face_selector_ew => get_face_selector_ew(mesh%get_id()) + face_selector_ns => get_face_selector_ns(mesh%get_id()) + call invoke(setval_c(self%inc(igh_u), 0.0_r_def), & tracer_viscosity_kernel_type( self%inc(igh_t), self%state(igh_t), 1, & dx_at_w2, viscosity_mu ), & momentum_viscosity_kernel_type( self%inc(igh_u), self%state(igh_u), 1, & - dx_at_w2, 1, viscosity_mu ), & + dx_at_w2, 1, face_selector_ew, & + face_selector_ns, viscosity_mu ), & setval_c(self%inc(igh_d), 0.0_r_def) ) call bundle_axpy( real(model_clock%get_seconds_per_step(), r_def), & self%inc, self%state, self%state, bundle_size ) diff --git a/science/gungho/source/kernel/diffusion/leonard_term_u_kernel_mod.F90 b/science/gungho/source/kernel/diffusion/leonard_term_u_kernel_mod.F90 index ebc5dae9a..c8e893359 100644 --- a/science/gungho/source/kernel/diffusion/leonard_term_u_kernel_mod.F90 +++ b/science/gungho/source/kernel/diffusion/leonard_term_u_kernel_mod.F90 @@ -14,10 +14,13 @@ module leonard_term_u_kernel_mod GH_READ, GH_WRITE, GH_WRITE, & CELL_COLUMN, STENCIL, REGION, & ANY_DISCONTINUOUS_SPACE_9, & + ANY_DISCONTINUOUS_SPACE_3, & + ANY_DISCONTINUOUS_SPACE_2, & GH_INTEGER use constants_mod, only : r_def, i_def use fs_continuity_mod, only : Wtheta, W2, W1 use kernel_mod, only : kernel_type + use reference_element_mod, only : N implicit none @@ -31,22 +34,25 @@ module leonard_term_u_kernel_mod type, public, extends(kernel_type) :: leonard_term_u_kernel_type private - type(arg_type) :: meta_args(13) = (/ & - arg_type(GH_FIELD, GH_REAL, GH_WRITE, W2), & - arg_type(GH_FIELD, GH_REAL, GH_READ, W2, STENCIL(REGION)), & - arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta, STENCIL(REGION)), & - arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & - arg_type(GH_FIELD, GH_REAL, GH_READ, W2), & - arg_type(GH_FIELD, GH_REAL, GH_READ, W1), & - arg_type(GH_FIELD, GH_REAL, GH_READ, W2), & - arg_type(GH_FIELD, GH_REAL, GH_READ, W2), & - arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_9, & - STENCIL(REGION)), & - arg_type(GH_SCALAR, GH_REAL, GH_READ), & - arg_type(GH_SCALAR, GH_REAL, GH_READ), & - arg_type(GH_SCALAR, GH_REAL, GH_READ), & - arg_type(GH_SCALAR, GH_INTEGER, GH_READ) & - /) + type(arg_type) :: meta_args(15) = (/ & + arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_2), & + arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_2, & + STENCIL(REGION)), & + arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta, STENCIL(REGION)), & + arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & + arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_2), & + arg_type(GH_FIELD, GH_REAL, GH_READ, W1), & + arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_2), & + arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_2), & + arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_9, & + STENCIL(REGION)), & + arg_type(GH_FIELD, GH_INTEGER, GH_READ, ANY_DISCONTINUOUS_SPACE_3), & + arg_type(GH_FIELD, GH_INTEGER, GH_READ, ANY_DISCONTINUOUS_SPACE_3), & + arg_type(GH_SCALAR, GH_REAL, GH_READ), & + arg_type(GH_SCALAR, GH_REAL, GH_READ), & + arg_type(GH_SCALAR, GH_REAL, GH_READ), & + arg_type(GH_SCALAR, GH_INTEGER, GH_READ) & + /) integer :: operates_on = CELL_COLUMN contains procedure, nopass :: leonard_term_u_code @@ -61,42 +67,49 @@ module leonard_term_u_kernel_mod !> @brief Calculates increment due to vertical Leonard term flux !> for variables held in wt space -!! @param[in] nlayers Number of layers in the mesh -!! @param[in,out] u_inc Leonard term increment on w2 -!! @param[in] u_n Wind on w2 -!! @param[in] map_w2_stencil_size Number of cells in the stencil at the base -!! of the column for w2 -!! @param[in] map_w2_stencil Array holding the dofmap for the stencil at the -!! base of the column for w2 -!! @param[in] velocity_w2v Velocity normal to the top face of the cell -!! @param[in] map_wt_stencil_size Number of cells in the stencil at the base -!! of the column for Wtheta -!! @param[in] map_wt_stencil Array holding the dofmap for the stencil at the -!! base of the column for Wtheta -!! @param[in] vel_w2v_inc Leonard term increment of velocity_w2v -!! @param[in] dtrdz_fd2 Array of dt/(r*dz) at FD2 points -!! @param[in] height_w1 Height of w1 space levels above the surface -!! @param[in] height_w2 Height of w2 space levels above the surface -!! @param[in] wetrho_in_w2 Density on w2 space levels -!! @param[in] panel_id The ID number of the current panel -!! @param[in] map_pid_stencil_size Size of the panel ID stencil -!! @param[in] map_pid_stencil Stencil map for the panel ID -!! @param[in] planet_radius The planet radius -!! @param[in] leonard_kl The user-specified Leonard term parameter -!! @param[in] dt The model timestep length -!! @param[in] bl_levels The number of boundary-layer levels -!! @param[in] ndf_w2 Number of degrees of freedom per cell for w2 space -!! @param[in] undf_w2 Number of unique degrees of freedom for w2 space -!! @param[in] map_w2 Cell dofmap for w2 space -!! @param[in] ndf_wt Number of degrees of freedom per cell for theta space -!! @param[in] undf_wt Number of unique degrees of freedom for theta space -!! @param[in] map_wt Cell dofmap for theta space -!! @param[in] ndf_w1 Number of degrees of freedom per cell for w1 space -!! @param[in] undf_w1 Number of unique degrees of freedom for w1 space -!! @param[in] map_w1 Cell dofmap for w1 space -!! @param[in] ndf_pid Number of degrees of freedom per cell for pid space -!! @param[in] undf_pid Number of unique degrees of freedom for pid space -!! @param[in] map_pid Cell dofmap for pid space +!> @param[in] nlayers Number of layers in the mesh +!> @param[in,out] u_inc Leonard term increment on w2 +!> @param[in] u_n Wind on w2 +!> @param[in] map_w2_stencil_size Number of cells in the stencil at the base +!> of the column for w2 +!> @param[in] map_w2_stencil Array holding the dofmap for the stencil at the +!> base of the column for w2 +!> @param[in] velocity_w2v Velocity normal to the top face of the cell +!> @param[in] map_wt_stencil_size Number of cells in the stencil at the base +!> of the column for Wtheta +!> @param[in] map_wt_stencil Array holding the dofmap for the stencil at the +!> base of the column for Wtheta +!> @param[in] vel_w2v_inc Leonard term increment of velocity_w2v +!> @param[in] dtrdz_fd2 Array of dt/(r*dz) at FD2 points +!> @param[in] height_w1 Height of w1 space levels above the surface +!> @param[in] height_w2 Height of w2 space levels above the surface +!> @param[in] wetrho_in_w2 Density on w2 space levels +!> @param[in] panel_id The ID number of the current panel +!> @param[in] map_pid_stencil_size Size of the panel ID stencil +!> @param[in] map_pid_stencil Stencil map for the panel ID +!> @param[in] face_selector_ew 2D field indicating which W/E faces to loop over +!! in this column +!> @param[in] face_selector_ns 2D field indicating which N/S faces to loop over +!! in this column +!> @param[in] planet_radius The planet radius +!> @param[in] leonard_kl The user-specified Leonard term parameter +!> @param[in] dt The model timestep length +!> @param[in] bl_levels The number of boundary-layer levels +!> @param[in] ndf_w2 Number of degrees of freedom per cell for w2 space +!> @param[in] undf_w2 Number of unique degrees of freedom for w2 space +!> @param[in] map_w2 Cell dofmap for w2 space +!> @param[in] ndf_wt Number of degrees of freedom per cell for theta space +!> @param[in] undf_wt Number of unique degrees of freedom for theta space +!> @param[in] map_wt Cell dofmap for theta space +!> @param[in] ndf_w1 Number of degrees of freedom per cell for w1 space +!> @param[in] undf_w1 Number of unique degrees of freedom for w1 space +!> @param[in] map_w1 Cell dofmap for w1 space +!> @param[in] ndf_pid Number of degrees of freedom per cell for pid space +!> @param[in] undf_pid Number of unique degrees of freedom for pid space +!> @param[in] map_pid Cell dofmap for pid space +!> @param[in] ndf_w3_2d Num of DoFs for 2D W3 per cell +!> @param[in] undf_w3_2d Num of DoFs for this partition for 2D W3 +!> @param[in] map_w3_2d Map for 2D W3 subroutine leonard_term_u_code( nlayers, & u_inc, & u_n, & @@ -110,13 +123,15 @@ subroutine leonard_term_u_code( nlayers, & wetrho_in_w2, & panel_id, & map_pid_stencil_size, map_pid_stencil, & + face_selector_ew, face_selector_ns, & planet_radius, & leonard_kl, & dt, bl_levels, & ndf_w2, undf_w2, map_w2, & ndf_wt, undf_wt, map_wt, & ndf_w1, undf_w1, map_w1, & - ndf_pid, undf_pid, map_pid & + ndf_pid, undf_pid, map_pid, & + ndf_w3_2d, undf_w3_2d, map_w3_2d & ) implicit none @@ -127,32 +142,37 @@ subroutine leonard_term_u_code( nlayers, & integer(kind=i_def), intent(in) :: ndf_w2, undf_w2 integer(kind=i_def), intent(in) :: ndf_wt, undf_wt integer(kind=i_def), intent(in) :: ndf_pid, undf_pid + integer(kind=i_def), intent(in) :: ndf_w3_2d, undf_w3_2d integer(kind=i_def), intent(in) :: map_w2_stencil_size integer(kind=i_def), dimension(ndf_w2,map_w2_stencil_size), intent(in) :: map_w2_stencil integer(kind=i_def), intent(in) :: map_wt_stencil_size integer(kind=i_def), dimension(ndf_wt,map_wt_stencil_size), intent(in) :: map_wt_stencil integer(kind=i_def), intent(in) :: map_pid_stencil_size integer(kind=i_def), dimension(ndf_pid,map_pid_stencil_size), intent(in) :: map_pid_stencil - integer(kind=i_def), dimension(ndf_w1), intent(in) :: map_w1 - integer(kind=i_def), dimension(ndf_w2), intent(in) :: map_w2 - integer(kind=i_def), dimension(ndf_wt), intent(in) :: map_wt - integer(kind=i_def), dimension(ndf_pid), intent(in) :: map_pid - - real(kind=r_def), dimension(undf_w2), intent(inout) :: u_inc - real(kind=r_def), dimension(undf_w2), intent(in) :: u_n - real(kind=r_def), dimension(undf_wt), intent(in) :: velocity_w2v - real(kind=r_def), dimension(undf_wt), intent(in) :: vel_w2v_inc - real(kind=r_def), dimension(undf_w2), intent(in) :: dtrdz_fd2 - real(kind=r_def), dimension(undf_w1), intent(in) :: height_w1 - real(kind=r_def), dimension(undf_w2), intent(in) :: height_w2 - real(kind=r_def), dimension(undf_w2), intent(in) :: wetrho_in_w2 - real(kind=r_def), dimension(undf_pid), intent(in) :: panel_id - real(kind=r_def), intent(in) :: planet_radius - real(kind=r_def), intent(in) :: leonard_kl - real(kind=r_def), intent(in) :: dt + integer(kind=i_def), dimension(ndf_w1), intent(in) :: map_w1 + integer(kind=i_def), dimension(ndf_w2), intent(in) :: map_w2 + integer(kind=i_def), dimension(ndf_wt), intent(in) :: map_wt + integer(kind=i_def), dimension(ndf_pid), intent(in) :: map_pid + integer(kind=i_def), dimension(ndf_w3_2d), intent(in) :: map_w3_2d + + real(kind=r_def), dimension(undf_w2), intent(inout) :: u_inc + real(kind=r_def), dimension(undf_w2), intent(in) :: u_n + real(kind=r_def), dimension(undf_wt), intent(in) :: velocity_w2v + real(kind=r_def), dimension(undf_wt), intent(in) :: vel_w2v_inc + real(kind=r_def), dimension(undf_w2), intent(in) :: dtrdz_fd2 + real(kind=r_def), dimension(undf_w1), intent(in) :: height_w1 + real(kind=r_def), dimension(undf_w2), intent(in) :: height_w2 + real(kind=r_def), dimension(undf_w2), intent(in) :: wetrho_in_w2 + real(kind=r_def), dimension(undf_pid), intent(in) :: panel_id + real(kind=r_def), intent(in) :: planet_radius + real(kind=r_def), intent(in) :: leonard_kl + real(kind=r_def), intent(in) :: dt + + integer(kind=i_def), dimension(undf_w3_2d), intent(in) :: face_selector_ew + integer(kind=i_def), dimension(undf_w3_2d), intent(in) :: face_selector_ns ! Internal variables - integer(kind=i_def) :: k, km, df + integer(kind=i_def) :: k, j, km, df integer(kind=i_def) :: df2, df2p1, df2p2, df2m1, df2m2, dfp2x2 ! interpolation weights @@ -232,206 +252,205 @@ subroutine leonard_term_u_code( nlayers, & end do ! Loop over horizontal faces of the cell - do df = 1,4 - - if (u_inc(map_w2(df)+1) == 0.0_r_def) then - - true_wt_stencil(:,1:map_wt_stencil_size) = map_wt_stencil - true_w2_stencil(:,1:map_w2_stencil_size) = rot_w2_stencil - vec_dir(:,1:map_w2_stencil_size) = rot_vec_dir - - ! To generate a 3x3 stencil, we need to pad the corners with something. We do this by - ! duplicating the 'missing' point from the adjacent point, such that when we're calculating - ! a difference, we find the difference across the corner. i.e. when computing the 'x' fluxes - ! (dofs 1 & 3), the stencil maps at the 4 corners should look like - ! b-lft b-rgt t-rgt t-lft - ! 8 7 6 8 7 6 8 7 7 8 8 7 - ! 2 1 5 2 1 5 2 1 6 2 1 6 - ! 3 3 4 3 4 4 3 4 5 3 4 5 - ! and for the 'y' fluxes (dofs 2 & 4), the corner maps should be - ! 8 7 6 8 7 6 8 7 6 2 8 7 - ! 2 1 5 2 1 5 2 1 6 2 1 6 - ! 2 3 4 3 4 5 3 4 5 3 4 5 - if (map_w2_stencil_size == 8_i_def) then - if (int(panel_id(map_pid_stencil(1, 2))) /= int(panel_id(map_pid_stencil(1, 3))) .and. & - int(panel_id(map_pid_stencil(1, 6))) == int(panel_id(map_pid_stencil(1, 7))) ) then - ! bottom left corner - if (df == 1 .or. df == 3) then - true_wt_stencil(1,4:9) = map_wt_stencil(1,3:8) - true_w2_stencil(:,4:9) = rot_w2_stencil(:,3:8) - vec_dir(:,4:9) = rot_vec_dir(:,3:8) - else if (df == 2 .or. df == 4) then - true_wt_stencil(1,3:9) = map_wt_stencil(1,2:8) - true_w2_stencil(:,3:9) = rot_w2_stencil(:,2:8) - vec_dir(:,3:9) = rot_vec_dir(:,2:8) - end if - else if (int(panel_id(map_pid_stencil(1, 4))) /= int(panel_id(map_pid_stencil(1, 5))) .and. & - int(panel_id(map_pid_stencil(1, 2))) == int(panel_id(map_pid_stencil(1, 8))) ) then - ! bottom right corner - if (df == 1 .or. df == 3) then - true_wt_stencil(1,5:9) = map_wt_stencil(1,4:8) - true_w2_stencil(:,5:9) = rot_w2_stencil(:,4:8) - vec_dir(:,5:9) = rot_vec_dir(:,4:8) - else if (df == 2 .or. df == 4) then - true_wt_stencil(1,6:9) = map_wt_stencil(1,5:8) - true_w2_stencil(:,6:9) = rot_w2_stencil(:,5:8) - vec_dir(:,6:9) = rot_vec_dir(:,5:8) - end if - else if (int(panel_id(map_pid_stencil(1, 6))) /= int(panel_id(map_pid_stencil(1, 7))) .and. & - int(panel_id(map_pid_stencil(1, 3))) == int(panel_id(map_pid_stencil(1, 4))) ) then - ! top right corner - if (df == 1 .or. df == 3) then - true_wt_stencil(1,8:9) = map_wt_stencil(1,7:8) - true_w2_stencil(:,8:9) = rot_w2_stencil(:,7:8) - vec_dir(:,8:9) = rot_vec_dir(:,7:8) - else if (df == 2 .or. df == 4) then - true_wt_stencil(1,7:9) = map_wt_stencil(1,6:8) - true_w2_stencil(:,7:9) = rot_w2_stencil(:,6:8) - vec_dir(:,7:9) = rot_vec_dir(:,6:8) - end if - else if (int(panel_id(map_pid_stencil(1, 2))) /= int(panel_id(map_pid_stencil(1, 8))) ) then - ! top left corner - if (df == 1 .or. df == 3) then - true_wt_stencil(1,9) = map_wt_stencil(1,8) - true_w2_stencil(:,9) = rot_w2_stencil(:,8) - vec_dir(:,9) = rot_vec_dir(:,8) - else if (df == 2 .or. df == 4) then - true_wt_stencil(1,9) = map_wt_stencil(1,2) - true_w2_stencil(:,9) = rot_w2_stencil(:,2) - vec_dir(:,9) = rot_vec_dir(:,2) - end if + do j = 1, face_selector_ew(map_w3_2d(1)) + face_selector_ns(map_w3_2d(1)) + df = j + if (j == 3 .and. face_selector_ns(map_w3_2d(1)) == 2 & + .and. face_selector_ew(map_w3_2d(1)) == 1) df = N + + true_wt_stencil(:,1:map_wt_stencil_size) = map_wt_stencil + true_w2_stencil(:,1:map_w2_stencil_size) = rot_w2_stencil + vec_dir(:,1:map_w2_stencil_size) = rot_vec_dir + + ! To generate a 3x3 stencil, we need to pad the corners with something. We do this by + ! duplicating the 'missing' point from the adjacent point, such that when we're calculating + ! a difference, we find the difference across the corner. i.e. when computing the 'x' fluxes + ! (dofs 1 & 3), the stencil maps at the 4 corners should look like + ! b-lft b-rgt t-rgt t-lft + ! 8 7 6 8 7 6 8 7 7 8 8 7 + ! 2 1 5 2 1 5 2 1 6 2 1 6 + ! 3 3 4 3 4 4 3 4 5 3 4 5 + ! and for the 'y' fluxes (dofs 2 & 4), the corner maps should be + ! 8 7 6 8 7 6 8 7 6 2 8 7 + ! 2 1 5 2 1 5 2 1 6 2 1 6 + ! 2 3 4 3 4 5 3 4 5 3 4 5 + if (map_w2_stencil_size == 8_i_def) then + if (int(panel_id(map_pid_stencil(1, 2))) /= int(panel_id(map_pid_stencil(1, 3))) .and. & + int(panel_id(map_pid_stencil(1, 6))) == int(panel_id(map_pid_stencil(1, 7))) ) then + ! bottom left corner + if (df == 1 .or. df == 3) then + true_wt_stencil(1,4:9) = map_wt_stencil(1,3:8) + true_w2_stencil(:,4:9) = rot_w2_stencil(:,3:8) + vec_dir(:,4:9) = rot_vec_dir(:,3:8) + else if (df == 2 .or. df == 4) then + true_wt_stencil(1,3:9) = map_wt_stencil(1,2:8) + true_w2_stencil(:,3:9) = rot_w2_stencil(:,2:8) + vec_dir(:,3:9) = rot_vec_dir(:,2:8) + end if + else if (int(panel_id(map_pid_stencil(1, 4))) /= int(panel_id(map_pid_stencil(1, 5))) .and. & + int(panel_id(map_pid_stencil(1, 2))) == int(panel_id(map_pid_stencil(1, 8))) ) then + ! bottom right corner + if (df == 1 .or. df == 3) then + true_wt_stencil(1,5:9) = map_wt_stencil(1,4:8) + true_w2_stencil(:,5:9) = rot_w2_stencil(:,4:8) + vec_dir(:,5:9) = rot_vec_dir(:,4:8) + else if (df == 2 .or. df == 4) then + true_wt_stencil(1,6:9) = map_wt_stencil(1,5:8) + true_w2_stencil(:,6:9) = rot_w2_stencil(:,5:8) + vec_dir(:,6:9) = rot_vec_dir(:,5:8) + end if + else if (int(panel_id(map_pid_stencil(1, 6))) /= int(panel_id(map_pid_stencil(1, 7))) .and. & + int(panel_id(map_pid_stencil(1, 3))) == int(panel_id(map_pid_stencil(1, 4))) ) then + ! top right corner + if (df == 1 .or. df == 3) then + true_wt_stencil(1,8:9) = map_wt_stencil(1,7:8) + true_w2_stencil(:,8:9) = rot_w2_stencil(:,7:8) + vec_dir(:,8:9) = rot_vec_dir(:,7:8) + else if (df == 2 .or. df == 4) then + true_wt_stencil(1,7:9) = map_wt_stencil(1,6:8) + true_w2_stencil(:,7:9) = rot_w2_stencil(:,6:8) + vec_dir(:,7:9) = rot_vec_dir(:,6:8) + end if + else if (int(panel_id(map_pid_stencil(1, 2))) /= int(panel_id(map_pid_stencil(1, 8))) ) then + ! top left corner + if (df == 1 .or. df == 3) then + true_wt_stencil(1,9) = map_wt_stencil(1,8) + true_w2_stencil(:,9) = rot_w2_stencil(:,8) + vec_dir(:,9) = rot_vec_dir(:,8) + else if (df == 2 .or. df == 4) then + true_wt_stencil(1,9) = map_wt_stencil(1,2) + true_w2_stencil(:,9) = rot_w2_stencil(:,2) + vec_dir(:,9) = rot_vec_dir(:,2) end if end if + end if - ! Calculate indices for finite differences between faces of cells - df2 = 2_i_def * df - df2m1 = df2 - 1_i_def - df2m2 = df2 - 2_i_def - df2p1 = df2 + 1_i_def - df2p2 = df2 + 2_i_def - dfp2x2 = 2_i_def * (df + 2_i_def) - - ! If index is less than 2 add 8 - if (df2m1 < 2_i_def) then - df2m1 = df2m1 + 8_i_def - end if - if (df2m2 < 2_i_def) then - df2m2 = df2m2 + 8_i_def - end if - ! If index is greater than 9 subtract 8 - if (df2p2 > 9_i_def) then - df2p2 = df2p2 - 8_i_def - end if - if (dfp2x2 > 9_i_def) then - dfp2x2 = dfp2x2 - 8_i_def - end if - - ! Calculate kl at FD1 points, - ! accounting for stability limit - do k = 1, bl_levels - kl_fd1(k,df) = MIN( leonard_kl, & - 6.0_r_def * ( height_w2(map_w2(df) + k) - & - height_w2(map_w2(df) + k-1) ) & - / (dt * MAX( & - ! Difference normal to face - ! (point lies exactly between these): - ABS( velocity_w2v(true_wt_stencil(1,df2) + k) - & - velocity_w2v(true_wt_stencil(1,1) + k) ), & - ! Terms from gradient parallel to face... - ! from difference to the left (when looking at face): - ABS( velocity_w2v(true_wt_stencil(1,1) + k) - & - velocity_w2v(true_wt_stencil(1,df2m2) + k) ), & - ABS( velocity_w2v(true_wt_stencil(1,df2) + k) - & - velocity_w2v(true_wt_stencil(1,df2m1) + k) ), & - ! from difference to the right (when looking at face): - ABS( velocity_w2v(true_wt_stencil(1,df2p2) + k) - & - velocity_w2v(true_wt_stencil(1,1) + k) ), & - ABS( velocity_w2v(true_wt_stencil(1,df2p1) + k) - & - velocity_w2v(true_wt_stencil(1,df2) + k) ), & - EPSILON( leonard_kl ) & - ) ) ) - end do - - ! Calculate rho * r^2 at FD1 - do k = 1, bl_levels - km = k - 1 + ! Calculate indices for finite differences between faces of cells + df2 = 2_i_def * df + df2m1 = df2 - 1_i_def + df2m2 = df2 - 2_i_def + df2p1 = df2 + 1_i_def + df2p2 = df2 + 2_i_def + dfp2x2 = 2_i_def * (df + 2_i_def) + + ! If index is less than 2 add 8 + if (df2m1 < 2_i_def) then + df2m1 = df2m1 + 8_i_def + end if + if (df2m2 < 2_i_def) then + df2m2 = df2m2 + 8_i_def + end if + ! If index is greater than 9 subtract 8 + if (df2p2 > 9_i_def) then + df2p2 = df2p2 - 8_i_def + end if + if (dfp2x2 > 9_i_def) then + dfp2x2 = dfp2x2 - 8_i_def + end if - ! Vertical interpolation weights: - weight_pl = (height_w1(map_w1(df) + k) - height_w2(map_w2(df) + km)) / & + ! Calculate kl at FD1 points, + ! accounting for stability limit + do k = 1, bl_levels + kl_fd1(k,df) = MIN( leonard_kl, & + 6.0_r_def * ( height_w2(map_w2(df) + k) - & + height_w2(map_w2(df) + k-1) ) & + / (dt * MAX( & + ! Difference normal to face + ! (point lies exactly between these): + ABS( velocity_w2v(true_wt_stencil(1,df2) + k) - & + velocity_w2v(true_wt_stencil(1,1) + k) ), & + ! Terms from gradient parallel to face... + ! from difference to the left (when looking at face): + ABS( velocity_w2v(true_wt_stencil(1,1) + k) - & + velocity_w2v(true_wt_stencil(1,df2m2) + k) ), & + ABS( velocity_w2v(true_wt_stencil(1,df2) + k) - & + velocity_w2v(true_wt_stencil(1,df2m1) + k) ), & + ! from difference to the right (when looking at face): + ABS( velocity_w2v(true_wt_stencil(1,df2p2) + k) - & + velocity_w2v(true_wt_stencil(1,1) + k) ), & + ABS( velocity_w2v(true_wt_stencil(1,df2p1) + k) - & + velocity_w2v(true_wt_stencil(1,df2) + k) ), & + EPSILON( leonard_kl ) & + ) ) ) + end do + + ! Calculate rho * r^2 at FD1 + do k = 1, bl_levels + km = k - 1 + + ! Vertical interpolation weights: + weight_pl = (height_w1(map_w1(df) + k) - height_w2(map_w2(df) + km)) / & + (height_w2(map_w2(df) + k) - height_w2(map_w2(df) + km)) + weight_min = (height_w2(map_w2(df) + k) - height_w1(map_w1(df) + k)) / & (height_w2(map_w2(df) + k) - height_w2(map_w2(df) + km)) - weight_min = (height_w2(map_w2(df) + k) - height_w1(map_w1(df) + k)) / & - (height_w2(map_w2(df) + k) - height_w2(map_w2(df) + km)) - - ! Vertical interpolation of wetrho_in_w2 from w2 to FD1 points - rho_fd1 = ( weight_min * wetrho_in_w2(map_w2(df) + km) ) + & - ( weight_pl * wetrho_in_w2(map_w2(df) + k) ) - - ! Calculate rho * r^2 at FD1 - rho_rsq_fd1(k,df) = rho_fd1 * & - ( height_w1(map_w1(df) + k) + planet_radius ) * & - ( height_w1(map_w1(df) + k) + planet_radius ) - - end do - - ! Calculate timestep / ( dz * rho * r^2 ) at w2 - do k = 0, bl_levels-1 - dtrdzrho_fd2(k,df) = dtrdz_fd2(map_w2(df) + k) / & - wetrho_in_w2(map_w2(df) + k) - end do - - ! Calculate vertical flux at FD1 points: - ! Leonard term sub-grid vertical flux is proportional to the - ! product of the grid-scale horizontal differences in u and w. - - ! Set flux and increment to zero at k=0 - k = 0 - flux(k,df) = 0.0_r_def - - do k = 1, bl_levels - km = k - 1 - flux(k,df) = ( kl_fd1(k,df) / 12.0_r_def ) & - ! 8 terms contribute to each direction, so scale by 1/8 - * ( 1.0_r_def / 8.0_r_def ) * ( & - ! Terms from gradient normal to face... - ! ( point is half-way between 2 w-points, - ! so only one dw_x contributes ): - 2.0_r_def * ( ( vec_dir(df,dfp2x2)*u_n(true_w2_stencil(df,dfp2x2) + km) - & - vec_dir(df,df2)*u_n(true_w2_stencil(df,df2) + km) ) & - + ( vec_dir(df,dfp2x2)*u_n(true_w2_stencil(df,dfp2x2) + k) - & - vec_dir(df,df2)*u_n(true_w2_stencil(df,df2) + k) ) ) & - * ( velocity_w2v(true_wt_stencil(1,1) + k) - & - velocity_w2v(true_wt_stencil(1,df2) + k) ) & - ! Terms from gradient parallel to face... - ! from difference to the left (when looking at face): - + ( ( u_n(true_w2_stencil(df,1) + km) - & - vec_dir(df,df2m2)*u_n(true_w2_stencil(df,df2m2) + km) ) & - + ( u_n(true_w2_stencil(df,1) + k) - & - vec_dir(df,df2m2)*u_n(true_w2_stencil(df,df2m2) + k) ) ) & - * ( ( velocity_w2v(true_wt_stencil(1,1) + k) - & - velocity_w2v(true_wt_stencil(1,df2m2) + k) ) & - + ( velocity_w2v(true_wt_stencil(1,df2) + k) - & - velocity_w2v(true_wt_stencil(1,df2m1) + k) ) ) & - ! from difference to the right (when looking at face): - + ( ( vec_dir(df,df2p2)*u_n(true_w2_stencil(df,df2p2) + km) - & - u_n(true_w2_stencil(df,1) + km) ) & - + ( vec_dir(df,df2p2)*u_n(true_w2_stencil(df,df2p2) + k) - & - u_n(true_w2_stencil(df,1) + k) ) ) & - * ( ( velocity_w2v(true_wt_stencil(1,df2p2) + k) - & - velocity_w2v(true_wt_stencil(1,1) + k) ) & - + ( velocity_w2v(true_wt_stencil(1,df2p1) + k) - & - velocity_w2v(true_wt_stencil(1,df2) + k) ) ) & - ) & - * rho_rsq_fd1(k,df) - ! Flux now includes density * r^2 factor - end do - - ! Difference the flux in the vertical to get increment on w2 - do k = 0, bl_levels-1 - u_inc(map_w2(df) + k) = -(flux(k+1,df) - flux(k,df)) & - * dtrdzrho_fd2(k,df) - end do - end if + ! Vertical interpolation of wetrho_in_w2 from w2 to FD1 points + rho_fd1 = ( weight_min * wetrho_in_w2(map_w2(df) + km) ) + & + ( weight_pl * wetrho_in_w2(map_w2(df) + k) ) + + ! Calculate rho * r^2 at FD1 + rho_rsq_fd1(k,df) = rho_fd1 * & + ( height_w1(map_w1(df) + k) + planet_radius ) * & + ( height_w1(map_w1(df) + k) + planet_radius ) + + end do + + ! Calculate timestep / ( dz * rho * r^2 ) at w2 + do k = 0, bl_levels-1 + dtrdzrho_fd2(k,df) = dtrdz_fd2(map_w2(df) + k) / & + wetrho_in_w2(map_w2(df) + k) + end do + + ! Calculate vertical flux at FD1 points: + ! Leonard term sub-grid vertical flux is proportional to the + ! product of the grid-scale horizontal differences in u and w. + + ! Set flux and increment to zero at k=0 + k = 0 + flux(k,df) = 0.0_r_def + + do k = 1, bl_levels + km = k - 1 + flux(k,df) = ( kl_fd1(k,df) / 12.0_r_def ) & + ! 8 terms contribute to each direction, so scale by 1/8 + * ( 1.0_r_def / 8.0_r_def ) * ( & + ! Terms from gradient normal to face... + ! ( point is half-way between 2 w-points, + ! so only one dw_x contributes ): + 2.0_r_def * ( ( vec_dir(df,dfp2x2)*u_n(true_w2_stencil(df,dfp2x2) + km) - & + vec_dir(df,df2)*u_n(true_w2_stencil(df,df2) + km) ) & + + ( vec_dir(df,dfp2x2)*u_n(true_w2_stencil(df,dfp2x2) + k) - & + vec_dir(df,df2)*u_n(true_w2_stencil(df,df2) + k) ) ) & + * ( velocity_w2v(true_wt_stencil(1,1) + k) - & + velocity_w2v(true_wt_stencil(1,df2) + k) ) & + ! Terms from gradient parallel to face... + ! from difference to the left (when looking at face): + + ( ( u_n(true_w2_stencil(df,1) + km) - & + vec_dir(df,df2m2)*u_n(true_w2_stencil(df,df2m2) + km) ) & + + ( u_n(true_w2_stencil(df,1) + k) - & + vec_dir(df,df2m2)*u_n(true_w2_stencil(df,df2m2) + k) ) ) & + * ( ( velocity_w2v(true_wt_stencil(1,1) + k) - & + velocity_w2v(true_wt_stencil(1,df2m2) + k) ) & + + ( velocity_w2v(true_wt_stencil(1,df2) + k) - & + velocity_w2v(true_wt_stencil(1,df2m1) + k) ) ) & + ! from difference to the right (when looking at face): + + ( ( vec_dir(df,df2p2)*u_n(true_w2_stencil(df,df2p2) + km) - & + u_n(true_w2_stencil(df,1) + km) ) & + + ( vec_dir(df,df2p2)*u_n(true_w2_stencil(df,df2p2) + k) - & + u_n(true_w2_stencil(df,1) + k) ) ) & + * ( ( velocity_w2v(true_wt_stencil(1,df2p2) + k) - & + velocity_w2v(true_wt_stencil(1,1) + k) ) & + + ( velocity_w2v(true_wt_stencil(1,df2p1) + k) - & + velocity_w2v(true_wt_stencil(1,df2) + k) ) ) & + ) & + * rho_rsq_fd1(k,df) + ! Flux now includes density * r^2 factor + end do + + ! Difference the flux in the vertical to get increment on w2 + do k = 0, bl_levels-1 + u_inc(map_w2(df) + k) = -(flux(k+1,df) - flux(k,df)) & + * dtrdzrho_fd2(k,df) + end do end do diff --git a/science/gungho/source/kernel/diffusion/momentum_smagorinsky_kernel_mod.F90 b/science/gungho/source/kernel/diffusion/momentum_smagorinsky_kernel_mod.F90 index bb39a4bd8..674c880ee 100644 --- a/science/gungho/source/kernel/diffusion/momentum_smagorinsky_kernel_mod.F90 +++ b/science/gungho/source/kernel/diffusion/momentum_smagorinsky_kernel_mod.F90 @@ -13,14 +13,18 @@ !> module momentum_smagorinsky_kernel_mod - use argument_mod, only : arg_type, & - GH_FIELD, GH_REAL, & - GH_READ, GH_WRITE, & - STENCIL, CROSS, CELL_COLUMN, & - ANY_DISCONTINUOUS_SPACE_9 - use constants_mod, only : r_def, i_def - use fs_continuity_mod, only : W2, W1, Wtheta - use kernel_mod, only : kernel_type + use argument_mod, only : arg_type, & + GH_FIELD, GH_REAL, & + GH_READ, GH_WRITE, & + STENCIL, CROSS, CELL_COLUMN, & + ANY_DISCONTINUOUS_SPACE_9, & + ANY_DISCONTINUOUS_SPACE_3, & + ANY_DISCONTINUOUS_SPACE_2, & + GH_INTEGER + use constants_mod, only : r_def, i_def + use fs_continuity_mod, only : W2, W1, Wtheta + use kernel_mod, only : kernel_type + use reference_element_mod, only : N implicit none private @@ -33,16 +37,20 @@ module momentum_smagorinsky_kernel_mod !> type, public, extends(kernel_type) :: momentum_smagorinsky_kernel_type private - type(arg_type) :: meta_args(7) = (/ & - arg_type(GH_FIELD, GH_REAL, GH_WRITE, W2), & - arg_type(GH_FIELD, GH_REAL, GH_READ, W2, STENCIL(CROSS)), & - arg_type(GH_FIELD, GH_REAL, GH_READ, W2, STENCIL(CROSS)), & - arg_type(GH_FIELD, GH_REAL, GH_READ, W2), & - arg_type(GH_FIELD, GH_REAL, GH_READ, W1), & - arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta, STENCIL(CROSS)), & - arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_9, & - STENCIL(CROSS)) & - /) + type(arg_type) :: meta_args(9) = (/ & + arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_2), & + arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_2, & + STENCIL(CROSS)), & + arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_2, & + STENCIL(CROSS)), & + arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_2), & + arg_type(GH_FIELD, GH_REAL, GH_READ, W1), & + arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta, STENCIL(CROSS)), & + arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_9, & + STENCIL(CROSS)), & + arg_type(GH_FIELD, GH_INTEGER, GH_READ, ANY_DISCONTINUOUS_SPACE_3), & + arg_type(GH_FIELD, GH_INTEGER, GH_READ, ANY_DISCONTINUOUS_SPACE_3) & + /) integer :: operates_on = CELL_COLUMN contains procedure, nopass :: momentum_smagorinsky_code @@ -56,34 +64,39 @@ module momentum_smagorinsky_kernel_mod contains !> @brief Calculates diffusion increment for wind field using horizontal Smagorinsky diffusion -!! @param[in] nlayers Number of layers in the mesh -!! @param[in,out] u_inc Diffusion increment for wind field -!! @param[in] u_n Input wind field -!! @param[in] map_w2_stencil_size Number of cells in the stencil at the base of the column for w2 -!! @param[in] map_w2_stencil Array holding the stencil dofmap for the cell at the base of the column for w2 -!! @param[in] dx_at_w2 Grid length at the w2 dofs -!! @param[in] map_dx_stencil_size Number of cells in the stencil at the base of the column for w2 -!! @param[in] map_dx_stencil Array holding the stencil dofmap for the cell at the base of the column for w2 -!! @param[in] height_w2 Height of w3 space levels above surface at cell faces -!! @param[in] height_w1 Height of wtheta levels above surface at cell faces -!! @param[in] visc_m Diffusion coefficient for momentum on wtheta points -!! @param[in] map_wt_stencil_size Number of cells in the stencil at the base of the column for wt -!! @param[in] map_wt_stencil Array holding the stencil dofmap for the cell at the base of the column for wt -!! @param[in] panel_id The ID number of the current panel -!! @param[in] map_pid_stencil_size Size of the panel ID stencil -!! @param[in] map_pid_stencil Stencil map for the panel ID -!! @param[in] ndf_w2 Number of degrees of freedom per cell for wind space -!! @param[in] undf_w2 Number of unique degrees of freedom for wind space -!! @param[in] map_w2 Array holding the dofmap for the cell at the base of the column for wind space -!! @param[in] ndf_w1 Number of degrees of freedom per cell for w1 -!! @param[in] undf_w1 Number of unique degrees of freedom for w1 -!! @param[in] map_w1 Array holding the dofmap for the cell at the base of the column for w1 -!! @param[in] ndf_wt Number of degrees of freedom per cell for theta space -!! @param[in] undf_wt Number of unique degrees of freedom for theta space -!! @param[in] map_wt Array holding the dofmap for the cell at the base of the column for theta space -!! @param[in] ndf_pid Number of degrees of freedom per cell for pid space -!! @param[in] undf_pid Number of unique degrees of freedom for pid space -!! @param[in] map_pid Cell dofmap for pid space +!> @param[in] nlayers Number of layers in the mesh +!> @param[in,out] u_inc Diffusion increment for wind field +!> @param[in] u_n Input wind field +!> @param[in] map_w2_stencil_size Number of cells in the stencil at the base of the column for w2 +!> @param[in] map_w2_stencil Array holding the stencil dofmap for the cell at the base of the column for w2 +!> @param[in] dx_at_w2 Grid length at the w2 dofs +!> @param[in] map_dx_stencil_size Number of cells in the stencil at the base of the column for w2 +!> @param[in] map_dx_stencil Array holding the stencil dofmap for the cell at the base of the column for w2 +!> @param[in] height_w2 Height of w3 space levels above surface at cell faces +!> @param[in] height_w1 Height of wtheta levels above surface at cell faces +!> @param[in] visc_m Diffusion coefficient for momentum on wtheta points +!> @param[in] map_wt_stencil_size Number of cells in the stencil at the base of the column for wt +!> @param[in] map_wt_stencil Array holding the stencil dofmap for the cell at the base of the column for wt +!> @param[in] panel_id The ID number of the current panel +!> @param[in] map_pid_stencil_size Size of the panel ID stencil +!> @param[in] map_pid_stencil Stencil map for the panel ID +!> @param[in] face_selector_ew 2D field indicating which W/E faces to loop over in this column +!> @param[in] face_selector_ns 2D field indicating which N/S faces to loop over in this column +!> @param[in] ndf_w2 Number of degrees of freedom per cell for wind space +!> @param[in] undf_w2 Number of unique degrees of freedom for wind space +!> @param[in] map_w2 Array holding the dofmap for the cell at the base of the column for wind space +!> @param[in] ndf_w1 Number of degrees of freedom per cell for w1 +!> @param[in] undf_w1 Number of unique degrees of freedom for w1 +!> @param[in] map_w1 Array holding the dofmap for the cell at the base of the column for w1 +!> @param[in] ndf_wt Number of degrees of freedom per cell for theta space +!> @param[in] undf_wt Number of unique degrees of freedom for theta space +!> @param[in] map_wt Array holding the dofmap for the cell at the base of the column for theta space +!> @param[in] ndf_pid Number of degrees of freedom per cell for pid space +!> @param[in] undf_pid Number of unique degrees of freedom for pid space +!> @param[in] map_pid Cell dofmap for pid space +!> @param[in] ndf_w3_2d Num of DoFs for 2D W3 per cell +!> @param[in] undf_w3_2d Num of DoFs for this partition for 2D W3 +!> @param[in] map_w3_2d Map for 2D W3 subroutine momentum_smagorinsky_code( nlayers, & u_inc, & u_n, & @@ -96,10 +109,12 @@ subroutine momentum_smagorinsky_code( nlayers, & map_wt_stencil_size, map_wt_stencil, & panel_id, & map_pid_stencil_size, map_pid_stencil, & + face_selector_ew, face_selector_ns, & ndf_w2, undf_w2, map_w2, & ndf_w1, undf_w1, map_w1, & ndf_wt, undf_wt, map_wt, & - ndf_pid, undf_pid, map_pid & + ndf_pid, undf_pid, map_pid, & + ndf_w3_2d, undf_w3_2d, map_w3_2d & ) implicit none @@ -108,15 +123,17 @@ subroutine momentum_smagorinsky_code( nlayers, & integer(kind=i_def), intent(in) :: nlayers integer(kind=i_def), intent(in) :: ndf_w2, ndf_w1, ndf_wt, ndf_pid integer(kind=i_def), intent(in) :: undf_w2, undf_w1, undf_wt, undf_pid + integer(kind=i_def), intent(in) :: ndf_w3_2d, undf_w3_2d integer(kind=i_def), intent(in) :: map_w2_stencil_size, map_dx_stencil_size, map_wt_stencil_size, map_pid_stencil_size integer(kind=i_def), dimension(ndf_w2,map_w2_stencil_size), intent(in) :: map_w2_stencil integer(kind=i_def), dimension(ndf_w2,map_dx_stencil_size), intent(in) :: map_dx_stencil integer(kind=i_def), dimension(ndf_wt,map_wt_stencil_size), intent(in) :: map_wt_stencil integer(kind=i_def), dimension(ndf_pid,map_pid_stencil_size), intent(in) :: map_pid_stencil - integer(kind=i_def), dimension(ndf_w2), intent(in) :: map_w2 - integer(kind=i_def), dimension(ndf_w1), intent(in) :: map_w1 - integer(kind=i_def), dimension(ndf_wt), intent(in) :: map_wt - integer(kind=i_def), dimension(ndf_pid), intent(in) :: map_pid + integer(kind=i_def), dimension(ndf_w2), intent(in) :: map_w2 + integer(kind=i_def), dimension(ndf_w1), intent(in) :: map_w1 + integer(kind=i_def), dimension(ndf_wt), intent(in) :: map_wt + integer(kind=i_def), dimension(ndf_pid), intent(in) :: map_pid + integer(kind=i_def), dimension(ndf_w3_2d), intent(in) :: map_w3_2d real(kind=r_def), dimension(undf_w2), intent(inout) :: u_inc real(kind=r_def), dimension(undf_w2), intent(in) :: u_n, dx_at_w2 @@ -125,8 +142,11 @@ subroutine momentum_smagorinsky_code( nlayers, & real(kind=r_def), dimension(undf_wt), intent(in) :: visc_m real(kind=r_def), dimension(undf_pid), intent(in) :: panel_id + integer(kind=i_def), dimension(undf_w3_2d), intent(in) :: face_selector_ew + integer(kind=i_def), dimension(undf_w3_2d), intent(in) :: face_selector_ns + ! Internal variables - integer(kind=i_def) :: k, kp, df + integer(kind=i_def) :: k, kp, df, j real(kind=r_def) :: d2dx, d2dy real(kind=r_def), dimension(1:nlayers-1) :: idx2, idy2 real(kind=r_def), dimension(0:nlayers-1,4) :: idx2_w2, idy2_w2 @@ -227,64 +247,63 @@ subroutine momentum_smagorinsky_code( nlayers, & end do ! Loop over horizontal faces of the cell for horizontal diffusion - do df = 1,4 - - ! Only calculate face if it's not already been done - if (u_inc(map_w2(df)+1) == 0.0_r_def) then - - ! Compute horizontal grid spacing - if (df == 1 .or. df == 3) then - ! For Dofs 1 & 3, dx is given by the input field, - ! whilst dy needs to be computed from the 4 neighbouring dy points - do k = 0, nlayers - 1 - idx2_w2(k,df) = 1.0_r_def/(dx_at_w2(true_stencil_map(df,1)+k))**2 - idy2_w2(k,df) = (4.0_r_def/(dx_at_w2(true_stencil_map(2,1)+k)+ & - dx_at_w2(true_stencil_map(4,1)+k)+ & - dx_at_w2(true_stencil_map(2,df+1)+k)+ & - dx_at_w2(true_stencil_map(4,df+1)+k)))**2 - end do - else - ! For Dofs 2 & 4, dy is given by the input field, - ! whilst dx needs to be computed from the 4 neighbouring dx points - do k = 0, nlayers - 1 - idx2_w2(k,df) = (4.0_r_def/(dx_at_w2(true_stencil_map(1,1)+k)+ & - dx_at_w2(true_stencil_map(3,1)+k)+ & - dx_at_w2(true_stencil_map(1,df+1)+k)+ & - dx_at_w2(true_stencil_map(3,df+1)+k)))**2 - idy2_w2(k,df) = 1.0_r_def/(dx_at_w2(true_stencil_map(df,1)+k))**2 - end do - end if - - ! Horizontal interpolation of visc_m to cell face - do k = 0, nlayers - visc_m_w2(k,df) = ( visc_m(map_wt_stencil(1,1) + k) + visc_m(map_wt_stencil(1,1+df) + k) ) / 2.0_r_def + do j = 1, face_selector_ew(map_w3_2d(1)) + face_selector_ns(map_w3_2d(1)) + df = j + if (j == 3 .and. face_selector_ns(map_w3_2d(1)) == 2 & + .and. face_selector_ew(map_w3_2d(1)) == 1) df = N + + ! Compute horizontal grid spacing + if (df == 1 .or. df == 3) then + ! For Dofs 1 & 3, dx is given by the input field, + ! whilst dy needs to be computed from the 4 neighbouring dy points + do k = 0, nlayers - 1 + idx2_w2(k,df) = 1.0_r_def/(dx_at_w2(true_stencil_map(df,1)+k))**2 + idy2_w2(k,df) = (4.0_r_def/(dx_at_w2(true_stencil_map(2,1)+k)+ & + dx_at_w2(true_stencil_map(4,1)+k)+ & + dx_at_w2(true_stencil_map(2,df+1)+k)+ & + dx_at_w2(true_stencil_map(4,df+1)+k)))**2 end do - - ! Horizontal velocity diffusion + else + ! For Dofs 2 & 4, dy is given by the input field, + ! whilst dx needs to be computed from the 4 neighbouring dx points do k = 0, nlayers - 1 - kp = k + 1 + idx2_w2(k,df) = (4.0_r_def/(dx_at_w2(true_stencil_map(1,1)+k)+ & + dx_at_w2(true_stencil_map(3,1)+k)+ & + dx_at_w2(true_stencil_map(1,df+1)+k)+ & + dx_at_w2(true_stencil_map(3,df+1)+k)))**2 + idy2_w2(k,df) = 1.0_r_def/(dx_at_w2(true_stencil_map(df,1)+k))**2 + end do + end if - ! Vertical interpolation weights: - weight_pl = (height_w2(map_w2(df) + k) - height_w1(map_w1(df) + k)) / & + ! Horizontal interpolation of visc_m to cell face + do k = 0, nlayers + visc_m_w2(k,df) = ( visc_m(map_wt_stencil(1,1) + k) + visc_m(map_wt_stencil(1,1+df) + k) ) / 2.0_r_def + end do + + ! Horizontal velocity diffusion + do k = 0, nlayers - 1 + kp = k + 1 + + ! Vertical interpolation weights: + weight_pl = (height_w2(map_w2(df) + k) - height_w1(map_w1(df) + k)) / & + (height_w1(map_w1(df) + kp) - height_w1(map_w1(df) + k)) + weight_min = (height_w1(map_w1(df) + kp) - height_w2(map_w2(df) + k)) / & (height_w1(map_w1(df) + kp) - height_w1(map_w1(df) + k)) - weight_min = (height_w1(map_w1(df) + kp) - height_w2(map_w2(df) + k)) /& - (height_w1(map_w1(df) + kp) - height_w1(map_w1(df) + k)) - ! Vertical interpolation of visc_m from wt to w3 levels - visc_m_w2_w3 = ( weight_min * visc_m_w2(k,df) ) + ( weight_pl * visc_m_w2(kp,df) ) + ! Vertical interpolation of visc_m from wt to w3 levels + visc_m_w2_w3 = ( weight_min * visc_m_w2(k,df) ) + ( weight_pl * visc_m_w2(kp,df) ) - ! horizontal diffusion: - d2dx = (vec_dir(df,2)*u_n(true_stencil_map(df,2) + k) & - - 2.0_r_def*u_n(true_stencil_map(df,1) + k) & - + vec_dir(df,4)*u_n(true_stencil_map(df,4) + k) ) * idx2_w2(k,df) - d2dy = (vec_dir(df,3)*u_n(true_stencil_map(df,3) + k) & - - 2.0_r_def*u_n(true_stencil_map(df,1) + k) & - + vec_dir(df,5)*u_n(true_stencil_map(df,5) + k) ) * idy2_w2(k,df) - u_inc(map_w2(df) + k) = visc_m_w2_w3 * (d2dx + d2dy) + ! horizontal diffusion: + d2dx = (vec_dir(df,2)*u_n(true_stencil_map(df,2) + k) & + - 2.0_r_def*u_n(true_stencil_map(df,1) + k) & + + vec_dir(df,4)*u_n(true_stencil_map(df,4) + k) ) * idx2_w2(k,df) + d2dy = (vec_dir(df,3)*u_n(true_stencil_map(df,3) + k) & + - 2.0_r_def*u_n(true_stencil_map(df,1) + k) & + + vec_dir(df,5)*u_n(true_stencil_map(df,5) + k) ) * idy2_w2(k,df) + u_inc(map_w2(df) + k) = visc_m_w2_w3 * (d2dx + d2dy) - end do + end do - end if end do diff --git a/science/gungho/source/kernel/diffusion/momentum_viscosity_kernel_mod.F90 b/science/gungho/source/kernel/diffusion/momentum_viscosity_kernel_mod.F90 index 40d8056ea..0a2a808b7 100644 --- a/science/gungho/source/kernel/diffusion/momentum_viscosity_kernel_mod.F90 +++ b/science/gungho/source/kernel/diffusion/momentum_viscosity_kernel_mod.F90 @@ -8,14 +8,18 @@ !> module momentum_viscosity_kernel_mod - use argument_mod, only : arg_type, & - GH_FIELD, GH_REAL, & - GH_READ, GH_WRITE, & - GH_SCALAR, & - STENCIL, CROSS, CELL_COLUMN - use constants_mod, only : r_def, i_def - use fs_continuity_mod, only : W2 - use kernel_mod, only : kernel_type + use argument_mod, only : arg_type, & + GH_FIELD, GH_REAL, & + GH_READ, GH_WRITE, & + GH_SCALAR, & + GH_INTEGER, & + ANY_DISCONTINUOUS_SPACE_3, & + ANY_DISCONTINUOUS_SPACE_2, & + STENCIL, CROSS, CELL_COLUMN + use constants_mod, only : r_def, i_def + use fs_continuity_mod, only : W2 + use kernel_mod, only : kernel_type + use reference_element_mod, only : N implicit none @@ -29,12 +33,16 @@ module momentum_viscosity_kernel_mod !> type, public, extends(kernel_type) :: momentum_viscosity_kernel_type private - type(arg_type) :: meta_args(4) = (/ & - arg_type(GH_FIELD, GH_REAL, GH_WRITE, W2), & - arg_type(GH_FIELD, GH_REAL, GH_READ, W2, STENCIL(CROSS)), & - arg_type(GH_FIELD, GH_REAL, GH_READ, W2, STENCIL(CROSS)), & - arg_type(GH_SCALAR, GH_REAL, GH_READ) & - /) + type(arg_type) :: meta_args(6) = (/ & + arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_2), & + arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_2, & + STENCIL(CROSS)), & + arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_2, & + STENCIL(CROSS)), & + arg_type(GH_FIELD, GH_INTEGER, GH_READ, ANY_DISCONTINUOUS_SPACE_3), & + arg_type(GH_FIELD, GH_INTEGER, GH_READ, ANY_DISCONTINUOUS_SPACE_3), & + arg_type(GH_SCALAR, GH_REAL, GH_READ) & + /) integer :: operates_on = CELL_COLUMN contains procedure, nopass :: momentum_viscosity_code @@ -48,25 +56,32 @@ module momentum_viscosity_kernel_mod contains !> @brief The subroutine which is called directly by the Psy layer -!! @param[in] nlayers Number of layers in the mesh -!! @param[in,out] u_inc Diffusion increment for wind field -!! @param[in] u_n Input wind field -!! @param[in] map_w2_size Number of cells in the stencil at the base of the column for w2 -!! @param[in] map_w2 Array holding the stencil dofmap for the cell at the base of the column for w2 -!! @param[in] dx_at_w2 Grid length at the w2 dofs -!! @param[in] map_dx_stencil_size Number of cells in the stencil at the base of the column for w2 -!! @param[in] map_dx_stencil Array holding the stencil dofmap for the cell at the base of the column for w2 -!! @param[in] viscosity_mu Viscosity constant -!! @param[in] ndf_w2 Number of degrees of freedom per cell for wind space -!! @param[in] undf_w2 Number of unique degrees of freedom for wind_space -!! @param[in] cell_map_w2 Array holding the dofmap for the cell at the base of the column for w2 +!> @param[in] nlayers Number of layers in the mesh +!> @param[in,out] u_inc Diffusion increment for wind field +!> @param[in] u_n Input wind field +!> @param[in] map_w2_size Number of cells in the stencil at the base of the column for w2 +!> @param[in] map_w2 Array holding the stencil dofmap for the cell at the base of the column for w2 +!> @param[in] dx_at_w2 Grid length at the w2 dofs +!> @param[in] map_dx_stencil_size Number of cells in the stencil at the base of the column for w2 +!> @param[in] map_dx_stencil Array holding the stencil dofmap for the cell at the base of the column for w2 +!> @param[in] face_selector_ew 2D field indicating which W/E faces to loop over in this column +!> @param[in] face_selector_ns 2D field indicating which N/S faces to loop over in this column +!> @param[in] viscosity_mu Viscosity constant +!> @param[in] ndf_w2 Number of degrees of freedom per cell for wind space +!> @param[in] undf_w2 Number of unique degrees of freedom for wind_space +!> @param[in] cell_map_w2 Array holding the dofmap for the cell at the base of the column for w2 +!> @param[in] ndf_w3_2d Num of DoFs for 2D W3 per cell +!> @param[in] undf_w3_2d Num of DoFs for this partition for 2D W3 +!> @param[in] map_w3_2d Map for 2D W3 subroutine momentum_viscosity_code(nlayers, & u_inc, u_n, & map_w2_size, map_w2, & dx_at_w2, & map_dx_stencil_size, map_dx_stencil, & + face_selector_ew, face_selector_ns, & viscosity_mu, & - ndf_w2, undf_w2, cell_map_w2) + ndf_w2, undf_w2, cell_map_w2, & + ndf_w3_2d, undf_w3_2d, map_w3_2d ) implicit none @@ -74,17 +89,20 @@ subroutine momentum_viscosity_code(nlayers, & integer(kind=i_def), intent(in) :: nlayers integer(kind=i_def), intent(in) :: ndf_w2, undf_w2 integer(kind=i_def), intent(in) :: map_w2_size, map_dx_stencil_size - integer(kind=i_def), dimension(ndf_w2,map_w2_size), intent(in) :: map_w2 - integer(kind=i_def), dimension(ndf_w2,map_dx_stencil_size), intent(in) :: map_dx_stencil - integer(kind=i_def), dimension(ndf_w2), intent(in) :: cell_map_w2 - - real(kind=r_def), dimension(undf_w2), intent(inout) :: u_inc - real(kind=r_def), dimension(undf_w2), intent(in) :: u_n, dx_at_w2 - - real(kind=r_def), intent(in) :: viscosity_mu + integer(kind=i_def), intent(in) :: ndf_w3_2d, undf_w3_2d + integer(kind=i_def), intent(in) :: map_w2(ndf_w2,map_w2_size) + integer(kind=i_def), intent(in) :: map_dx_stencil(ndf_w2,map_dx_stencil_size) + integer(kind=i_def), intent(in) :: cell_map_w2(ndf_w2) + integer(kind=i_def), intent(in) :: map_w3_2d(ndf_w3_2d) + + real(kind=r_def), intent(inout) :: u_inc(undf_w2) + real(kind=r_def), intent(in) :: u_n(undf_w2), dx_at_w2(undf_w2) + real(kind=r_def), intent(in) :: viscosity_mu + integer(kind=i_def), intent(in) :: face_selector_ew(undf_w3_2d) + integer(kind=i_def), intent(in) :: face_selector_ns(undf_w3_2d) ! Internal variables - integer(kind=i_def) :: k, km, kp, df + integer(kind=i_def) :: k, j, km, kp, df real(kind=r_def) :: d2dx, d2dy, d2dz real(kind=r_def), dimension(0:nlayers-1) :: idx2, idy2, idz2 real(kind=r_def), dimension(0:nlayers-1,4) :: idx2_w2, idy2_w2, idz2_w2 @@ -102,63 +120,60 @@ subroutine momentum_viscosity_code(nlayers, & end if ! Loop over horizontal faces of the cell for horizontal diffusion - do df = 1,4 - - ! Only calculate face if it's not already been done - if (u_inc(cell_map_w2(df)+1) == 0.0_r_def) then - - ! Compute horizontal grid spacing - if (df == 1 .or. df == 3) then - ! For Dofs 1 & 3, dx is given by the input field, - ! whilst dy needs to be computed from the 4 neighbouring dy points - do k = 0, nlayers - 1 - idx2_w2(k,df) = 1.0_r_def/(dx_at_w2(map_dx_stencil(df,1)+k))**2 - idy2_w2(k,df) = (4.0_r_def/(dx_at_w2(map_dx_stencil(2,1)+k)+dx_at_w2(map_dx_stencil(4,1)+k)+ & - dx_at_w2(map_dx_stencil(2,df+1)+k)+dx_at_w2(map_dx_stencil(4,df+1)+k)))**2 - idz2_w2(k,df) = (2.0_r_def/(dx_at_w2(map_dx_stencil(5,1)+k)+dx_at_w2(map_dx_stencil(5,df+1)+k)))**2 - end do - else - ! For Dofs 2 & 4, dy is given by the input field, - ! whilst dx needs to be computed from the 4 neighbouring dx points - do k = 0, nlayers - 1 - idx2_w2(k,df) = (4.0_r_def/(dx_at_w2(map_dx_stencil(1,1)+k)+dx_at_w2(map_dx_stencil(3,1)+k)+ & - dx_at_w2(map_dx_stencil(1,df+1)+k)+dx_at_w2(map_dx_stencil(3,df+1)+k)))**2 - idy2_w2(k,df) = 1.0_r_def/(dx_at_w2(map_dx_stencil(df,1)+k))**2 - idz2_w2(k,df) = (2.0_r_def/(dx_at_w2(map_dx_stencil(5,1)+k)+dx_at_w2(map_dx_stencil(5,df+1)+k)))**2 - end do - end if - - ! Horizontal Velocity diffusion - k = 0 - km = k - kp = k+1 - - d2dx = (u_n(map_w2(df,2) + k) - 2.0_r_def*u_n(map_w2(df,1) + k) + u_n(map_w2(df,4) + k) )*idx2_w2(k,df) - d2dy = (u_n(map_w2(df,3) + k) - 2.0_r_def*u_n(map_w2(df,1) + k) + u_n(map_w2(df,5) + k) )*idy2_w2(k,df) - d2dz = (u_n(map_w2(df,1) + kp) - 2.0_r_def*u_n(map_w2(df,1) + k) + u_n(map_w2(df,1) + km))*idz2_w2(k,df) - u_inc(cell_map_w2(df)+k) = viscosity_mu*(d2dx + d2dy + d2dz) + do j = 1, face_selector_ew(map_w3_2d(1)) + face_selector_ns(map_w3_2d(1)) + df = j + if (j == 3 .and. face_selector_ns(map_w3_2d(1)) == 2 & + .and. face_selector_ew(map_w3_2d(1)) == 1) df = N + + ! Compute horizontal grid spacing + if (df == 1 .or. df == 3) then + ! For Dofs 1 & 3, dx is given by the input field, + ! whilst dy needs to be computed from the 4 neighbouring dy points + do k = 0, nlayers - 1 + idx2_w2(k,df) = 1.0_r_def/(dx_at_w2(map_dx_stencil(df,1)+k))**2 + idy2_w2(k,df) = (4.0_r_def/(dx_at_w2(map_dx_stencil(2,1)+k)+dx_at_w2(map_dx_stencil(4,1)+k)+ & + dx_at_w2(map_dx_stencil(2,df+1)+k)+dx_at_w2(map_dx_stencil(4,df+1)+k)))**2 + idz2_w2(k,df) = (2.0_r_def/(dx_at_w2(map_dx_stencil(5,1)+k)+dx_at_w2(map_dx_stencil(5,df+1)+k)))**2 + end do + else + ! For Dofs 2 & 4, dy is given by the input field, + ! whilst dx needs to be computed from the 4 neighbouring dx points + do k = 0, nlayers - 1 + idx2_w2(k,df) = (4.0_r_def/(dx_at_w2(map_dx_stencil(1,1)+k)+dx_at_w2(map_dx_stencil(3,1)+k)+ & + dx_at_w2(map_dx_stencil(1,df+1)+k)+dx_at_w2(map_dx_stencil(3,df+1)+k)))**2 + idy2_w2(k,df) = 1.0_r_def/(dx_at_w2(map_dx_stencil(df,1)+k))**2 + idz2_w2(k,df) = (2.0_r_def/(dx_at_w2(map_dx_stencil(5,1)+k)+dx_at_w2(map_dx_stencil(5,df+1)+k)))**2 + end do + end if - do k = 1, nlayers-2 - km = k - 1 - kp = k + 1 + ! Horizontal Velocity diffusion + k = 0 + km = k + kp = k+1 - d2dx = (u_n(map_w2(df,2) + k) - 2.0_r_def*u_n(map_w2(df,1) + k) + u_n(map_w2(df,4) + k) )*idx2_w2(k,df) - d2dy = (u_n(map_w2(df,3) + k) - 2.0_r_def*u_n(map_w2(df,1) + k) + u_n(map_w2(df,5) + k) )*idy2_w2(k,df) - d2dz = (u_n(map_w2(df,1) + kp) - 2.0_r_def*u_n(map_w2(df,1) + k) + u_n(map_w2(df,1) + km))*idz2_w2(k,df) - u_inc(cell_map_w2(df)+k) = viscosity_mu*(d2dx + d2dy + d2dz) - end do + d2dx = (u_n(map_w2(df,2) + k) - 2.0_r_def*u_n(map_w2(df,1) + k) + u_n(map_w2(df,4) + k) )*idx2_w2(k,df) + d2dy = (u_n(map_w2(df,3) + k) - 2.0_r_def*u_n(map_w2(df,1) + k) + u_n(map_w2(df,5) + k) )*idy2_w2(k,df) + d2dz = (u_n(map_w2(df,1) + kp) - 2.0_r_def*u_n(map_w2(df,1) + k) + u_n(map_w2(df,1) + km))*idz2_w2(k,df) + u_inc(cell_map_w2(df)+k) = viscosity_mu*(d2dx + d2dy + d2dz) - k = nlayers-1 + do k = 1, nlayers-2 km = k - 1 - kp = k + kp = k + 1 d2dx = (u_n(map_w2(df,2) + k) - 2.0_r_def*u_n(map_w2(df,1) + k) + u_n(map_w2(df,4) + k) )*idx2_w2(k,df) d2dy = (u_n(map_w2(df,3) + k) - 2.0_r_def*u_n(map_w2(df,1) + k) + u_n(map_w2(df,5) + k) )*idy2_w2(k,df) d2dz = (u_n(map_w2(df,1) + kp) - 2.0_r_def*u_n(map_w2(df,1) + k) + u_n(map_w2(df,1) + km))*idz2_w2(k,df) u_inc(cell_map_w2(df)+k) = viscosity_mu*(d2dx + d2dy + d2dz) + end do - end if + k = nlayers-1 + km = k - 1 + kp = k + d2dx = (u_n(map_w2(df,2) + k) - 2.0_r_def*u_n(map_w2(df,1) + k) + u_n(map_w2(df,4) + k) )*idx2_w2(k,df) + d2dy = (u_n(map_w2(df,3) + k) - 2.0_r_def*u_n(map_w2(df,1) + k) + u_n(map_w2(df,5) + k) )*idy2_w2(k,df) + d2dz = (u_n(map_w2(df,1) + kp) - 2.0_r_def*u_n(map_w2(df,1) + k) + u_n(map_w2(df,1) + km))*idz2_w2(k,df) + u_inc(cell_map_w2(df)+k) = viscosity_mu*(d2dx + d2dy + d2dz) end do ! Vertical Velocity diffusion diff --git a/science/gungho/unit-test/kernel/diffusion/leonard_term_u_kernel_mod_test.pf b/science/gungho/unit-test/kernel/diffusion/leonard_term_u_kernel_mod_test.pf index 0feecaaec..ed473de0a 100644 --- a/science/gungho/unit-test/kernel/diffusion/leonard_term_u_kernel_mod_test.pf +++ b/science/gungho/unit-test/kernel/diffusion/leonard_term_u_kernel_mod_test.pf @@ -78,6 +78,8 @@ contains integer(i_def), allocatable :: true_stencil_map_w2(:,:,:) integer(i_def), allocatable :: stencil_map_wt(:,:,:) integer(i_def), allocatable :: stencil_map_pid(:,:,:) + integer(i_def), allocatable :: face_selector_ew(:) + integer(i_def), allocatable :: face_selector_ns(:) real(r_def), allocatable :: u_inc(:) real(r_def), allocatable :: u(:) @@ -160,6 +162,8 @@ contains allocate(dtrdz(undf_w2)) allocate(rho(undf_w2)) allocate(panel_id(undf_pid)) + allocate(face_selector_ew(undf_pid)) + allocate(face_selector_ns(undf_pid)) u_inc(:) = 0.0_r_def u(:) = 1.0_r_def @@ -168,6 +172,8 @@ contains dtrdz(:) = 0.001_r_def rho(:) = 1.0_r_def panel_id = 2.0_r_def + face_selector_ew(:) = 2 + face_selector_ns(:) = 2 cell = 5 k = 1 @@ -203,34 +209,41 @@ contains rho, & panel_id, & 9, stencil_map_pid(:,:,cell), & + face_selector_ew, & + face_selector_ns, & planet_radius, & leonard_kl, & dt, nlayers, & ndf_w2, undf_w2, map_w2(:,cell), & ndf_wt, undf_wt, map_wt(:,cell), & ndf_w1, undf_w1, map_w1(:,cell), & + ndf_pid, undf_pid, map_pid(:,cell), & ndf_pid, undf_pid, map_pid(:,cell) & ) - answer = 0.1152_r_def - @assertEqual(answer, u_inc(map_w2(1,cell) + k), tol) - - deallocate(map_w2) - deallocate(stencil_map_w2) - deallocate(true_map_w2) - deallocate(true_stencil_map_w2) - deallocate(map_wt) - deallocate(stencil_map_wt) - deallocate(map_w1) - deallocate(height_w1) - deallocate(height_w2) - deallocate(u) - deallocate(u_inc) - deallocate(w) - deallocate(w_inc) - deallocate(dtrdz) - deallocate(rho) - deallocate(panel_id) + answer = 0.1152_r_def + @assertEqual(answer, u_inc(map_w2(1,cell) + k), tol) + + deallocate(map_w2) + deallocate(stencil_map_w2) + deallocate(true_map_w2) + deallocate(true_stencil_map_w2) + deallocate(map_wt) + deallocate(stencil_map_wt) + deallocate(map_w1) + deallocate(height_w1) + deallocate(height_w2) + deallocate(u) + deallocate(u_inc) + deallocate(w) + deallocate(w_inc) + deallocate(dtrdz) + deallocate(rho) + deallocate(panel_id) + deallocate(stencil_map_pid) + deallocate(map_pid) + deallocate(face_selector_ew) + deallocate(face_selector_ns) end subroutine test_all diff --git a/science/gungho/unit-test/kernel/diffusion/momentum_smagorinsky_kernel_mod_test.pf b/science/gungho/unit-test/kernel/diffusion/momentum_smagorinsky_kernel_mod_test.pf index 1500f8366..f055c6ba8 100644 --- a/science/gungho/unit-test/kernel/diffusion/momentum_smagorinsky_kernel_mod_test.pf +++ b/science/gungho/unit-test/kernel/diffusion/momentum_smagorinsky_kernel_mod_test.pf @@ -81,6 +81,8 @@ contains integer(i_def), allocatable :: stencil_map_w1(:,:,:) integer(i_def), allocatable :: stencil_map_wt(:,:,:) integer(i_def), allocatable :: stencil_map_pid(:,:,:) + integer(i_def), allocatable :: face_selector_ew(:) + integer(i_def), allocatable :: face_selector_ns(:) real(r_def), allocatable :: u(:), dx_at_w2(:) real(r_def), allocatable :: u_inc(:) @@ -157,11 +159,15 @@ contains allocate(u_inc(undf_w2)) allocate(visc_m(undf_wt)) allocate(panel_id(undf_pid)) + allocate(face_selector_ew(undf_pid)) + allocate(face_selector_ns(undf_pid)) u(:) = 1.0_r_def u_inc(:) = 0.0_r_def visc_m(:) = 1.0_r_def dx_at_w2(:) = dx panel_id = 1.0_r_def + face_selector_ew(:) = 2 + face_selector_ns(:) = 2 ! Create non-uniform wind to generate some diffusion u(stencil_map_w2(2,1,cell)+1) = 2.0_r_def @@ -183,22 +189,24 @@ contains u(true_stencil_map_w2(3,3,cell)+2) = -1.0_r_def ! Call the kernel - call momentum_smagorinsky_code( nlayers, & - u_inc, & - u, & - 5, stencil_map_w2(:,:,cell), & - dx_at_w2, & - 5, stencil_map_w2(:,:,cell), & - height_w2, & - height_w1, & - visc_m, & - 5, stencil_map_wt(:,:,cell), & - panel_id, & - 5, stencil_map_pid(:,:,cell), & - ndf_w2, undf_w2, map_w2(:,cell), & - ndf_w1, undf_w1, map_w1(:,cell), & - ndf_wt, undf_wt, map_wt(:,cell), & - ndf_pid, undf_pid, map_pid(:,cell)& + call momentum_smagorinsky_code( nlayers, & + u_inc, & + u, & + 5, stencil_map_w2(:,:,cell), & + dx_at_w2, & + 5, stencil_map_w2(:,:,cell), & + height_w2, & + height_w1, & + visc_m, & + 5, stencil_map_wt(:,:,cell), & + panel_id, & + 5, stencil_map_pid(:,:,cell), & + face_selector_ew, face_selector_ns, & + ndf_w2, undf_w2, map_w2(:,cell), & + ndf_w1, undf_w1, map_w1(:,cell), & + ndf_wt, undf_wt, map_wt(:,cell), & + ndf_pid, undf_pid, map_pid(:,cell), & + ndf_pid, undf_pid, map_pid(:,cell) & ) answer = -0.4e-5_r_def @@ -220,6 +228,10 @@ contains deallocate(height_w1) deallocate(height_w2) deallocate(panel_id) + deallocate(stencil_map_pid) + deallocate(map_pid) + deallocate(face_selector_ew) + deallocate(face_selector_ns) end subroutine test_all diff --git a/science/gungho/unit-test/kernel/diffusion/momentum_viscosity_kernel_mod_test.pf b/science/gungho/unit-test/kernel/diffusion/momentum_viscosity_kernel_mod_test.pf index 31f0bb9cf..4e2dd034f 100644 --- a/science/gungho/unit-test/kernel/diffusion/momentum_viscosity_kernel_mod_test.pf +++ b/science/gungho/unit-test/kernel/diffusion/momentum_viscosity_kernel_mod_test.pf @@ -96,6 +96,10 @@ contains integer(i_def) :: ndf_w2, undf_w2 integer(i_def) :: dim_space, dim_space_diff integer(i_def) :: nqp_h, nqp_v + integer(i_def) :: ndf_w3_2d, undf_w3_2d + integer(i_def) :: map_w3_2d(1) + integer(i_def) :: face_selector_ew(1) + integer(i_def) :: face_selector_ns(1) integer(i_def), allocatable :: map_w2(:,:) integer(i_def), allocatable :: stencil_map_w2(:,:,:) @@ -112,6 +116,12 @@ contains call get_w2_m3x3_dofmap(map_w2) call get_m3x3_stencil_dofmap_cross(stencil_map_w2, map_w2) + ndf_w3_2d = 1 + undf_w3_2d = 1 + map_w3_2d(1) = 1 + face_selector_ew(1) = 2 + face_selector_ns(1) = 2 + ! Compute coordinates allocate(dx_at_w2(undf_w2)) icell = 1 @@ -134,16 +144,20 @@ contains u_inc(:) = 0.0_r_def cell = 5 - call momentum_viscosity_code( & - nlayers, & - u_inc, & - u, & + call momentum_viscosity_code( & + nlayers, & + u_inc, & + u, & 5, stencil_map_w2(:,:,cell), & - dx_at_w2, & + dx_at_w2, & 5, stencil_map_w2(:,:,cell), & - viscosity_mu, & - ndf_w2, undf_w2, & - map_w2(:,cell) & + face_selector_ew, & + face_selector_ns, & + viscosity_mu, & + ndf_w2, undf_w2, & + map_w2(:,cell), & + ndf_w3_2d, undf_w3_2d, & + map_w3_2d & ) answer = 0.0_r_def diff --git a/science/linear/source/algorithm/timestepping/tl_rk_alg_timestep_mod.x90 b/science/linear/source/algorithm/timestepping/tl_rk_alg_timestep_mod.x90 index 30d391f99..f91137e2b 100644 --- a/science/linear/source/algorithm/timestepping/tl_rk_alg_timestep_mod.x90 +++ b/science/linear/source/algorithm/timestepping/tl_rk_alg_timestep_mod.x90 @@ -25,9 +25,11 @@ module tl_rk_alg_timestep_mod use derived_config_mod, only: bundle_size use sci_fem_constants_mod, only: get_qr_fe, & get_inverse_mass_matrix_fe - use sci_geometric_constants_mod, only: get_coordinates, & - get_panel_id, & - get_dx_at_w2 + use sci_geometric_constants_mod, only: get_coordinates, & + get_panel_id, & + get_dx_at_w2, & + get_face_selector_ew, & + get_face_selector_ns use transport_config_mod, only: operators, & operators_fv use mixing_config_mod, only: viscosity, viscosity_mu @@ -65,6 +67,7 @@ module tl_rk_alg_timestep_mod ! Derived Types use field_mod, only: field_type use field_collection_mod, only: field_collection_type + use integer_field_mod, only: integer_field_type use quadrature_xyoz_mod, only: quadrature_xyoz_type use quadrature_face_mod, only: quadrature_face_type use quadrature_rule_gaussian_mod, only: quadrature_rule_gaussian_type @@ -246,6 +249,8 @@ contains type( field_type ), pointer :: dx_at_w2 => null() type(operator_type), pointer :: m3_inv => null(), & coriolis => null() + type(integer_field_type), pointer :: face_selector_ew + type(integer_field_type), pointer :: face_selector_ns type(quadrature_rule_gaussian_type) :: gaussian_quadrature @@ -493,7 +498,10 @@ contains tracer_viscosity_kernel_type( inc(igh_t), state(igh_t), 1, & dx_at_w2, viscosity_mu ), & momentum_viscosity_kernel_type( inc(igh_u), state(igh_u), 1, & - dx_at_w2, 1, viscosity_mu ), & + dx_at_w2, 1, & + face_selector_ew, & + face_selector_ns, & + viscosity_mu ), & setval_c(inc(igh_d), 0.0_r_def) ) call bundle_axpy(cast_dt, inc, state, state, bundle_size) end if