diff --git a/config_src/drivers/timing_tests/time_MOM_remapping.F90 b/config_src/drivers/timing_tests/time_MOM_remapping.F90 index 1e5aece7bf..e4bea9d94f 100644 --- a/config_src/drivers/timing_tests/time_MOM_remapping.F90 +++ b/config_src/drivers/timing_tests/time_MOM_remapping.F90 @@ -61,11 +61,12 @@ program time_MOM_remapping do isamp = 1, nsamp ! Time reconstruction + remapping do ischeme = 1, nschemes - call initialize_remapping(CS, remapping_scheme=trim(scheme_labels(ischeme))) + call initialize_remapping(CS, remapping_scheme=trim(scheme_labels(ischeme)), & + h_neglect=h_neglect, h_neglect_edge=h_neglect) call cpu_time(start) do iter = 1, nits ! Make many passes to reduce sampling error do ij = 1, nij ! Calling nij times to make similar to cost in MOM_ALE() - call remapping_core_h(CS, nk, h0(:,ij), u0(:,ij), nk, h1(:,ij), u1(:,ij), h_neglect) + call remapping_core_h(CS, nk, h0(:,ij), u0(:,ij), nk, h1(:,ij), u1(:,ij)) enddo enddo call cpu_time(finish) diff --git a/src/ALE/MOM_ALE.F90 b/src/ALE/MOM_ALE.F90 index 0ae5fb1e92..6e9f1c7df9 100644 --- a/src/ALE/MOM_ALE.F90 +++ b/src/ALE/MOM_ALE.F90 @@ -181,6 +181,7 @@ subroutine ALE_init( param_file, GV, US, max_depth, CS) logical :: om4_remap_via_sub_cells type(hybgen_regrid_CS), pointer :: hybgen_regridCS => NULL() ! Control structure for hybgen regridding ! for sharing parameters. + real :: h_neglect, h_neglect_edge ! small thicknesses [H ~> m or kg m-2] if (associated(CS)) then call MOM_error(WARNING, "ALE_init called with an associated "// & @@ -248,20 +249,30 @@ subroutine ALE_init( param_file, GV, US, max_depth, CS) default=default_answer_date, do_not_log=.not.GV%Boussinesq) if (.not.GV%Boussinesq) CS%answer_date = max(CS%answer_date, 20230701) + if (CS%answer_date >= 20190101) then + h_neglect = GV%H_subroundoff ; h_neglect_edge = GV%H_subroundoff + elseif (GV%Boussinesq) then + h_neglect = GV%m_to_H * 1.0e-30 ; h_neglect_edge = GV%m_to_H * 1.0e-10 + else + h_neglect = GV%kg_m2_to_H * 1.0e-30 ; h_neglect_edge = GV%kg_m2_to_H * 1.0e-10 + endif + call initialize_remapping( CS%remapCS, string, & boundary_extrapolation=init_boundary_extrap, & check_reconstruction=check_reconstruction, & check_remapping=check_remapping, & force_bounds_in_subcell=force_bounds_in_subcell, & om4_remap_via_sub_cells=om4_remap_via_sub_cells, & - answer_date=CS%answer_date) + answer_date=CS%answer_date, & + h_neglect=h_neglect, h_neglect_edge=h_neglect_edge) call initialize_remapping( CS%vel_remapCS, vel_string, & boundary_extrapolation=init_boundary_extrap, & check_reconstruction=check_reconstruction, & check_remapping=check_remapping, & force_bounds_in_subcell=force_bounds_in_subcell, & om4_remap_via_sub_cells=om4_remap_via_sub_cells, & - answer_date=CS%answer_date) + answer_date=CS%answer_date, & + h_neglect=h_neglect, h_neglect_edge=h_neglect_edge) call get_param(param_file, mdl, "PARTIAL_CELL_VELOCITY_REMAP", CS%partial_cell_vel_remap, & "If true, use partial cell thicknesses at velocity points that are masked out "//& @@ -653,7 +664,6 @@ subroutine ALE_regrid_accelerated(CS, G, GV, US, h, tv, n_itt, u, v, OBC, Reg, d real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: dzInterface ! Interface height changes within ! an iteration [H ~> m or kg m-2] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: dzIntTotal ! Cumulative interface position changes [H ~> m or kg m-2] - real :: h_neglect, h_neglect_edge ! small thicknesses [H ~> m or kg m-2] nz = GV%ke @@ -680,14 +690,6 @@ subroutine ALE_regrid_accelerated(CS, G, GV, US, h, tv, n_itt, u, v, OBC, Reg, d if (present(dt)) & call ALE_update_regrid_weights(dt, CS) - if (CS%answer_date >= 20190101) then - h_neglect = GV%H_subroundoff ; h_neglect_edge = GV%H_subroundoff - elseif (GV%Boussinesq) then - h_neglect = GV%m_to_H * 1.0e-30 ; h_neglect_edge = GV%m_to_H * 1.0e-10 - else - h_neglect = GV%kg_m2_to_H * 1.0e-30 ; h_neglect_edge = GV%kg_m2_to_H * 1.0e-10 - endif - do itt = 1, n_itt call do_group_pass(pass_T_S_h, G%domain) @@ -704,10 +706,8 @@ subroutine ALE_regrid_accelerated(CS, G, GV, US, h, tv, n_itt, u, v, OBC, Reg, d ! remap from original grid onto new grid do j = G%jsc-1,G%jec+1 ; do i = G%isc-1,G%iec+1 - call remapping_core_h(CS%remapCS, nz, h_orig(i,j,:), tv%S(i,j,:), nz, h(i,j,:), tv_local%S(i,j,:), & - h_neglect, h_neglect_edge) - call remapping_core_h(CS%remapCS, nz, h_orig(i,j,:), tv%T(i,j,:), nz, h(i,j,:), tv_local%T(i,j,:), & - h_neglect, h_neglect_edge) + call remapping_core_h(CS%remapCS, nz, h_orig(i,j,:), tv%S(i,j,:), nz, h(i,j,:), tv_local%S(i,j,:)) + call remapping_core_h(CS%remapCS, nz, h_orig(i,j,:), tv%T(i,j,:), nz, h(i,j,:), tv_local%T(i,j,:)) enddo ; enddo ! starting grid for next iteration @@ -763,7 +763,6 @@ subroutine ALE_remap_tracers(CS, G, GV, h_old, h_new, Reg, debug, dt, PCM_cell) real :: Idt ! The inverse of the timestep [T-1 ~> s-1] real :: h1(GV%ke) ! A column of source grid layer thicknesses [H ~> m or kg m-2] real :: h2(GV%ke) ! A column of target grid layer thicknesses [H ~> m or kg m-2] - real :: h_neglect, h_neglect_edge ! Tiny thicknesses used in remapping [H ~> m or kg m-2] logical :: show_call_tree type(tracer_type), pointer :: Tr => NULL() integer :: i, j, k, m, nz, ntr @@ -771,14 +770,6 @@ subroutine ALE_remap_tracers(CS, G, GV, h_old, h_new, Reg, debug, dt, PCM_cell) show_call_tree = .false. if (present(debug)) show_call_tree = debug - if (CS%answer_date >= 20190101) then - h_neglect = GV%H_subroundoff ; h_neglect_edge = GV%H_subroundoff - elseif (GV%Boussinesq) then - h_neglect = GV%m_to_H*1.0e-30 ; h_neglect_edge = GV%m_to_H*1.0e-10 - else - h_neglect = GV%kg_m2_to_H*1.0e-30 ; h_neglect_edge = GV%kg_m2_to_H*1.0e-10 - endif - if (show_call_tree) call callTree_enter("ALE_remap_tracers(), MOM_ALE.F90") nz = GV%ke @@ -803,11 +794,9 @@ subroutine ALE_remap_tracers(CS, G, GV, h_old, h_new, Reg, debug, dt, PCM_cell) h2(:) = h_new(i,j,:) if (present(PCM_cell)) then PCM(:) = PCM_cell(i,j,:) - call remapping_core_h(CS%remapCS, nz, h1, Tr%t(i,j,:), nz, h2, tr_column, & - h_neglect, h_neglect_edge, PCM_cell=PCM) + call remapping_core_h(CS%remapCS, nz, h1, Tr%t(i,j,:), nz, h2, tr_column, PCM_cell=PCM) else - call remapping_core_h(CS%remapCS, nz, h1, Tr%t(i,j,:), nz, h2, tr_column, & - h_neglect, h_neglect_edge) + call remapping_core_h(CS%remapCS, nz, h1, Tr%t(i,j,:), nz, h2, tr_column) endif ! Possibly underflow any very tiny tracer concentrations to 0. Note that this is not conservative! @@ -1091,7 +1080,6 @@ subroutine ALE_remap_velocities(CS, G, GV, h_old_u, h_old_v, h_new_u, h_new_v, u real :: v_tgt(GV%ke) ! A column of v-velocities on the target grid [L T-1 ~> m s-1] real :: h1(GV%ke) ! A column of source grid layer thicknesses [H ~> m or kg m-2] real :: h2(GV%ke) ! A column of target grid layer thicknesses [H ~> m or kg m-2] - real :: h_neglect, h_neglect_edge ! Tiny thicknesses used in remapping [H ~> m or kg m-2] logical :: show_call_tree integer :: i, j, k, nz @@ -1099,14 +1087,6 @@ subroutine ALE_remap_velocities(CS, G, GV, h_old_u, h_old_v, h_new_u, h_new_v, u if (present(debug)) show_call_tree = debug if (show_call_tree) call callTree_enter("ALE_remap_velocities()") - if (CS%answer_date >= 20190101) then - h_neglect = GV%H_subroundoff ; h_neglect_edge = GV%H_subroundoff - elseif (GV%Boussinesq) then - h_neglect = GV%m_to_H*1.0e-30 ; h_neglect_edge = GV%m_to_H*1.0e-10 - else - h_neglect = GV%kg_m2_to_H*1.0e-30 ; h_neglect_edge = GV%kg_m2_to_H*1.0e-10 - endif - nz = GV%ke ! --- Remap u profiles from the source vertical grid onto the new target grid. @@ -1120,8 +1100,7 @@ subroutine ALE_remap_velocities(CS, G, GV, h_old_u, h_old_v, h_new_u, h_new_v, u u_src(k) = u(I,j,k) enddo - call remapping_core_h(CS%vel_remapCS, nz, h1, u_src, nz, h2, u_tgt, & - h_neglect, h_neglect_edge) + call remapping_core_h(CS%vel_remapCS, nz, h1, u_src, nz, h2, u_tgt) if ((CS%BBL_h_vel_mask > 0.0) .and. (CS%h_vel_mask > 0.0)) & call mask_near_bottom_vel(u_tgt, h2, CS%BBL_h_vel_mask, CS%h_vel_mask, nz) @@ -1146,8 +1125,7 @@ subroutine ALE_remap_velocities(CS, G, GV, h_old_u, h_old_v, h_new_u, h_new_v, u v_src(k) = v(i,J,k) enddo - call remapping_core_h(CS%vel_remapCS, nz, h1, v_src, nz, h2, v_tgt, & - h_neglect, h_neglect_edge) + call remapping_core_h(CS%vel_remapCS, nz, h1, v_src, nz, h2, v_tgt) if ((CS%BBL_h_vel_mask > 0.0) .and. (CS%h_vel_mask > 0.0)) then call mask_near_bottom_vel(v_tgt, h2, CS%BBL_h_vel_mask, CS%h_vel_mask, nz) @@ -1301,7 +1279,7 @@ end subroutine mask_near_bottom_vel !! h_dst must be dimensioned as a model array with GV%ke layers while h_src can !! have an arbitrary number of layers specified by nk_src. subroutine ALE_remap_scalar(CS, G, GV, nk_src, h_src, s_src, h_dst, s_dst, all_cells, old_remap, & - answers_2018, answer_date, h_neglect, h_neglect_edge) + answers_2018, answer_date) type(remapping_CS), intent(in) :: CS !< Remapping control structure type(ocean_grid_type), intent(in) :: G !< Ocean grid structure type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure @@ -1325,16 +1303,9 @@ subroutine ALE_remap_scalar(CS, G, GV, nk_src, h_src, s_src, h_dst, s_dst, all_c !! use more robust forms of the same expressions. integer, optional, intent(in) :: answer_date !< The vintage of the expressions to use !! for remapping - real, optional, intent(in) :: h_neglect !< A negligibly small thickness used in - !! remapping cell reconstructions, in the same - !! units as h_src, often [H ~> m or kg m-2] - real, optional, intent(in) :: h_neglect_edge !< A negligibly small thickness used in - !! remapping edge value calculations, in the same - !! units as h_src, often [H ~> m or kg m-2] ! Local variables integer :: i, j, k, n_points real :: dx(GV%ke+1) ! Change in interface position [H ~> m or kg m-2] - real :: h_neg, h_neg_edge ! Tiny thicknesses used in remapping [H ~> m or kg m-2] logical :: ignore_vanished_layers, use_remapping_core_w, use_2018_remap ignore_vanished_layers = .false. @@ -1345,19 +1316,6 @@ subroutine ALE_remap_scalar(CS, G, GV, nk_src, h_src, s_src, h_dst, s_dst, all_c use_2018_remap = .true. ; if (present(answers_2018)) use_2018_remap = answers_2018 if (present(answer_date)) use_2018_remap = (answer_date < 20190101) - if (present(h_neglect)) then - h_neg = h_neglect - h_neg_edge = h_neg ; if (present(h_neglect_edge)) h_neg_edge = h_neglect_edge - else - if (.not.use_2018_remap) then - h_neg = GV%H_subroundoff ; h_neg_edge = GV%H_subroundoff - elseif (GV%Boussinesq) then - h_neg = GV%m_to_H*1.0e-30 ; h_neg_edge = GV%m_to_H*1.0e-10 - else - h_neg = GV%kg_m2_to_H*1.0e-30 ; h_neg_edge = GV%kg_m2_to_H*1.0e-10 - endif - endif - !$OMP parallel do default(shared) firstprivate(n_points,dx) do j = G%jsc,G%jec ; do i = G%isc,G%iec if (G%mask2dT(i,j) > 0.) then @@ -1371,10 +1329,10 @@ subroutine ALE_remap_scalar(CS, G, GV, nk_src, h_src, s_src, h_dst, s_dst, all_c if (use_remapping_core_w) then call dzFromH1H2( n_points, h_src(i,j,1:n_points), GV%ke, h_dst(i,j,:), dx ) call remapping_core_w(CS, n_points, h_src(i,j,1:n_points), s_src(i,j,1:n_points), & - GV%ke, dx, s_dst(i,j,:), h_neg, h_neg_edge) + GV%ke, dx, s_dst(i,j,:)) else call remapping_core_h(CS, n_points, h_src(i,j,1:n_points), s_src(i,j,1:n_points), & - GV%ke, h_dst(i,j,:), s_dst(i,j,:), h_neg, h_neg_edge) + GV%ke, h_dst(i,j,:), s_dst(i,j,:)) endif else s_dst(i,j,:) = 0. diff --git a/src/ALE/MOM_remapping.F90 b/src/ALE/MOM_remapping.F90 index f2f0841fa0..68a5973ba9 100644 --- a/src/ALE/MOM_remapping.F90 +++ b/src/ALE/MOM_remapping.F90 @@ -40,6 +40,11 @@ module MOM_remapping !> If true, use the OM4 version of the remapping algorithm that makes poor assumptions !! about the reconstructions in top and bottom layers of the source grid logical :: om4_remap_via_sub_cells = .false. + + !> A negligibly small width for the purpose of cell reconstructions in the same units as h0 [H] + real :: h_neglect + !> A negligibly small width for the purpose of edge value calculations in the same units as h0 [H]. + real :: h_neglect_edge end type !> Class to assist in unit tests @@ -114,7 +119,8 @@ module MOM_remapping !> Set parameters within remapping object subroutine remapping_set_param(CS, remapping_scheme, boundary_extrapolation, & check_reconstruction, check_remapping, force_bounds_in_subcell, & - om4_remap_via_sub_cells, answers_2018, answer_date) + om4_remap_via_sub_cells, answers_2018, answer_date, & + h_neglect, h_neglect_edge) type(remapping_CS), intent(inout) :: CS !< Remapping control structure character(len=*), optional, intent(in) :: remapping_scheme !< Remapping scheme to use logical, optional, intent(in) :: boundary_extrapolation !< Indicate to extrapolate in boundary cells @@ -124,6 +130,10 @@ subroutine remapping_set_param(CS, remapping_scheme, boundary_extrapolation, & logical, optional, intent(in) :: om4_remap_via_sub_cells !< If true, use OM4 remapping algorithm logical, optional, intent(in) :: answers_2018 !< If true use older, less accurate expressions. integer, optional, intent(in) :: answer_date !< The vintage of the expressions to use + real, optional, intent(in) :: h_neglect !< A negligibly small width for the purpose of cell + !! reconstructions in the same units as h0 [H] + real, optional, intent(in) :: h_neglect_edge !< A negligibly small width for the purpose of edge + !! value calculations in the same units as h0 [H]. if (present(remapping_scheme)) then call setReconstructionType( remapping_scheme, CS ) @@ -153,6 +163,12 @@ subroutine remapping_set_param(CS, remapping_scheme, boundary_extrapolation, & if (present(answer_date)) then CS%answer_date = answer_date endif + if (present(h_neglect)) then + CS%h_neglect = h_neglect + endif + if (present(h_neglect_edge)) then + CS%h_neglect_edge = h_neglect_edge + endif end subroutine remapping_set_param @@ -182,7 +198,7 @@ end subroutine extract_member_remapping_CS !! \todo Remove h_neglect argument by moving into remapping_CS !! \todo Remove PCM_cell argument by adding new method in Recon1D class !! \todo Inline the two versions of remap_via_sub_cells() in remapping_core_h() to eliminate remap_via_sub_cells() -subroutine remapping_core_h(CS, n0, h0, u0, n1, h1, u1, h_neglect, h_neglect_edge, PCM_cell) +subroutine remapping_core_h(CS, n0, h0, u0, n1, h1, u1, PCM_cell) type(remapping_CS), intent(in) :: CS !< Remapping control structure integer, intent(in) :: n0 !< Number of cells on source grid real, dimension(n0), intent(in) :: h0 !< Cell widths on source grid [H] @@ -190,12 +206,6 @@ subroutine remapping_core_h(CS, n0, h0, u0, n1, h1, u1, h_neglect, h_neglect_edg integer, intent(in) :: n1 !< Number of cells on target grid real, dimension(n1), intent(in) :: h1 !< Cell widths on target grid [H] real, dimension(n1), intent(out) :: u1 !< Cell averages on target grid [A] - real, intent(in) :: h_neglect !< A negligibly small width for the - !! purpose of cell reconstructions - !! in the same units as h0 [H] - real, optional, intent(in) :: h_neglect_edge !< A negligibly small width for the purpose - !! of edge value calculations in the same units as h0 [H]. - !! The default is h_neglect. logical, dimension(n0), optional, intent(in) :: PCM_cell !< If present, use PCM remapping for !! cells in the source grid where this is true. @@ -204,20 +214,18 @@ subroutine remapping_core_h(CS, n0, h0, u0, n1, h1, u1, h_neglect, h_neglect_edg real, dimension(n0,2) :: ppoly_r_S ! Edge slope of polynomial [A H-1] real, dimension(n0,CS%degree+1) :: ppoly_r_coefs ! Coefficients of polynomial reconstructions [A] real :: uh_err ! A bound on the error in the sum of u*h, as estimated by the remapping code [A H] - real :: hNeglect_edge ! Negligibly small cell widths in the same units as h0 [H] integer :: iMethod ! An integer indicating the integration method used - hNeglect_edge = h_neglect ; if (present(h_neglect_edge)) hNeglect_edge = h_neglect_edge call build_reconstructions_1d( CS, n0, h0, u0, ppoly_r_coefs, ppoly_r_E, ppoly_r_S, iMethod, & - h_neglect, hNeglect_edge, PCM_cell ) + CS%h_neglect, CS%h_neglect_edge, PCM_cell ) if (CS%om4_remap_via_sub_cells) then if (CS%check_reconstruction) call check_reconstructions_1d(n0, h0, u0, CS%degree, & CS%boundary_extrapolation, ppoly_r_coefs, ppoly_r_E) - ! This calls the OM4 version of the remapmping algorithms + ! This calls the OM4 version of the remapping algorithms call remap_via_sub_cells_om4( n0, h0, u0, ppoly_r_E, ppoly_r_coefs, n1, h1, iMethod, & CS%force_bounds_in_subcell, u1, uh_err ) @@ -235,7 +243,7 @@ end subroutine remapping_core_h !> Remaps column of values u0 on grid h0 to implied grid h1 !! where the interfaces of h1 differ from those of h0 by dx. -subroutine remapping_core_w( CS, n0, h0, u0, n1, dx, u1, h_neglect, h_neglect_edge ) +subroutine remapping_core_w( CS, n0, h0, u0, n1, dx, u1) type(remapping_CS), intent(in) :: CS !< Remapping control structure integer, intent(in) :: n0 !< Number of cells on source grid real, dimension(n0), intent(in) :: h0 !< Cell widths on source grid [H] @@ -243,26 +251,17 @@ subroutine remapping_core_w( CS, n0, h0, u0, n1, dx, u1, h_neglect, h_neglect_ed integer, intent(in) :: n1 !< Number of cells on target grid real, dimension(n1+1), intent(in) :: dx !< Cell widths on target grid [H] real, dimension(n1), intent(out) :: u1 !< Cell averages on target grid [A] - real, intent(in) :: h_neglect !< A negligibly small width for the - !! purpose of cell reconstructions - !! in the same units as h0 [H]. - real, optional, intent(in) :: h_neglect_edge !< A negligibly small width for the purpose - !! of edge value calculations in the same units as h0 [H]. - !! The default is h_neglect. ! Local variables real, dimension(n0,2) :: ppoly_r_E ! Edge value of polynomial [A] real, dimension(n0,2) :: ppoly_r_S ! Edge slope of polynomial [A H-1] real, dimension(n0,CS%degree+1) :: ppoly_r_coefs ! Coefficients of polynomial reconstructions [A] real, dimension(n1) :: h1 !< Cell widths on target grid [H] real :: uh_err ! A bound on the error in the sum of u*h, as estimated by the remapping code [A H] - real :: hNeglect_edge ! Negligibly small thicknesses [H] integer :: iMethod ! An integer indicating the integration method used integer :: k - hNeglect_edge = h_neglect ; if (present(h_neglect_edge)) hNeglect_edge = h_neglect_edge - call build_reconstructions_1d( CS, n0, h0, u0, ppoly_r_coefs, ppoly_r_E, ppoly_r_S, iMethod,& - h_neglect, hNeglect_edge ) + CS%h_neglect, CS%h_neglect_edge ) if (CS%check_reconstruction) call check_reconstructions_1d(n0, h0, u0, CS%degree, & CS%boundary_extrapolation, ppoly_r_coefs, ppoly_r_E) @@ -1520,7 +1519,8 @@ end subroutine dzFromH1H2 !> Constructor for remapping control structure subroutine initialize_remapping( CS, remapping_scheme, boundary_extrapolation, & check_reconstruction, check_remapping, force_bounds_in_subcell, & - om4_remap_via_sub_cells, answers_2018, answer_date) + om4_remap_via_sub_cells, answers_2018, answer_date, & + h_neglect, h_neglect_edge) ! Arguments type(remapping_CS), intent(inout) :: CS !< Remapping control structure character(len=*), intent(in) :: remapping_scheme !< Remapping scheme to use @@ -1531,12 +1531,17 @@ subroutine initialize_remapping( CS, remapping_scheme, boundary_extrapolation, & logical, optional, intent(in) :: om4_remap_via_sub_cells !< If true, use OM4 remapping algorithm logical, optional, intent(in) :: answers_2018 !< If true use older, less accurate expressions. integer, optional, intent(in) :: answer_date !< The vintage of the expressions to use + real, optional, intent(in) :: h_neglect !< A negligibly small width for the purpose of cell + !! reconstructions in the same units as h0 [H] + real, optional, intent(in) :: h_neglect_edge !< A negligibly small width for the purpose of edge + !! value calculations in the same units as h0 [H]. ! Note that remapping_scheme is mandatory for initialize_remapping() call remapping_set_param(CS, remapping_scheme=remapping_scheme, boundary_extrapolation=boundary_extrapolation, & check_reconstruction=check_reconstruction, check_remapping=check_remapping, & force_bounds_in_subcell=force_bounds_in_subcell, & - om4_remap_via_sub_cells=om4_remap_via_sub_cells, answers_2018=answers_2018, answer_date=answer_date) + om4_remap_via_sub_cells=om4_remap_via_sub_cells, answers_2018=answers_2018, answer_date=answer_date, & + h_neglect=h_neglect, h_neglect_edge=h_neglect_edge) end subroutine initialize_remapping @@ -1636,7 +1641,8 @@ logical function remapping_unit_tests(verbose) if (verbose) write(test%stdout,*) ' - - - - - 1st generation tests - - - - -' - call initialize_remapping(CS, 'PPM_H4', answer_date=answer_date) + call initialize_remapping(CS, 'PPM_H4', answer_date=answer_date, & + h_neglect=h_neglect, h_neglect_edge=h_neglect_edge) ! Profile 0: 4 layers of thickness 0.75 and total depth 3, with du/dz=8 n0 = 4 @@ -1656,7 +1662,7 @@ logical function remapping_unit_tests(verbose) ! Mapping u1 from h1 to h2 call dzFromH1H2( n0, h0, n1, h1, dx1 ) - call remapping_core_w( CS, n0, h0, u0, n1, dx1, u1, h_neglect, h_neglect_edge) + call remapping_core_w( CS, n0, h0, u0, n1, dx1, u1) call test%real_arr(3, u1, (/8.,0.,-8./), 'remapping_core_w() PPM_H4') allocate(ppoly0_E(n0,2), ppoly0_S(n0,2), ppoly0_coefs(n0,CS%degree+1)) @@ -2120,7 +2126,7 @@ logical function remapping_unit_tests(verbose) u0 = (/1.0, 1.5, 2.5, 3.5, 4.5, 5.5, 6.0, 6.0/) allocate( u1(8) ) - call initialize_remapping(CS, 'PLM', answer_date=99990101) + call initialize_remapping(CS, 'PLM', answer_date=99990101, h_neglect=1.e-17, h_neglect_edge=1.e-2) do om4 = 0, 1 if ( om4 == 0 ) then @@ -2132,27 +2138,27 @@ logical function remapping_unit_tests(verbose) endif ! Unchanged grid - call remapping_core_h( CS, n0, h0, u0, 8, [0.,1.,1.,1.,1.,1.,0.,0.], u1, 1.e-17, 1.e-2) + call remapping_core_h( CS, n0, h0, u0, 8, [0.,1.,1.,1.,1.,1.,0.,0.], u1) call test%real_arr(8, u1, (/1.0,1.5,2.5,3.5,4.5,5.5,6.0,6.0/), 'PLM: remapped h=01111100->h=01111100'//om4_tag) ! Removing vanished layers (unchanged values for non-vanished layers, layer centers 0.5, 1.5, 2.5, 3.5, 4.5) - call remapping_core_h( CS, n0, h0, u0, 5, [1.,1.,1.,1.,1.], u1, 1.e-17, 1.e-2) + call remapping_core_h( CS, n0, h0, u0, 5, [1.,1.,1.,1.,1.], u1) call test%real_arr(5, u1, (/1.5,2.5,3.5,4.5,5.5/), 'PLM: remapped h=01111100->h=11111'//om4_tag) ! Remapping to variable thickness layers (layer centers 0.25, 1.0, 2.25, 4.0) - call remapping_core_h( CS, n0, h0, u0, 4, [0.5,1.,1.5,2.], u1, 1.e-17, 1.e-2) + call remapping_core_h( CS, n0, h0, u0, 4, [0.5,1.,1.5,2.], u1) call test%real_arr(4, u1, (/1.25,2.,3.25,5./), 'PLM: remapped h=01111100->h=h1t2'//om4_tag) ! Remapping to variable thickness + vanished layers (layer centers 0.25, 1.0, 1.5, 2.25, 4.0) - call remapping_core_h( CS, n0, h0, u0, 6, [0.5,1.,0.,1.5,2.,0.], u1, 1.e-17, 1.e-2) + call remapping_core_h( CS, n0, h0, u0, 6, [0.5,1.,0.,1.5,2.,0.], u1) call test%real_arr(6, u1, (/1.25,2.,2.5,3.25,5.,6./), 'PLM: remapped h=01111100->h=h10t20'//om4_tag) ! Remapping to deeper water column (layer centers 0.75, 2.25, 3., 5., 8.) - call remapping_core_h( CS, n0, h0, u0, 5, [1.5,1.5,0.,4.,2.], u1, 1.e-17, 1.e-2) + call remapping_core_h( CS, n0, h0, u0, 5, [1.5,1.5,0.,4.,2.], u1) call test%real_arr(5, u1, (/1.75,3.25,4.,5.5,6./), 'PLM: remapped h=01111100->h=tt02'//om4_tag) ! Remapping to slightly shorter water column (layer centers 0.5, 1.5, 2.5,, 3.5, 4.25) - call remapping_core_h( CS, n0, h0, u0, 5, [1.,1.,1.,1.,0.5], u1, 1.e-17, 1.e-2) + call remapping_core_h( CS, n0, h0, u0, 5, [1.,1.,1.,1.,0.5], u1) if ( om4 == 0 ) then call test%real_arr(5, u1, (/1.5,2.5,3.5,4.5,5.25/), 'PLM: remapped h=01111100->h=1111h') else @@ -2160,7 +2166,7 @@ logical function remapping_unit_tests(verbose) endif ! Remapping to much shorter water column (layer centers 0.25, 0.5, 1.) - call remapping_core_h( CS, n0, h0, u0, 3, [0.5,0.,1.], u1, 1.e-17, 1.e-2) + call remapping_core_h( CS, n0, h0, u0, 3, [0.5,0.,1.], u1) if ( om4 == 0 ) then call test%real_arr(3, u1, (/1.25,1.5,2./), 'PLM: remapped h=01111100->h=h01') else diff --git a/src/ALE/coord_hycom.F90 b/src/ALE/coord_hycom.F90 index af0f999064..1e5474770a 100644 --- a/src/ALE/coord_hycom.F90 +++ b/src/ALE/coord_hycom.F90 @@ -175,8 +175,8 @@ subroutine build_hycom1_column(CS, remapCS, eqn_of_state, nz, depth, h, T, S, p_ ( p_col(nz) - p_col(1) ) enddo ! Remap from original h and T,S to get T,S_col_new - call remapping_core_h(remapCS, nz, h(:), T, CS%nk, h_col_new, T_col_new, h_neglect, h_neglect_edge) - call remapping_core_h(remapCS, nz, h(:), S, CS%nk, h_col_new, S_col_new, h_neglect, h_neglect_edge) + call remapping_core_h(remapCS, nz, h(:), T, CS%nk, h_col_new, T_col_new) + call remapping_core_h(remapCS, nz, h(:), S, CS%nk, h_col_new, S_col_new) call build_hycom1_target_anomaly(CS, remapCS, eqn_of_state, CS%nk, depth, & h_col_new, T_col_new, S_col_new, p_col_new, r_col_new, RiA_new, h_neglect, h_neglect_edge) do k= 2,CS%nk diff --git a/src/ALE/coord_rho.F90 b/src/ALE/coord_rho.F90 index 391e4a5dc9..c967687dc8 100644 --- a/src/ALE/coord_rho.F90 +++ b/src/ALE/coord_rho.F90 @@ -272,9 +272,9 @@ subroutine build_rho_column_iteratively(CS, remapCS, nz, depth, h, T, S, eqn_of_ h1(k) = x1(k+1) - x1(k) enddo - call remapping_core_h(remapCS, nz, h0, S, nz, h1, S_tmp, h_neglect, h_neglect_edge) + call remapping_core_h(remapCS, nz, h0, S, nz, h1, S_tmp) - call remapping_core_h(remapCS, nz, h0, T, nz, h1, T_tmp, h_neglect, h_neglect_edge) + call remapping_core_h(remapCS, nz, h0, T, nz, h1, T_tmp) ! Compute the deviation between two successive grids deviation = 0.0 diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 8d45114a39..acd1b945de 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -113,6 +113,7 @@ module MOM use MOM_open_boundary, only : ocean_OBC_type, OBC_registry_type use MOM_open_boundary, only : register_temp_salt_segments, update_segment_tracer_reservoirs use MOM_open_boundary, only : open_boundary_register_restarts, remap_OBC_fields +use MOM_open_boundary, only : open_boundary_setup_vert use MOM_open_boundary, only : rotate_OBC_config, rotate_OBC_init use MOM_porous_barriers, only : porous_widths_layer, porous_widths_interface, porous_barriers_init use MOM_porous_barriers, only : porous_barrier_CS @@ -2652,6 +2653,9 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, & endif CS%HFrz = (US%Z_to_m * GV%m_to_H) * HFrz_z + ! Finish OBC configuration that depend on the vertical grid + call open_boundary_setup_vert(GV, US, OBC_in) + ! Shift from using the temporary dynamic grid type to using the final (potentially static) ! and properly rotated ocean-specific grid type and horizontal index type. if (CS%rotate_index) then diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index 9b8d26cb09..680d80fb1d 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -41,6 +41,7 @@ module MOM_open_boundary public open_boundary_apply_normal_flow public open_boundary_config +public open_boundary_setup_vert public open_boundary_init public open_boundary_query public open_boundary_end @@ -382,6 +383,11 @@ module MOM_open_boundary !! for remapping. Values below 20190101 recover the remapping !! answers from 2018, while higher values use more robust !! forms of the same remapping expressions. + logical :: check_reconstruction !< Flag for remapping to run checks on reconstruction + logical :: check_remapping !< Flag for remapping to run internal checks + logical :: force_bounds_in_subcell !< Flag for remapping to hide overshoot using bounds + logical :: om4_remap_via_sub_cells !< If true, use the OM4 remapping algorithm + character(40) :: remappingScheme !< String selecting the vertical remapping scheme type(group_pass_type) :: pass_oblique !< Structure for group halo pass end type ocean_OBC_type @@ -425,7 +431,6 @@ module MOM_open_boundary !> and ALE_init. Therefore segment data are not fully initialized !> here. The remainder of the segment data are initialized in a !> later call to update_open_boundary_data - subroutine open_boundary_config(G, US, param_file, OBC) type(dyn_horgrid_type), intent(inout) :: G !< Ocean grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type @@ -439,6 +444,7 @@ subroutine open_boundary_config(G, US, param_file, OBC) character(len=1024) :: segment_str ! The contents (rhs) for parameter "segment_param_str" character(len=200) :: config1 ! String for OBC_USER_CONFIG real :: Lscale_in, Lscale_out ! parameters controlling tracer values at the boundaries [L ~> m] + real :: dz_neglect, dz_neglect_edge ! Small thicknesses [Z ~> m] integer :: default_answer_date ! The default setting for the various ANSWER_DATE flags. logical :: check_reconstruction, check_remapping, force_bounds_in_subcell logical :: om4_remap_via_sub_cells ! If true, use the OM4 remapping algorithm @@ -667,23 +673,23 @@ subroutine open_boundary_config(G, US, param_file, OBC) if (Lscale_out>0.) OBC%segment(l)%Tr_InvLscale_out = 1.0/Lscale_out enddo - call get_param(param_file, mdl, "REMAPPING_SCHEME", remappingScheme, & + call get_param(param_file, mdl, "REMAPPING_SCHEME", OBC%remappingScheme, & "This sets the reconstruction scheme used "//& "for vertical remapping for all variables. "//& "It can be one of the following schemes: \n"//& trim(remappingSchemesDoc), default=remappingDefaultScheme,do_not_log=.true.) - call get_param(param_file, mdl, "FATAL_CHECK_RECONSTRUCTIONS", check_reconstruction, & + call get_param(param_file, mdl, "FATAL_CHECK_RECONSTRUCTIONS", OBC%check_reconstruction, & "If true, cell-by-cell reconstructions are checked for "//& "consistency and if non-monotonicity or an inconsistency is "//& "detected then a FATAL error is issued.", default=.false.,do_not_log=.true.) - call get_param(param_file, mdl, "FATAL_CHECK_REMAPPING", check_remapping, & + call get_param(param_file, mdl, "FATAL_CHECK_REMAPPING", OBC%check_remapping, & "If true, the results of remapping are checked for "//& "conservation and new extrema and if an inconsistency is "//& "detected then a FATAL error is issued.", default=.false.,do_not_log=.true.) call get_param(param_file, mdl, "BRUSHCUTTER_MODE", OBC%brushcutter_mode, & "If true, read external OBC data on the supergrid.", & default=.false.) - call get_param(param_file, mdl, "REMAP_BOUND_INTERMEDIATE_VALUES", force_bounds_in_subcell, & + call get_param(param_file, mdl, "REMAP_BOUND_INTERMEDIATE_VALUES", OBC%force_bounds_in_subcell, & "If true, the values on the intermediate grid used for remapping "//& "are forced to be bounded, which might not be the case due to "//& "round off.", default=.false.,do_not_log=.true.) @@ -696,17 +702,11 @@ subroutine open_boundary_config(G, US, param_file, OBC) "that were in use at the end of 2018. Higher values result in the use of more "//& "robust and accurate forms of mathematically equivalent expressions.", & default=default_answer_date) - call get_param(param_file, mdl, "OBC_REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, & + call get_param(param_file, mdl, "OBC_REMAPPING_USE_OM4_SUBCELLS", OBC%om4_remap_via_sub_cells, & "If true, use the OM4 remapping-via-subcells algorithm for neutral diffusion. "//& "See REMAPPING_USE_OM4_SUBCELLS for more details. "//& "We recommend setting this option to false.", default=.true.) - allocate(OBC%remap_CS) - call initialize_remapping(OBC%remap_CS, remappingScheme, boundary_extrapolation = .false., & - check_reconstruction=check_reconstruction, check_remapping=check_remapping, & - om4_remap_via_sub_cells=om4_remap_via_sub_cells, & - force_bounds_in_subcell=force_bounds_in_subcell, answer_date=OBC%remap_answer_date) - endif ! OBC%number_of_segments > 0 ! Safety check @@ -729,6 +729,35 @@ subroutine open_boundary_config(G, US, param_file, OBC) end subroutine open_boundary_config +!> Setup vertical remaping for open boundaries +subroutine open_boundary_setup_vert(GV, US, OBC) + type(verticalGrid_type), intent(in) :: GV !< Container for vertical grid information + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + type(ocean_OBC_type), pointer :: OBC !< Open boundary control structure + + ! Local variables + real :: dz_neglect, dz_neglect_edge ! Small thicknesses [Z ~> m] + + if (associated(OBC)) then + if (OBC%number_of_segments > 0) then + if (GV%Boussinesq .and. (OBC%remap_answer_date < 20190101)) then + dz_neglect = US%m_to_Z * 1.0e-30 ; dz_neglect_edge = US%m_to_Z * 1.0e-10 + elseif (GV%semi_Boussinesq .and. (OBC%remap_answer_date < 20190101)) then + dz_neglect = GV%kg_m2_to_H*GV%H_to_Z * 1.0e-30 ; dz_neglect_edge = GV%kg_m2_to_H*GV%H_to_Z * 1.0e-10 + else + dz_neglect = GV%dZ_subroundoff ; dz_neglect_edge = GV%dZ_subroundoff + endif + allocate(OBC%remap_CS) + call initialize_remapping(OBC%remap_CS, OBC%remappingScheme, boundary_extrapolation = .false., & + check_reconstruction=OBC%check_reconstruction, check_remapping=OBC%check_remapping, & + om4_remap_via_sub_cells=OBC%om4_remap_via_sub_cells, & + force_bounds_in_subcell=OBC%force_bounds_in_subcell, answer_date=OBC%remap_answer_date, & + h_neglect=dz_neglect, h_neglect_edge=dz_neglect_edge) + endif + endif + +end subroutine open_boundary_setup_vert + !> Allocate space for reading OBC data from files. It sets up the required vertical !! remapping. In the process, it does funky stuff with the MPI processes. subroutine initialize_segment_data(G, GV, US, OBC, PF) @@ -3867,7 +3896,6 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) real, allocatable :: normal_trans_bt(:,:) ! barotropic transport [H L2 T-1 ~> m3 s-1] integer :: turns ! Number of index quarter turns real :: time_delta ! Time since tidal reference date [T ~> s] - real :: dz_neglect, dz_neglect_edge ! Small thicknesses [Z ~> m] is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed @@ -3880,14 +3908,6 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) if (OBC%add_tide_constituents) time_delta = US%s_to_T * time_type_to_real(Time - OBC%time_ref) - if (GV%Boussinesq .and. (OBC%remap_answer_date < 20190101)) then - dz_neglect = US%m_to_Z * 1.0e-30 ; dz_neglect_edge = US%m_to_Z * 1.0e-10 - elseif (GV%semi_Boussinesq .and. (OBC%remap_answer_date < 20190101)) then - dz_neglect = GV%kg_m2_to_H*GV%H_to_Z * 1.0e-30 ; dz_neglect_edge = GV%kg_m2_to_H*GV%H_to_Z * 1.0e-10 - else - dz_neglect = GV%dZ_subroundoff ; dz_neglect_edge = GV%dZ_subroundoff - endif - if (OBC%number_of_segments >= 1) then call thickness_to_dz(h, tv, dz, G, GV, US) call pass_var(dz, G%Domain) @@ -4179,22 +4199,19 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) call remapping_core_h(OBC%remap_CS, & segment%field(m)%nk_src, segment%field(m)%dz_src(I,J,:), & segment%field(m)%buffer_src(I,J,:), & - GV%ke, dz_stack, segment%field(m)%buffer_dst(I,J,:), & - dz_neglect, dz_neglect_edge) + GV%ke, dz_stack, segment%field(m)%buffer_dst(I,J,:)) elseif (G%mask2dCu(I,j)>0.) then dz_stack(:) = dz(i+ishift,j,:) call remapping_core_h(OBC%remap_CS, & segment%field(m)%nk_src, segment%field(m)%dz_src(I,J,:), & segment%field(m)%buffer_src(I,J,:), & - GV%ke, dz_stack, segment%field(m)%buffer_dst(I,J,:), & - dz_neglect, dz_neglect_edge) + GV%ke, dz_stack, segment%field(m)%buffer_dst(I,J,:)) elseif (G%mask2dCu(I,j+1)>0.) then dz_stack(:) = dz(i+ishift,j+1,:) call remapping_core_h(OBC%remap_CS, & segment%field(m)%nk_src, segment%field(m)%dz_src(I,j,:), & segment%field(m)%buffer_src(I,J,:), & - GV%ke, dz_stack, segment%field(m)%buffer_dst(I,J,:), & - dz_neglect, dz_neglect_edge) + GV%ke, dz_stack, segment%field(m)%buffer_dst(I,J,:)) endif enddo else @@ -4209,8 +4226,7 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) call remapping_core_h(OBC%remap_CS, & segment%field(m)%nk_src, scl_fac*segment%field(m)%dz_src(I,j,:), & segment%field(m)%buffer_src(I,j,:), & - GV%ke, dz(i+ishift,j,:), segment%field(m)%buffer_dst(I,j,:), & - dz_neglect, dz_neglect_edge) + GV%ke, dz(i+ishift,j,:), segment%field(m)%buffer_dst(I,j,:)) endif enddo endif @@ -4229,22 +4245,19 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) call remapping_core_h(OBC%remap_CS, & segment%field(m)%nk_src, segment%field(m)%dz_src(I,J,:), & segment%field(m)%buffer_src(I,J,:), & - GV%ke, dz_stack, segment%field(m)%buffer_dst(I,J,:), & - dz_neglect, dz_neglect_edge) + GV%ke, dz_stack, segment%field(m)%buffer_dst(I,J,:)) elseif (G%mask2dCv(i,J)>0.) then dz_stack(:) = dz(i,j+jshift,:) call remapping_core_h(OBC%remap_CS, & segment%field(m)%nk_src, segment%field(m)%dz_src(I,J,:), & segment%field(m)%buffer_src(I,J,:), & - GV%ke, dz_stack, segment%field(m)%buffer_dst(I,J,:), & - dz_neglect, dz_neglect_edge) + GV%ke, dz_stack, segment%field(m)%buffer_dst(I,J,:)) elseif (G%mask2dCv(i+1,J)>0.) then dz_stack(:) = dz(i+1,j+jshift,:) call remapping_core_h(OBC%remap_CS, & segment%field(m)%nk_src, segment%field(m)%dz_src(I,J,:), & segment%field(m)%buffer_src(I,J,:), & - GV%ke, dz_stack, segment%field(m)%buffer_dst(I,J,:), & - dz_neglect, dz_neglect_edge) + GV%ke, dz_stack, segment%field(m)%buffer_dst(I,J,:)) endif enddo else @@ -4259,8 +4272,7 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) call remapping_core_h(OBC%remap_CS, & segment%field(m)%nk_src, scl_fac* segment%field(m)%dz_src(i,J,:), & segment%field(m)%buffer_src(i,J,:), & - GV%ke, dz(i,j+jshift,:), segment%field(m)%buffer_dst(i,J,:), & - dz_neglect, dz_neglect_edge) + GV%ke, dz(i,j+jshift,:), segment%field(m)%buffer_dst(i,J,:)) endif enddo endif @@ -5567,10 +5579,9 @@ subroutine remap_OBC_fields(G, GV, h_old, h_new, OBC, PCM_cell) if (present(PCM_cell)) then call remapping_core_h(OBC%remap_CS, nz, h1, segment%tr_Reg%Tr(m)%tres(I,j,:), nz, h2, tr_column, & - h_neglect, h_neglect, PCM_cell=PCM) + PCM_cell=PCM) else - call remapping_core_h(OBC%remap_CS, nz, h1, segment%tr_Reg%Tr(m)%tres(I,j,:), nz, h2, tr_column, & - h_neglect, h_neglect) + call remapping_core_h(OBC%remap_CS, nz, h1, segment%tr_Reg%Tr(m)%tres(I,j,:), nz, h2, tr_column) endif ! Possibly underflow any very tiny tracer concentrations to 0? @@ -5585,7 +5596,7 @@ subroutine remap_OBC_fields(G, GV, h_old, h_new, OBC, PCM_cell) if (segment%radiation .and. (OBC%gamma_uv < 1.0)) then call remapping_core_h(OBC%remap_CS, nz, h1, segment%rx_norm_rad(I,j,:), nz, h2, r_norm_col, & - h_neglect, h_neglect, PCM_cell=PCM) + PCM_cell=PCM) do k=1,nz segment%rx_norm_rad(I,j,k) = r_norm_col(k) @@ -5595,13 +5606,13 @@ subroutine remap_OBC_fields(G, GV, h_old, h_new, OBC, PCM_cell) if (segment%oblique .and. (OBC%gamma_uv < 1.0)) then call remapping_core_h(OBC%remap_CS, nz, h1, segment%rx_norm_obl(I,j,:), nz, h2, rxy_col, & - h_neglect, h_neglect, PCM_cell=PCM) + PCM_cell=PCM) segment%rx_norm_obl(I,j,:) = rxy_col(:) call remapping_core_h(OBC%remap_CS, nz, h1, segment%ry_norm_obl(I,j,:), nz, h2, rxy_col, & - h_neglect, h_neglect, PCM_cell=PCM) + PCM_cell=PCM) segment%ry_norm_obl(I,j,:) = rxy_col(:) call remapping_core_h(OBC%remap_CS, nz, h1, segment%cff_normal(I,j,:), nz, h2, rxy_col, & - h_neglect, h_neglect, PCM_cell=PCM) + PCM_cell=PCM) segment%cff_normal(I,j,:) = rxy_col(:) do k=1,nz @@ -5635,10 +5646,9 @@ subroutine remap_OBC_fields(G, GV, h_old, h_new, OBC, PCM_cell) if (present(PCM_cell)) then call remapping_core_h(OBC%remap_CS, nz, h1, segment%tr_Reg%Tr(m)%tres(i,J,:), nz, h2, tr_column, & - h_neglect, h_neglect, PCM_cell=PCM) + PCM_cell=PCM) else - call remapping_core_h(OBC%remap_CS, nz, h1, segment%tr_Reg%Tr(m)%tres(i,J,:), nz, h2, tr_column, & - h_neglect, h_neglect) + call remapping_core_h(OBC%remap_CS, nz, h1, segment%tr_Reg%Tr(m)%tres(i,J,:), nz, h2, tr_column) endif ! Possibly underflow any very tiny tracer concentrations to 0? @@ -5653,7 +5663,7 @@ subroutine remap_OBC_fields(G, GV, h_old, h_new, OBC, PCM_cell) if (segment%radiation .and. (OBC%gamma_uv < 1.0)) then call remapping_core_h(OBC%remap_CS, nz, h1, segment%ry_norm_rad(i,J,:), nz, h2, r_norm_col, & - h_neglect, h_neglect, PCM_cell=PCM) + PCM_cell=PCM) do k=1,nz segment%ry_norm_rad(i,J,k) = r_norm_col(k) @@ -5663,13 +5673,13 @@ subroutine remap_OBC_fields(G, GV, h_old, h_new, OBC, PCM_cell) if (segment%oblique .and. (OBC%gamma_uv < 1.0)) then call remapping_core_h(OBC%remap_CS, nz, h1, segment%rx_norm_obl(i,J,:), nz, h2, rxy_col, & - h_neglect, h_neglect, PCM_cell=PCM) + PCM_cell=PCM) segment%rx_norm_obl(i,J,:) = rxy_col(:) call remapping_core_h(OBC%remap_CS, nz, h1, segment%ry_norm_obl(i,J,:), nz, h2, rxy_col, & - h_neglect, h_neglect, PCM_cell=PCM) + PCM_cell=PCM) segment%ry_norm_obl(i,J,:) = rxy_col(:) call remapping_core_h(OBC%remap_CS, nz, h1, segment%cff_normal(i,J,:), nz, h2, rxy_col, & - h_neglect, h_neglect, PCM_cell=PCM) + PCM_cell=PCM) segment%cff_normal(i,J,:) = rxy_col(:) do k=1,nz diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90 index b819c39ef1..cb01cc030c 100644 --- a/src/diagnostics/MOM_diagnostics.F90 +++ b/src/diagnostics/MOM_diagnostics.F90 @@ -1859,7 +1859,7 @@ subroutine MOM_diagnostics_init(MIS, ADp, CDp, Time, G, GV, US, param_file, diag if ((CS%id_cg1>0) .or. (CS%id_Rd1>0) .or. (CS%id_cfl_cg1>0) .or. & (CS%id_cfl_cg1_x>0) .or. (CS%id_cfl_cg1_y>0) .or. & (CS%id_cg_ebt>0) .or. (CS%id_Rd_ebt>0) .or. (CS%id_p_ebt>0)) then - call wave_speed_init(CS%wave_speed, remap_answer_date=remap_answer_date, & + call wave_speed_init(CS%wave_speed, GV, remap_answer_date=remap_answer_date, & better_speed_est=better_speed_est, min_speed=wave_speed_min, & wave_speed_tol=wave_speed_tol, om4_remap_via_sub_cells=om4_remap_via_sub_cells) endif diff --git a/src/diagnostics/MOM_wave_speed.F90 b/src/diagnostics/MOM_wave_speed.F90 index 8ee271f315..b520fafc1a 100644 --- a/src/diagnostics/MOM_wave_speed.F90 +++ b/src/diagnostics/MOM_wave_speed.F90 @@ -673,15 +673,8 @@ subroutine wave_speed(h, tv, G, GV, US, cg1, CS, halo_size, use_ebt_mode, mono_N mode_struct(1:kc)=0. endif - if (CS%remap_answer_date < 20190101) then - call remapping_core_h(CS%remapping_CS, kc, Hc(:), mode_struct, & - nz, h(i,j,:), modal_structure(i,j,:), & - 1.0e-30*GV%m_to_H, 1.0e-10*GV%m_to_H) - else - call remapping_core_h(CS%remapping_CS, kc, Hc(:), mode_struct, & - nz, h(i,j,:), modal_structure(i,j,:), & - GV%H_subroundoff, GV%H_subroundoff) - endif + call remapping_core_h(CS%remapping_CS, kc, Hc(:), mode_struct, & + nz, h(i,j,:), modal_structure(i,j,:)) endif else cg1(i,j) = 0.0 @@ -1358,8 +1351,7 @@ subroutine wave_speeds(h, tv, G, GV, US, nmodes, cn, CS, w_struct, u_struct, u_s ! for u (remap) onto all layers call remapping_core_h(CS%remapping_CS, kc, Hc(1:kc), mode_struct_fder(1:kc), & - nz, h(i,j,:), modal_structure_fder(:), & - GV%H_subroundoff, GV%H_subroundoff) + nz, h(i,j,:), modal_structure_fder(:)) ! write the wave structure do k=1,nz+1 @@ -1534,8 +1526,7 @@ subroutine wave_speeds(h, tv, G, GV, US, nmodes, cn, CS, w_struct, u_struct, u_s ! for u (remap) onto all layers call remapping_core_h(CS%remapping_CS, kc, Hc(1:kc), mode_struct_fder(1:kc), & - nz, h(i,j,:), modal_structure_fder(:), & - GV%H_subroundoff, GV%H_subroundoff) + nz, h(i,j,:), modal_structure_fder(:)) ! write the wave structure ! note that m=1 solves for 2nd mode,... @@ -1610,10 +1601,11 @@ subroutine tridiag_det(a, c, ks, ke, lam, det, ddet, row_scale) end subroutine tridiag_det !> Initialize control structure for MOM_wave_speed -subroutine wave_speed_init(CS, use_ebt_mode, mono_N2_column_fraction, mono_N2_depth, remap_answers_2018, & +subroutine wave_speed_init(CS, GV, use_ebt_mode, mono_N2_column_fraction, mono_N2_depth, remap_answers_2018, & remap_answer_date, better_speed_est, om4_remap_via_sub_cells, & min_speed, wave_speed_tol, c1_thresh) type(wave_speed_CS), intent(inout) :: CS !< Wave speed control struct + type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure logical, optional, intent(in) :: use_ebt_mode !< If true, use the equivalent !! barotropic mode instead of the first baroclinic mode. real, optional, intent(in) :: mono_N2_column_fraction !< The lower fraction of water column over @@ -1658,9 +1650,17 @@ subroutine wave_speed_init(CS, use_ebt_mode, mono_N2_column_fraction, mono_N2_de c1_thresh=c1_thresh) ! The remap_answers_2018 argument here is irrelevant, because remapping is hard-coded to use PLM. - call initialize_remapping(CS%remapping_CS, 'PLM', boundary_extrapolation=.false., & - om4_remap_via_sub_cells=om4_remap_via_sub_cells, & - answer_date=CS%remap_answer_date) + if (CS%remap_answer_date < 20190101) then + call initialize_remapping(CS%remapping_CS, 'PLM', boundary_extrapolation=.false., & + om4_remap_via_sub_cells=om4_remap_via_sub_cells, & + answer_date=CS%remap_answer_date, & + h_neglect=1.0e-30*GV%m_to_H, h_neglect_edge=1.0e-10*GV%m_to_H) + else + call initialize_remapping(CS%remapping_CS, 'PLM', boundary_extrapolation=.false., & + om4_remap_via_sub_cells=om4_remap_via_sub_cells, & + answer_date=CS%remap_answer_date, & + h_neglect=GV%H_subroundoff, h_neglect_edge=GV%H_subroundoff) + endif end subroutine wave_speed_init diff --git a/src/framework/MOM_diag_remap.F90 b/src/framework/MOM_diag_remap.F90 index e8e6a756e9..1151cd04b2 100644 --- a/src/framework/MOM_diag_remap.F90 +++ b/src/framework/MOM_diag_remap.F90 @@ -313,7 +313,8 @@ subroutine diag_remap_update(remap_cs, G, GV, US, h, T, S, eqn_of_state, h_targe ! Initialize remapping and regridding on the first call call initialize_remapping(remap_cs%remap_cs, 'PPM_IH4', boundary_extrapolation=.false., & om4_remap_via_sub_cells=remap_cs%om4_remap_via_sub_cells, & - answer_date=remap_cs%answer_date) + answer_date=remap_cs%answer_date, & + h_neglect=h_neglect, h_neglect_edge=h_neglect_edge) remap_cs%initialized = .true. endif @@ -432,16 +433,9 @@ subroutine do_remap(remap_cs, G, GV, US, isdf, jsdf, h, staggered_in_x, staggere ! Local variables real, dimension(remap_cs%nz) :: h_dest ! Destination thicknesses [H ~> m or kg m-2] or [Z ~> m] real, dimension(size(h,3)) :: h_src ! A column of source thicknesses [H ~> m or kg m-2] or [Z ~> m] - real :: h_neglect, h_neglect_edge ! Negligible thicknesses [H ~> m or kg m-2] or [Z ~> m] integer :: nz_src, nz_dest ! The number of layers on the native and remapped grids integer :: i, j ! Grid index - if (remap_cs%Z_based_coord) then - h_neglect = set_dz_neglect(GV, US, remap_cs%answer_date, h_neglect_edge) - else - h_neglect = set_h_neglect(GV, remap_cs%answer_date, h_neglect_edge) - endif - nz_src = size(field,3) nz_dest = remap_cs%nz remapped_field(:,:,:) = 0. @@ -453,14 +447,14 @@ subroutine do_remap(remap_cs, G, GV, US, isdf, jsdf, h, staggered_in_x, staggere h_src(:) = 0.5 * (h(i,j,:) + h(i+1,j,:)) h_dest(:) = 0.5 * (remap_cs%h(i,j,:) + remap_cs%h(i+1,j,:)) call remapping_core_h(remap_cs%remap_cs, nz_src, h_src(:), field(I,j,:), & - nz_dest, h_dest(:), remapped_field(I,j,:), h_neglect, h_neglect_edge) + nz_dest, h_dest(:), remapped_field(I,j,:)) endif ; enddo ; enddo else do j=G%jsc,G%jec ; do I=G%IscB,G%IecB h_src(:) = 0.5 * (h(i,j,:) + h(i+1,j,:)) h_dest(:) = 0.5 * (remap_cs%h(i,j,:) + remap_cs%h(i+1,j,:)) call remapping_core_h(remap_cs%remap_cs, nz_src, h_src(:), field(I,j,:), & - nz_dest, h_dest(:), remapped_field(I,j,:), h_neglect, h_neglect_edge) + nz_dest, h_dest(:), remapped_field(I,j,:)) enddo ; enddo endif elseif (staggered_in_y .and. .not. staggered_in_x) then @@ -470,14 +464,14 @@ subroutine do_remap(remap_cs, G, GV, US, isdf, jsdf, h, staggered_in_x, staggere h_src(:) = 0.5 * (h(i,j,:) + h(i,j+1,:)) h_dest(:) = 0.5 * (remap_cs%h(i,j,:) + remap_cs%h(i,j+1,:)) call remapping_core_h(remap_cs%remap_cs, nz_src, h_src(:), field(i,J,:), & - nz_dest, h_dest(:), remapped_field(i,J,:), h_neglect, h_neglect_edge) + nz_dest, h_dest(:), remapped_field(i,J,:)) endif ; enddo ; enddo else do J=G%jscB,G%jecB ; do i=G%isc,G%iec h_src(:) = 0.5 * (h(i,j,:) + h(i,j+1,:)) h_dest(:) = 0.5 * (remap_cs%h(i,j,:) + remap_cs%h(i,j+1,:)) call remapping_core_h(remap_cs%remap_cs, nz_src, h_src(:), field(i,J,:), & - nz_dest, h_dest(:), remapped_field(i,J,:), h_neglect, h_neglect_edge) + nz_dest, h_dest(:), remapped_field(i,J,:)) enddo ; enddo endif elseif ((.not. staggered_in_x) .and. (.not. staggered_in_y)) then @@ -485,14 +479,12 @@ subroutine do_remap(remap_cs, G, GV, US, isdf, jsdf, h, staggered_in_x, staggere if (present(mask)) then do j=G%jsc,G%jec ; do i=G%isc,G%iec ; if (mask(i,j) > 0.) then call remapping_core_h(remap_cs%remap_cs, nz_src, h(i,j,:), field(i,j,:), & - nz_dest, remap_cs%h(i,j,:), remapped_field(i,j,:), & - h_neglect, h_neglect_edge) + nz_dest, remap_cs%h(i,j,:), remapped_field(i,j,:)) endif ; enddo ; enddo else do j=G%jsc,G%jec ; do i=G%isc,G%iec call remapping_core_h(remap_cs%remap_cs, nz_src, h(i,j,:), field(i,j,:), & - nz_dest, remap_cs%h(i,j,:), remapped_field(i,j,:), & - h_neglect, h_neglect_edge) + nz_dest, remap_cs%h(i,j,:), remapped_field(i,j,:)) enddo ; enddo endif else diff --git a/src/initialization/MOM_fixed_initialization.F90 b/src/initialization/MOM_fixed_initialization.F90 index 10236994e0..e8ae4ac0f4 100644 --- a/src/initialization/MOM_fixed_initialization.F90 +++ b/src/initialization/MOM_fixed_initialization.F90 @@ -26,6 +26,7 @@ module MOM_fixed_initialization use MOM_shared_initialization, only : set_subgrid_topo_at_vel_from_file use MOM_shared_initialization, only : compute_global_grid_integrals, write_ocean_geometry_file use MOM_unit_scaling, only : unit_scale_type +use MOM_verticalGrid, only : verticalGrid_type use user_initialization, only : user_initialize_topography use DOME_initialization, only : DOME_initialize_topography diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index 1467cdaaad..c18a37f277 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -1189,7 +1189,13 @@ subroutine trim_for_ice(PF, G, GV, US, ALE_CSp, tv, h, just_read) if (use_remapping) then allocate(remap_CS) - call initialize_remapping(remap_CS, 'PLM', boundary_extrapolation=.true.) + if (remap_answer_date < 20190101) then + call initialize_remapping(remap_CS, 'PLM', boundary_extrapolation=.true., & + h_neglect=1.0e-30*GV%m_to_H, h_neglect_edge=1.0e-10*GV%m_to_H) + else + call initialize_remapping(remap_CS, 'PLM', boundary_extrapolation=.true., & + h_neglect=GV%H_subroundoff, h_neglect_edge=GV%H_subroundoff) + endif endif ! Find edge values of T and S used in reconstructions @@ -1332,11 +1338,8 @@ subroutine cut_off_column_top(nk, tv, GV, US, G_earth, depth, min_thickness, T, real :: z_out, e_top ! Interface height positions [Z ~> m] real :: min_dz ! The minimum thickness in depth units [Z ~> m] real :: dh_surf_rem ! The remaining thickness to remove in non-Bousinesq mode [H ~> kg m-2] - logical :: answers_2018 integer :: k - answers_2018 = .true. ; if (present(remap_answer_date)) answers_2018 = (remap_answer_date < 20190101) - ! Keep a copy of the initial thicknesses in reverse order to use in remapping do k=1,nk ; h0(k) = h(nk+1-k) ; enddo @@ -1407,13 +1410,8 @@ subroutine cut_off_column_top(nk, tv, GV, US, G_earth, depth, min_thickness, T, T0(k) = T(nk+1-k) h1(k) = h(nk+1-k) enddo - if (answers_2018) then - call remapping_core_h(remap_CS, nk, h0, T0, nk, h1, T1, 1.0e-30*GV%m_to_H, 1.0e-10*GV%m_to_H) - call remapping_core_h(remap_CS, nk, h0, S0, nk, h1, S1, 1.0e-30*GV%m_to_H, 1.0e-10*GV%m_to_H) - else - call remapping_core_h(remap_CS, nk, h0, T0, nk, h1, T1, GV%H_subroundoff, GV%H_subroundoff) - call remapping_core_h(remap_CS, nk, h0, S0, nk, h1, S1, GV%H_subroundoff, GV%H_subroundoff) - endif + call remapping_core_h(remap_CS, nk, h0, T0, nk, h1, T1) + call remapping_core_h(remap_CS, nk, h0, S0, nk, h1, S1) do k=1,nk S(k) = S1(nk+1-k) T(k) = T1(nk+1-k) @@ -2758,8 +2756,10 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just ! Build the target grid (and set the model thickness to it) call ALE_initRegridding( GV, US, G%max_depth, PF, mdl, regridCS ) ! sets regridCS + dz_neglect = set_dz_neglect(GV, US, remap_answer_date, dz_neglect_edge) call initialize_remapping( remapCS, remappingScheme, boundary_extrapolation=.false., & - om4_remap_via_sub_cells=om4_remap_via_sub_cells, answer_date=remap_answer_date ) + om4_remap_via_sub_cells=om4_remap_via_sub_cells, answer_date=remap_answer_date, & + h_neglect=dz_neglect, h_neglect_edge=dz_neglect_edge ) ! Now remap from source grid to target grid, first setting reconstruction parameters if (remap_general) then @@ -2797,13 +2797,10 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just enddo ; enddo deallocate( hTarget ) - dz_neglect = set_dz_neglect(GV, US, remap_answer_date, dz_neglect_edge) call ALE_remap_scalar(remapCS, G, GV, nkd, dz1, tmpT1dIn, dz, tv%T, all_cells=remap_full_column, & - old_remap=remap_old_alg, answer_date=remap_answer_date, & - H_neglect=dz_neglect, H_neglect_edge=dz_neglect_edge) + old_remap=remap_old_alg, answer_date=remap_answer_date) call ALE_remap_scalar(remapCS, G, GV, nkd, dz1, tmpS1dIn, dz, tv%S, all_cells=remap_full_column, & - old_remap=remap_old_alg, answer_date=remap_answer_date, & - H_neglect=dz_neglect, H_neglect_edge=dz_neglect_edge) + old_remap=remap_old_alg, answer_date=remap_answer_date) if (GV%Boussinesq .or. GV%semi_Boussinesq) then ! This is a simple conversion of the target grid to thickness units that is not diff --git a/src/initialization/MOM_tracer_initialization_from_Z.F90 b/src/initialization/MOM_tracer_initialization_from_Z.F90 index bafa5d8c36..57a60de8b8 100644 --- a/src/initialization/MOM_tracer_initialization_from_Z.F90 +++ b/src/initialization/MOM_tracer_initialization_from_Z.F90 @@ -179,8 +179,10 @@ subroutine MOM_initialize_tracer_from_Z(h, tr, G, GV, US, PF, src_file, src_var_ allocate( dzSrc(isd:ied,jsd:jed,kd) ) allocate( hSrc(isd:ied,jsd:jed,kd) ) ! Set parameters for reconstructions + dz_neglect = set_dz_neglect(GV, US, remap_answer_date, dz_neglect_edge) call initialize_remapping( remapCS, remapScheme, boundary_extrapolation=.false., & - om4_remap_via_sub_cells=om4_remap_via_sub_cells, answer_date=remap_answer_date ) + om4_remap_via_sub_cells=om4_remap_via_sub_cells, answer_date=remap_answer_date, & + h_neglect=dz_neglect, h_neglect_edge=dz_neglect_edge ) ! Next we initialize the regridding package so that it knows about the target grid do j = js, je ; do i = is, ie @@ -208,9 +210,7 @@ subroutine MOM_initialize_tracer_from_Z(h, tr, G, GV, US, PF, src_file, src_var_ if (h_is_in_Z_units) then ! Because h is in units of [Z ~> m], dzSrc is already in the right units, but we need to ! specify negligible thickness values with the right units. - dz_neglect = set_dz_neglect(GV, US, remap_answer_date, dz_neglect_edge) - call ALE_remap_scalar(remapCS, G, GV, kd, dzSrc, tr_z, h, tr, all_cells=.false., answer_date=remap_answer_date, & - H_neglect=dz_neglect, H_neglect_edge=dz_neglect_edge) + call ALE_remap_scalar(remapCS, G, GV, kd, dzSrc, tr_z, h, tr, all_cells=.false., answer_date=remap_answer_date) else ! Equation of state data is not available, so a simpler rescaling will have to suffice, ! but it might be problematic in non-Boussinesq mode. diff --git a/src/ocean_data_assim/MOM_oda_driver.F90 b/src/ocean_data_assim/MOM_oda_driver.F90 index 9275555afc..832fcdd203 100644 --- a/src/ocean_data_assim/MOM_oda_driver.F90 +++ b/src/ocean_data_assim/MOM_oda_driver.F90 @@ -184,6 +184,7 @@ subroutine init_oda(Time, G, GV, US, diag_CS, CS) character(len=80) :: bias_correction_file, inc_file integer :: default_answer_date ! The default setting for the various ANSWER_DATE flags. logical :: om4_remap_via_sub_cells ! If true, use the OM4 remapping algorithm + real :: h_neglect, h_neglect_edge ! small thicknesses [H ~> m or kg m-2] if (associated(CS)) call MOM_error(FATAL, 'Calling oda_init with associated control structure') allocate(CS) @@ -323,6 +324,14 @@ subroutine init_oda(Time, G, GV, US, diag_CS, CS) default="ZSTAR", fail_if_missing=.false.) call get_param(PF, mdl, "REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, & do_not_log=.true., default=.true.) + if (CS%answer_date >= 20190101) then + h_neglect = GV%H_subroundoff ; h_neglect_edge = GV%H_subroundoff + elseif (GV%Boussinesq) then + h_neglect = GV%m_to_H * 1.0e-30 ; h_neglect_edge = GV%m_to_H * 1.0e-10 + else + h_neglect = GV%kg_m2_to_H * 1.0e-30 ; h_neglect_edge = GV%kg_m2_to_H * 1.0e-10 + endif + call initialize_regridding(CS%regridCS, CS%GV, CS%US, dG%max_depth,PF,'oda_driver',coord_mode,'','') call initialize_remapping(CS%remapCS, remap_scheme, om4_remap_via_sub_cells=om4_remap_via_sub_cells) call set_regrid_params(CS%regridCS, min_thickness=0.) @@ -415,7 +424,6 @@ subroutine set_prior_tracer(Time, G, GV, h, tv, CS) real, dimension(SZI_(G),SZJ_(G),CS%nk) :: S ! Salinity on the analysis grid [S ~> ppt] integer :: i, j, m integer :: isc, iec, jsc, jec - real :: h_neglect, h_neglect_edge ! small thicknesses [H ~> m or kg m-2] ! return if not time for analysis if (Time < CS%Time) return @@ -427,14 +435,6 @@ subroutine set_prior_tracer(Time, G, GV, h, tv, CS) call set_PElist(CS%filter_pelist) !call MOM_mesg('Setting prior') - if (CS%answer_date >= 20190101) then - h_neglect = GV%H_subroundoff ; h_neglect_edge = GV%H_subroundoff - elseif (GV%Boussinesq) then - h_neglect = GV%m_to_H * 1.0e-30 ; h_neglect_edge = GV%m_to_H * 1.0e-10 - else - h_neglect = GV%kg_m2_to_H * 1.0e-30 ; h_neglect_edge = GV%kg_m2_to_H * 1.0e-10 - endif - ! computational domain for the analysis grid isc=CS%Grid%isc;iec=CS%Grid%iec;jsc=CS%Grid%jsc;jec=CS%Grid%jec ! array extents for the ensemble member @@ -443,9 +443,9 @@ subroutine set_prior_tracer(Time, G, GV, h, tv, CS) ! remap temperature and salinity from the ensemble member to the analysis grid do j=G%jsc,G%jec ; do i=G%isc,G%iec call remapping_core_h(CS%remapCS, GV%ke, h(i,j,:), tv%T(i,j,:), & - CS%nk, CS%h(i,j,:), T(i,j,:), h_neglect, h_neglect_edge) + CS%nk, CS%h(i,j,:), T(i,j,:)) call remapping_core_h(CS%remapCS, GV%ke, h(i,j,:), tv%S(i,j,:), & - CS%nk, CS%h(i,j,:), S(i,j,:), h_neglect, h_neglect_edge) + CS%nk, CS%h(i,j,:), S(i,j,:)) enddo ; enddo ! cast ensemble members to the analysis domain do m=1,CS%ensemble_size @@ -711,9 +711,9 @@ subroutine apply_oda_tracer_increments(dt, Time_end, G, GV, tv, h, CS) isc=G%isc; iec=G%iec; jsc=G%jsc; jec=G%jec do j=jsc,jec; do i=isc,iec call remapping_core_h(CS%remapCS, CS%nk, CS%h(i,j,:), T_tend(i,j,:), & - G%ke, h(i,j,:), T_tend_inc(i,j,:), h_neglect, h_neglect_edge) + G%ke, h(i,j,:), T_tend_inc(i,j,:)) call remapping_core_h(CS%remapCS, CS%nk, CS%h(i,j,:), S_tend(i,j,:), & - G%ke, h(i,j,:), S_tend_inc(i,j,:), h_neglect, h_neglect_edge) + G%ke, h(i,j,:), S_tend_inc(i,j,:)) enddo; enddo diff --git a/src/ocean_data_assim/MOM_oda_incupd.F90 b/src/ocean_data_assim/MOM_oda_incupd.F90 index 94d09554c2..cb57139824 100644 --- a/src/ocean_data_assim/MOM_oda_incupd.F90 +++ b/src/ocean_data_assim/MOM_oda_incupd.F90 @@ -144,6 +144,7 @@ subroutine initialize_oda_incupd( G, GV, US, param_file, CS, data_h, nz_data, re character(len=256) :: mesg character(len=64) :: remapScheme logical :: om4_remap_via_sub_cells ! If true, use the OM4 remapping algorithm + real :: h_neglect, h_neglect_edge ! Negligible thicknesses [H ~> m or kg m-2] if (.not.associated(CS)) then call MOM_error(WARNING, "initialize_oda_incupd called without an associated "// & @@ -239,8 +240,15 @@ subroutine initialize_oda_incupd( G, GV, US, param_file, CS, data_h, nz_data, re ! Call the constructor for remapping control structure !### Revisit this hard-coded answer_date. + if (GV%Boussinesq) then + h_neglect = GV%m_to_H*1.0e-30 ; h_neglect_edge = GV%m_to_H*1.0e-10 + else + h_neglect = GV%H_subroundoff ; h_neglect_edge = GV%H_subroundoff + endif + call initialize_remapping(CS%remap_cs, remapScheme, boundary_extrapolation=bndExtrapolation, & - om4_remap_via_sub_cells=om4_remap_via_sub_cells, answer_date=20190101) + om4_remap_via_sub_cells=om4_remap_via_sub_cells, answer_date=20190101, & + h_neglect=h_neglect, h_neglect_edge=h_neglect_edge) end subroutine initialize_oda_incupd @@ -347,7 +355,6 @@ subroutine calc_oda_increments(h, tv, u, v, G, GV, US, CS) integer :: i, j, k, is, ie, js, je, nz, nz_data integer :: isB, ieB, jsB, jeB - real :: h_neglect, h_neglect_edge ! Negligible thicknesses [H ~> m or kg m-2] real :: sum_h1, sum_h2 ! vertical sums of h's [H ~> m or kg m-2] is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke @@ -359,13 +366,6 @@ subroutine calc_oda_increments(h, tv, u, v, G, GV, US, CS) if (CS%ncount /= 0.0) call MOM_error(FATAL,'calc_oda_increments: '// & 'CS%ncount should be 0.0 to get accurate increments.') - - if (GV%Boussinesq) then - h_neglect = GV%m_to_H*1.0e-30 ; h_neglect_edge = GV%m_to_H*1.0e-10 - else - h_neglect = GV%H_subroundoff ; h_neglect_edge = GV%H_subroundoff - endif - ! get h_obs nz_data = CS%Inc(1)%nz_data allocate(h_obs(G%isd:G%ied,G%jsd:G%jed,nz_data), source=0.0) @@ -404,8 +404,7 @@ subroutine calc_oda_increments(h, tv, u, v, G, GV, US, CS) enddo ! remap tracer on h_obs call remapping_core_h(CS%remap_cs, nz, h(i,j,1:nz), tmp_val1, & - nz_data, tmp_h(1:nz_data), tmp_val2, & - h_neglect, h_neglect_edge) + nz_data, tmp_h(1:nz_data), tmp_val2) ! get increment from full field on h_obs do k=1,nz_data CS%Inc(1)%p(i,j,k) = CS%Inc(1)%p(i,j,k) - tmp_val2(k) @@ -417,8 +416,7 @@ subroutine calc_oda_increments(h, tv, u, v, G, GV, US, CS) enddo ! remap tracer on h_obs call remapping_core_h(CS%remap_cs, nz, h(i,j,1:nz), tmp_val1, & - nz_data, tmp_h(1:nz_data), tmp_val2, & - h_neglect, h_neglect_edge) + nz_data, tmp_h(1:nz_data), tmp_val2) ! get increment from full field on h_obs do k=1,nz_data CS%Inc(2)%p(i,j,k) = CS%Inc(2)%p(i,j,k) - tmp_val2(k) @@ -456,8 +454,7 @@ subroutine calc_oda_increments(h, tv, u, v, G, GV, US, CS) enddo ! remap model u on hu_obs call remapping_core_h(CS%remap_cs, nz, hu(1:nz), tmp_val1, & - nz_data, hu_obs(1:nz_data), tmp_val2, & - h_neglect, h_neglect_edge) + nz_data, hu_obs(1:nz_data), tmp_val2) ! get increment from full field on h_obs do k=1,nz_data CS%Inc_u%p(i,j,k) = CS%Inc_u%p(i,j,k) - tmp_val2(k) @@ -492,8 +489,7 @@ subroutine calc_oda_increments(h, tv, u, v, G, GV, US, CS) enddo ! remap model v on hv_obs call remapping_core_h(CS%remap_cs, nz, hv(1:nz), tmp_val1, & - nz_data, hv_obs(1:nz_data), tmp_val2, & - h_neglect, h_neglect_edge) + nz_data, hv_obs(1:nz_data), tmp_val2) ! get increment from full field on h_obs do k=1,nz_data CS%Inc_v%p(i,j,k) = CS%Inc_v%p(i,j,k) - tmp_val2(k) @@ -621,7 +617,7 @@ subroutine apply_oda_incupd(h, tv, u, v, dt, G, GV, US, CS) enddo ! remap increment profile on model h call remapping_core_h(CS%remap_cs, nz_data, tmp_h(1:nz_data), tmp_val2, & - nz, h(i,j,1:nz),tmp_val1, h_neglect, h_neglect_edge) + nz, h(i,j,1:nz), tmp_val1) do k=1,nz ! add increment to tracer on model h tv%T(i,j,k) = tv%T(i,j,k) + inc_wt * tmp_val1(k) @@ -634,7 +630,7 @@ subroutine apply_oda_incupd(h, tv, u, v, dt, G, GV, US, CS) enddo ! remap increment profile on model h call remapping_core_h(CS%remap_cs, nz_data, tmp_h(1:nz_data),tmp_val2,& - nz, h(i,j,1:nz),tmp_val1, h_neglect, h_neglect_edge) + nz, h(i,j,1:nz), tmp_val1) ! add increment to tracer on model h do k=1,nz tv%S(i,j,k) = tv%S(i,j,k) + inc_wt * tmp_val1(k) @@ -680,7 +676,7 @@ subroutine apply_oda_incupd(h, tv, u, v, dt, G, GV, US, CS) enddo ! remap increment profile on hu call remapping_core_h(CS%remap_cs, nz_data, hu_obs(1:nz_data), tmp_val2, & - nz, hu(1:nz), tmp_val1, h_neglect, h_neglect_edge) + nz, hu(1:nz), tmp_val1) ! add increment to u-velocity on hu do k=1,nz u(i,j,k) = u(i,j,k) + inc_wt * tmp_val1(k) @@ -718,7 +714,7 @@ subroutine apply_oda_incupd(h, tv, u, v, dt, G, GV, US, CS) enddo ! remap increment profile on hv call remapping_core_h(CS%remap_cs, nz_data, hv_obs(1:nz_data), tmp_val2, & - nz, hv(1:nz), tmp_val1, h_neglect, h_neglect_edge) + nz, hv(1:nz), tmp_val1) ! add increment to v-velocity on hv do k=1,nz v(i,j,k) = v(i,j,k) + inc_wt * tmp_val1(k) diff --git a/src/parameterizations/lateral/MOM_internal_tides.F90 b/src/parameterizations/lateral/MOM_internal_tides.F90 index 8b3a2c2dbe..97f9a0ba6d 100644 --- a/src/parameterizations/lateral/MOM_internal_tides.F90 +++ b/src/parameterizations/lateral/MOM_internal_tides.F90 @@ -3112,7 +3112,7 @@ subroutine internal_tides_init(Time, G, GV, US, param_file, diag, CS) enddo ! Initialize the module that calculates the wave speeds. - call wave_speed_init(CS%wave_speed, c1_thresh=IGW_c1_thresh, & + call wave_speed_init(CS%wave_speed, GV, c1_thresh=IGW_c1_thresh, & om4_remap_via_sub_cells=om4_remap_via_sub_cells) end subroutine internal_tides_init diff --git a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 index eb1cf7b1cf..4d2d4b44ff 100644 --- a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 +++ b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 @@ -1598,7 +1598,7 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS) "If true, use the OM4 remapping-via-subcells algorithm for calculating EBT structure. "//& "See REMAPPING_USE_OM4_SUBCELLS for details. "//& "We recommend setting this option to false.", default=.true.) - call wave_speed_init(CS%wave_speed, use_ebt_mode=CS%Resoln_use_ebt, & + call wave_speed_init(CS%wave_speed, GV, use_ebt_mode=CS%Resoln_use_ebt, & mono_N2_depth=N2_filter_depth, remap_answer_date=remap_answer_date, & better_speed_est=better_speed_est, min_speed=wave_speed_min, & om4_remap_via_sub_cells=om4_remap_via_sub_cells, wave_speed_tol=wave_speed_tol) diff --git a/src/parameterizations/vertical/MOM_ALE_sponge.F90 b/src/parameterizations/vertical/MOM_ALE_sponge.F90 index 170265d27a..b69e800254 100644 --- a/src/parameterizations/vertical/MOM_ALE_sponge.F90 +++ b/src/parameterizations/vertical/MOM_ALE_sponge.F90 @@ -468,6 +468,7 @@ subroutine initialize_ALE_sponge_varying(Iresttime, G, GV, US, param_file, CS, I character(len=40) :: mdl = "MOM_sponge" ! This module's name. real, allocatable, dimension(:,:) :: Iresttime_u !< inverse of the restoring time at u points [T-1 ~> s-1] real, allocatable, dimension(:,:) :: Iresttime_v !< inverse of the restoring time at v points [T-1 ~> s-1] + real :: dz_neglect, dz_neglect_edge ! Negligible layer extents [Z ~> m] ! This include declares and sets the variable "version". # include "version_variable.h" character(len=64) :: remapScheme @@ -559,9 +560,19 @@ subroutine initialize_ALE_sponge_varying(Iresttime, G, GV, US, param_file, CS, I call sum_across_PEs(total_sponge_cols) ! Call the constructor for remapping control structure + if (CS%remap_answer_date >= 20190101) then + dz_neglect = GV%dZ_subroundoff ; dz_neglect_edge = GV%dZ_subroundoff + elseif (GV%Boussinesq) then + dz_neglect = US%m_to_Z*1.0e-30 ; dz_neglect_edge = US%m_to_Z*1.0e-10 + elseif (GV%semi_Boussinesq) then + dz_neglect = GV%kg_m2_to_H*GV%H_to_Z*1.0e-30 ; dz_neglect_edge = GV%kg_m2_to_H*GV%H_to_Z*1.0e-10 + else + dz_neglect = GV%dZ_subroundoff ; dz_neglect_edge = GV%dZ_subroundoff + endif call initialize_remapping(CS%remap_cs, remapScheme, boundary_extrapolation=bndExtrapolation, & om4_remap_via_sub_cells=om4_remap_via_sub_cells, & - answer_date=CS%remap_answer_date) + answer_date=CS%remap_answer_date, & + h_neglect=dz_neglect, h_neglect_edge=dz_neglect_edge) call log_param(param_file, mdl, "!Total sponge columns at h points", total_sponge_cols, & "The total number of columns where sponges are applied at h points.", like_default=.true.) if (CS%sponge_uv) then @@ -950,7 +961,6 @@ subroutine apply_ALE_sponge(h, tv, dt, G, GV, US, CS, Time) ! edges in the input file [Z ~> m] real :: missing_value ! The missing value in the input data field [various] real :: Idt ! The inverse of the timestep [T-1 ~> s-1] - real :: dz_neglect, dz_neglect_edge ! Negligible layer extents [Z ~> m] real :: zTopOfCell, zBottomOfCell ! Interface heights (positive upward) in the input dataset [Z ~> m]. real :: sp_val_u ! Interpolation of sp_val to u-points, often a velocity in [L T-1 ~> m s-1] real :: sp_val_v ! Interpolation of sp_val to v-points, often a velocity in [L T-1 ~> m s-1] @@ -961,16 +971,6 @@ subroutine apply_ALE_sponge(h, tv, dt, G, GV, US, CS, Time) Idt = 1.0/dt - if (CS%remap_answer_date >= 20190101) then - dz_neglect = GV%dZ_subroundoff ; dz_neglect_edge = GV%dZ_subroundoff - elseif (GV%Boussinesq) then - dz_neglect = US%m_to_Z*1.0e-30 ; dz_neglect_edge = US%m_to_Z*1.0e-10 - elseif (GV%semi_Boussinesq) then - dz_neglect = GV%kg_m2_to_H*GV%H_to_Z*1.0e-30 ; dz_neglect_edge = GV%kg_m2_to_H*GV%H_to_Z*1.0e-10 - else - dz_neglect = GV%dZ_subroundoff ; dz_neglect_edge = GV%dZ_subroundoff - endif - if (CS%time_varying_sponges) then do m=1,CS%fldno nz_data = CS%Ref_val(m)%nz_data @@ -1038,12 +1038,11 @@ subroutine apply_ALE_sponge(h, tv, dt, G, GV, US, CS, Time) enddo endif if (CS%time_varying_sponges) then - call remapping_core_h(CS%remap_cs, nz_data, CS%Ref_val(m)%dz(1:nz_data,c), tmp_val2, & - CS%nz, dz_col, tmp_val1, dz_neglect, dz_neglect_edge) + CS%nz, dz_col, tmp_val1) else call remapping_core_h(CS%remap_cs, nz_data, CS%Ref_dz%p(1:nz_data,c), tmp_val2, & - CS%nz, dz_col, tmp_val1, dz_neglect, dz_neglect_edge) + CS%nz, dz_col, tmp_val1) endif !Backward Euler method if (CS%id_sp_tendency(m) > 0) tmp(i,j,1:nz) = CS%var(m)%p(i,j,1:nz) @@ -1189,10 +1188,10 @@ subroutine apply_ALE_sponge(h, tv, dt, G, GV, US, CS, Time) enddo if (CS%time_varying_sponges) then call remapping_core_h(CS%remap_cs, nz_data, CS%Ref_val_u%dz(1:nz_data,c), tmp_val2, & - CS%nz, dz_col, tmp_val1, dz_neglect, dz_neglect_edge) + CS%nz, dz_col, tmp_val1) else call remapping_core_h(CS%remap_cs, nz_data, CS%Ref_dzu%p(1:nz_data,c), tmp_val2, & - CS%nz, dz_col, tmp_val1, dz_neglect, dz_neglect_edge) + CS%nz, dz_col, tmp_val1) endif if (CS%id_sp_u_tendency > 0) tmp_u(i,j,1:nz) = CS%var_u%p(i,j,1:nz) !Backward Euler method @@ -1222,10 +1221,10 @@ subroutine apply_ALE_sponge(h, tv, dt, G, GV, US, CS, Time) enddo if (CS%time_varying_sponges) then call remapping_core_h(CS%remap_cs, nz_data, CS%Ref_val_v%dz(1:nz_data,c), tmp_val2, & - CS%nz, dz_col, tmp_val1, dz_neglect, dz_neglect_edge) + CS%nz, dz_col, tmp_val1) else call remapping_core_h(CS%remap_cs, nz_data, CS%Ref_dzv%p(1:nz_data,c), tmp_val2, & - CS%nz, dz_col, tmp_val1, dz_neglect, dz_neglect_edge) + CS%nz, dz_col, tmp_val1) endif if (CS%id_sp_v_tendency > 0) tmp_v(i,j,1:nz) = CS%var_v%p(i,j,1:nz) !Backward Euler method diff --git a/src/parameterizations/vertical/MOM_tidal_mixing.F90 b/src/parameterizations/vertical/MOM_tidal_mixing.F90 index 31f90cdcb1..ff53a15548 100644 --- a/src/parameterizations/vertical/MOM_tidal_mixing.F90 +++ b/src/parameterizations/vertical/MOM_tidal_mixing.F90 @@ -602,7 +602,7 @@ logical function tidal_mixing_init(Time, G, GV, US, param_file, int_tide_CSp, di local_mixing_frac = CS%Gamma_itides, & depth_cutoff = CS%min_zbot_itides*US%Z_to_m) - call read_tidal_energy(G, US, tidal_energy_type, param_file, CS) + call read_tidal_energy(G, GV, US, tidal_energy_type, param_file, CS) !call closeParameterBlock(param_file) @@ -912,7 +912,7 @@ subroutine calculate_CVMix_tidal(dz, j, N2_int, G, GV, US, CS, Kv, Kd_lay, Kd_in ! remap from input z coordinate to model coordinate: tidal_qe_md(:) = 0.0 call remapping_core_h(CS%remap_cs, size(CS%h_src), CS%h_src, CS%tidal_qe_3d_in(i,j,:), & - GV%ke, h_m, tidal_qe_md, GV%H_subroundoff, GV%H_subroundoff) + GV%ke, h_m, tidal_qe_md) ! form the Schmittner coefficient that is based on 3D q*E, which is formed from ! summing q_i*TidalConstituent_i over the number of constituents. @@ -1571,8 +1571,9 @@ end subroutine tidal_mixing_h_amp ! TODO: move this subroutine to MOM_internal_tide_input module (?) !> This subroutine read tidal energy inputs from a file. -subroutine read_tidal_energy(G, US, tidal_energy_type, param_file, CS) +subroutine read_tidal_energy(G, GV, US, tidal_energy_type, param_file, CS) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure. type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type character(len=20), intent(in) :: tidal_energy_type !< The type of tidal energy inputs to read type(param_file_type), intent(in) :: param_file !< Run-time parameter file handle @@ -1606,7 +1607,7 @@ subroutine read_tidal_energy(G, US, tidal_energy_type, param_file, CS) enddo ; enddo deallocate(tidal_energy_flux_2d) case ('ER03') ! Egbert & Ray 2003 - call read_tidal_constituents(G, US, tidal_energy_file, param_file, CS) + call read_tidal_constituents(G, GV, US, tidal_energy_file, param_file, CS) case default call MOM_error(FATAL, "read_tidal_energy: Unknown tidal energy file type.") end select @@ -1614,8 +1615,9 @@ subroutine read_tidal_energy(G, US, tidal_energy_type, param_file, CS) end subroutine read_tidal_energy !> This subroutine reads tidal input energy from a file by constituent. -subroutine read_tidal_constituents(G, US, tidal_energy_file, param_file, CS) +subroutine read_tidal_constituents(G, GV, US, tidal_energy_file, param_file, CS) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure. type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type character(len=200), intent(in) :: tidal_energy_file !< The file from which to read tidal energy inputs type(param_file_type), intent(in) :: param_file !< Run-time parameter file handle @@ -1700,7 +1702,8 @@ subroutine read_tidal_constituents(G, US, tidal_energy_file, param_file, CS) ! initialize input remapping: call initialize_remapping(CS%remap_cs, remapping_scheme="PLM", & boundary_extrapolation=.false., check_remapping=CS%debug, & - answer_date=CS%remap_answer_date) + answer_date=CS%remap_answer_date, & + h_neglect=GV%H_subroundoff, h_neglect_edge=GV%H_subroundoff) deallocate(tc_m2) deallocate(tc_s2) diff --git a/src/tracer/MOM_hor_bnd_diffusion.F90 b/src/tracer/MOM_hor_bnd_diffusion.F90 index b6714148ea..1d2a50d2c3 100644 --- a/src/tracer/MOM_hor_bnd_diffusion.F90 +++ b/src/tracer/MOM_hor_bnd_diffusion.F90 @@ -151,7 +151,8 @@ logical function hor_bnd_diffusion_init(Time, G, GV, US, param_file, diag, diaba ! GMM, TODO: add HBD params to control optional arguments in initialize_remapping. call initialize_remapping( CS%remap_CS, string, boundary_extrapolation=boundary_extrap, & om4_remap_via_sub_cells=om4_remap_via_sub_cells, & - check_reconstruction=.false., check_remapping=.false.) + check_reconstruction=.false., check_remapping=.false., & + h_neglect=CS%H_subroundoff, h_neglect_edge=CS%H_subroundoff) call extract_member_remapping_CS(CS%remap_CS, degree=CS%deg) call get_param(param_file, mdl, "DEBUG", debug, default=.false., do_not_log=.true.) call get_param(param_file, mdl, "HBD_DEBUG", CS%debug, & @@ -739,10 +740,8 @@ subroutine fluxes_layer_method(boundary, ke, hbl_L, hbl_R, h_L, h_R, phi_L, phi_ allocate(khtr_ul_z(nk), source=0.0) ! remap tracer to dz_top - call remapping_core_h(CS%remap_cs, ke, h_L(:), phi_L(:), nk, dz_top(:), phi_L_z(:), & - CS%H_subroundoff, CS%H_subroundoff) - call remapping_core_h(CS%remap_cs, ke, h_R(:), phi_R(:), nk, dz_top(:), phi_R_z(:), & - CS%H_subroundoff, CS%H_subroundoff) + call remapping_core_h(CS%remap_cs, ke, h_L(:), phi_L(:), nk, dz_top(:), phi_L_z(:)) + call remapping_core_h(CS%remap_cs, ke, h_R(:), phi_R(:), nk, dz_top(:), phi_R_z(:)) ! thicknesses at velocity points & khtr_u at layer centers do k = 1,ke @@ -753,8 +752,7 @@ subroutine fluxes_layer_method(boundary, ke, hbl_L, hbl_R, h_L, h_R, phi_L, phi_ enddo ! remap khtr_ul to khtr_ul_z - call remapping_core_h(CS%remap_cs, ke, h_vel(:), khtr_ul(:), nk, dz_top(:), khtr_ul_z(:), & - CS%H_subroundoff, CS%H_subroundoff) + call remapping_core_h(CS%remap_cs, ke, h_vel(:), khtr_ul(:), nk, dz_top(:), khtr_ul_z(:)) ! Calculate vertical indices containing the boundary layer in dz_top call boundary_k_range(boundary, nk, dz_top, hbl_L, k_top_L, zeta_top_L, k_bot_L, zeta_bot_L) @@ -855,15 +853,16 @@ logical function near_boundary_unit_tests( verbose ) allocate(CS) ! fill required fields in CS CS%linear=.false. - call initialize_remapping( CS%remap_CS, 'PLM', boundary_extrapolation=.true., & - om4_remap_via_sub_cells=.true., & ! ### see fail below when using fixed remapping alg. - check_reconstruction=.true., check_remapping=.true.) call extract_member_remapping_CS(CS%remap_CS, degree=CS%deg) CS%H_subroundoff = 1.0E-20 CS%debug=.false. CS%limiter=.false. CS%limiter_remap=.false. CS%hbd_nk = 2 + (2*2) + call initialize_remapping( CS%remap_CS, 'PLM', boundary_extrapolation=.true., & + om4_remap_via_sub_cells=.true., & ! ### see fail below when using fixed remapping alg. + check_reconstruction=.true., check_remapping=.true., & + h_neglect=CS%H_subroundoff, h_neglect_edge=CS%H_subroundoff) allocate(CS%hbd_grd_u(1,1,CS%hbd_nk), source=0.0) allocate(CS%hbd_u_kmax(1,1), source=0) near_boundary_unit_tests = .false.