From 3d1b6cd80561a4de2358f3516cf5b2ab7c77e350 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Tue, 27 Feb 2024 04:35:02 -0500 Subject: [PATCH] +Eliminated h_neglect argument to remapping_core_h Eliminated the previously optional h_neglect and h_neglect_edge arguments to remapping_core_h and remapping_core_w, and added optional h_neglect and h_neglect_edge arguments to initialize_remapping for storage in the remapping control structures. The answer_date, h_neglect and h_neglect_edge arguments to ALE_remap_scalar were also eliminated, as they are no longer used. The new routine open_boundary_setup_vert was added to manage these changes. A new verticalGrid_type argument was added to wave_speed_init. The h_neglect argument to 29 routines was made non-optional, because there is no way to provide a consistent hard-coded default for a dimensional parameter. The (mostly low-level) routines where this change to a non-optional h_neglect was made include build_reconstructions_1d, P1M_interpolation, P3M_interpolation, P3M_boundary_extrapolation, PLM_reconstruction, PLM_boundary_extrapolation, PPM_reconstruction, PPM_limiter_standard, PPM_boundary_extrapolation, PQM_reconstruction, PQM_limiter, PQM_boundary_extrapolation_v1, build_hycom1_column, build_rho_column, bound_edge_values, edge_values_explicit_h4, edge_values_explicit_h4cw, edge_values_implicit_h4, edge_slopes_implicit_h3, edge_slopes_implicit_h5, edge_values_implicit_h6, regridding_set_ppolys, build_and_interpolate_grid, remapByProjection, remapByDeltaZ, and integrateReconOnInterval. In two modules (MOM_open_boundary and MOM_wave_speed), two separate copies of remapping_CS variables were needed to accommodate different negligible thicknesses or units in different remapping calls. In those cases that also have an optional h_neglect_edge argument the default is now set to h_neglect. Improperly hard-coded dimensional default values for h_neglect or h_neglect_edge were eliminated in 16 places as a result of this change, but they were added back in 2 unit testing routines (one of these is archaic and unused) that do not use exercise dimensional rescaling tests. Because the previously optional arguments were already present in all calls from the dev/gfdl or main branches of the MOM6 code, and because the places where arguments have been removed they are balanced by equivalent arguments added to initialize_remapping, all answers are bitwise identical in the regression tests, but this change in interfaces could impact other code that is not in the main branch of MOM6. --- .../timing_tests/time_MOM_remapping.F90 | 7 +- src/ALE/MOM_ALE.F90 | 105 +++-------- src/ALE/MOM_remapping.F90 | 114 ++++++------ src/ALE/P1M_functions.F90 | 2 +- src/ALE/P3M_functions.F90 | 42 ++--- src/ALE/PLM_functions.F90 | 18 +- src/ALE/PPM_functions.F90 | 20 +- src/ALE/PQM_functions.F90 | 56 +++--- src/ALE/coord_hycom.F90 | 8 +- src/ALE/coord_rho.F90 | 8 +- src/ALE/regrid_edge_values.F90 | 112 ++++------- src/ALE/regrid_interp.F90 | 28 +-- src/ALE/remapping_attic.F90 | 38 ++-- src/core/MOM.F90 | 4 + src/core/MOM_open_boundary.F90 | 174 ++++++++++-------- src/diagnostics/MOM_diagnostics.F90 | 2 +- src/diagnostics/MOM_wave_speed.F90 | 41 +++-- src/framework/MOM_diag_remap.F90 | 24 +-- .../MOM_state_initialization.F90 | 61 +++--- .../MOM_tracer_initialization_from_Z.F90 | 27 +-- src/ocean_data_assim/MOM_oda_driver.F90 | 35 ++-- src/ocean_data_assim/MOM_oda_incupd.F90 | 47 ++--- .../lateral/MOM_internal_tides.F90 | 2 +- .../lateral/MOM_lateral_mixing_coeffs.F90 | 2 +- .../vertical/MOM_ALE_sponge.F90 | 37 ++-- .../vertical/MOM_tidal_mixing.F90 | 15 +- src/tracer/MOM_hor_bnd_diffusion.F90 | 21 +-- 27 files changed, 474 insertions(+), 576 deletions(-) diff --git a/config_src/drivers/timing_tests/time_MOM_remapping.F90 b/config_src/drivers/timing_tests/time_MOM_remapping.F90 index 5f4c0258ca..e4bea9d94f 100644 --- a/config_src/drivers/timing_tests/time_MOM_remapping.F90 +++ b/config_src/drivers/timing_tests/time_MOM_remapping.F90 @@ -17,8 +17,9 @@ program time_MOM_remapping real, dimension(nschemes) :: tmin ! Shortest time for a call [s] real, dimension(nschemes) :: tmax ! Longest time for a call [s] real, dimension(:,:), allocatable :: u0, u1 ! Source/target values [arbitrary but same units as each other] -real, dimension(:,:), allocatable :: h0, h1 ! Source target thicknesses [0..1] +real, dimension(:,:), allocatable :: h0, h1 ! Source target thicknesses [0..1] [nondim] real :: start, finish ! Times [s] +real :: h_neglect ! A negligible thickness [nondim] real :: h0sum, h1sum ! Totals of h0 and h1 [nondim] integer :: ij, k, isamp, iter, ischeme ! Indices and counters integer :: seed_size ! Number of integers used by seed @@ -50,6 +51,7 @@ program time_MOM_remapping h0(:,ij) = h0(:,ij) / h0sum h1(:,ij) = h1(:,ij) / h1sum enddo +h_neglect = 1.0-30 ! Loop over many samples of timing loop to collect statistics tmean(:) = 0. @@ -59,7 +61,8 @@ 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() diff --git a/src/ALE/MOM_ALE.F90 b/src/ALE/MOM_ALE.F90 index 0ae5fb1e92..bc3099d68d 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 "//& @@ -345,7 +356,7 @@ subroutine ALE_set_OM4_remap_algorithm( CS, om4_remap_via_sub_cells ) type(ALE_CS), pointer :: CS !< Module control structure logical, intent(in) :: om4_remap_via_sub_cells !< If true, use OM4 remapping algorithm - call remapping_set_param(CS%remapCS, om4_remap_via_sub_cells =om4_remap_via_sub_cells ) + call remapping_set_param(CS%remapCS, om4_remap_via_sub_cells=om4_remap_via_sub_cells ) end subroutine ALE_set_OM4_remap_algorithm @@ -591,8 +602,8 @@ subroutine ALE_offline_inputs(CS, G, GV, US, h, tv, Reg, uhtr, vhtr, Kd, debug, endif enddo ; enddo - call ALE_remap_scalar(CS%remapCS, G, GV, nk, h, tv%T, h_new, tv%T, answer_date=CS%answer_date) - call ALE_remap_scalar(CS%remapCS, G, GV, nk, h, tv%S, h_new, tv%S, answer_date=CS%answer_date) + call ALE_remap_scalar(CS%remapCS, G, GV, nk, h, tv%T, h_new, tv%T) + call ALE_remap_scalar(CS%remapCS, G, GV, nk, h, tv%S, h_new, tv%S) if (debug) call MOM_tracer_chkinv("After ALE_offline_inputs", G, GV, h_new, Reg%Tr, Reg%ntr) @@ -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) @@ -1300,8 +1278,7 @@ end subroutine mask_near_bottom_vel !> Remaps a single scalar between grids described by thicknesses h_src and h_dst. !! 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) +subroutine ALE_remap_scalar(CS, G, GV, nk_src, h_src, s_src, h_dst, s_dst, all_cells, old_remap) 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 @@ -1319,44 +1296,16 @@ subroutine ALE_remap_scalar(CS, G, GV, nk_src, h_src, s_src, h_dst, s_dst, all_c !! layers otherwise (default). logical, optional, intent(in) :: old_remap !< If true, use the old "remapping_core_w" !! method, otherwise use "remapping_core_h". - logical, optional, intent(in) :: answers_2018 !< If true, use the order of arithmetic - !! and expressions that recover the answers for - !! remapping from the end of 2018. Otherwise, - !! 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 + ! 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 + logical :: ignore_vanished_layers, use_remapping_core_w ignore_vanished_layers = .false. if (present(all_cells)) ignore_vanished_layers = .not. all_cells use_remapping_core_w = .false. if (present(old_remap)) use_remapping_core_w = old_remap n_points = nk_src - 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 @@ -1371,10 +1320,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 4495e4a170..3c2c0af6df 100644 --- a/src/ALE/MOM_remapping.F90 +++ b/src/ALE/MOM_remapping.F90 @@ -40,6 +40,13 @@ 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 the h0 argument to remapping_core_h [H] + real :: h_neglect + !> A negligibly small width for the purpose of edge value calculations in the same units + !! as the h0 argument to remapping_core_h [H] + real :: h_neglect_edge end type !> Class to assist in unit tests @@ -114,7 +121,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 +132,12 @@ 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 the h0 argument + !! to remapping_core_h [H] + real, optional, intent(in) :: h_neglect_edge !< A negligibly small width for the purpose of edge + !! value calculations in the same units as as the h0 + !! argument to remapping_core_h [H] if (present(remapping_scheme)) then call setReconstructionType( remapping_scheme, CS ) @@ -153,6 +167,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 @@ -181,7 +201,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 -subroutine remapping_core_h(CS, n0, h0, u0, n1, h1, u1, h_neglect, h_neglect_edge, net_err, PCM_cell) +subroutine remapping_core_h(CS, n0, h0, u0, n1, h1, u1, net_err, 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] @@ -189,12 +209,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, 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] real, optional, intent(out) :: net_err !< Error in total column [A H] logical, dimension(n0), optional, intent(in) :: PCM_cell !< If present, use PCM remapping for !! cells in the source grid where this is true. @@ -217,14 +231,10 @@ 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, hNeglect_edge ! Negligibly small cell widths in the same units as h0 [H] integer :: iMethod ! An integer indicating the integration method used - hNeglect = 1.0e-30 ; if (present(h_neglect)) hNeglect = h_neglect - hNeglect_edge = 1.0e-10 ; 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, & - hNeglect, hNeglect_edge, PCM_cell ) + CS%h_neglect, CS%h_neglect_edge, PCM_cell ) if (CS%om4_remap_via_sub_cells) then @@ -284,7 +294,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] @@ -292,12 +302,7 @@ 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, 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]. + ! Local variables real, dimension(n0+n1+1) :: h_sub ! Width of each each sub-cell [H] real, dimension(n0+n1+1) :: uh_sub ! Integral of u*h over each sub-cell [A H] @@ -317,15 +322,11 @@ subroutine remapping_core_w( CS, n0, h0, u0, n1, dx, u1, h_neglect, h_neglect_ed 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, hNeglect_edge ! Negligibly small thicknesses [H] integer :: iMethod ! An integer indicating the integration method used integer :: k - hNeglect = 1.0e-30 ; if (present(h_neglect)) hNeglect = h_neglect - hNeglect_edge = 1.0e-10 ; 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,& - hNeglect, 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) @@ -377,19 +378,23 @@ subroutine build_reconstructions_1d( CS, n0, h0, u0, ppoly_r_coefs, & real, dimension(n0,2), intent(out) :: ppoly_r_E !< Edge value of polynomial [A] real, dimension(n0,2), intent(out) :: ppoly_r_S !< Edge slope of polynomial [A H-1] integer, intent(out) :: iMethod !< Integration method - real, optional, intent(in) :: h_neglect !< A negligibly small width for the + 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] + 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, optional, intent(in) :: PCM_cell(n0) !< If present, use PCM remapping for !! cells from the source grid where this is true. ! Local variables + real :: h_neg_edge ! A negligibly small width for the purpose of edge value + ! calculations in the same units as h0 [H] integer :: local_remapping_scheme integer :: k, n + h_neg_edge = h_neglect ; if (present(h_neglect_edge)) h_neg_edge = h_neglect_edge + ! Reset polynomial ppoly_r_E(:,:) = 0.0 ppoly_r_S(:,:) = 0.0 @@ -426,7 +431,7 @@ subroutine build_reconstructions_1d( CS, n0, h0, u0, ppoly_r_coefs, & iMethod = INTEGRATION_PLM case ( REMAPPING_PPM_CW ) ! identical to REMAPPING_PPM_HYBGEN - call edge_values_explicit_h4cw( n0, h0, u0, ppoly_r_E, h_neglect_edge ) + call edge_values_explicit_h4cw( n0, h0, u0, ppoly_r_E, h_neg_edge ) call PPM_monotonicity( n0, u0, ppoly_r_E ) call PPM_reconstruction( n0, h0, u0, ppoly_r_E, ppoly_r_coefs, h_neglect, answer_date=CS%answer_date ) if ( CS%boundary_extrapolation ) then @@ -434,14 +439,14 @@ subroutine build_reconstructions_1d( CS, n0, h0, u0, ppoly_r_coefs, & endif iMethod = INTEGRATION_PPM case ( REMAPPING_PPM_H4 ) - call edge_values_explicit_h4( n0, h0, u0, ppoly_r_E, h_neglect_edge, answer_date=CS%answer_date ) + call edge_values_explicit_h4( n0, h0, u0, ppoly_r_E, h_neg_edge, answer_date=CS%answer_date ) call PPM_reconstruction( n0, h0, u0, ppoly_r_E, ppoly_r_coefs, h_neglect, answer_date=CS%answer_date ) if ( CS%boundary_extrapolation ) then call PPM_boundary_extrapolation( n0, h0, u0, ppoly_r_E, ppoly_r_coefs, h_neglect ) endif iMethod = INTEGRATION_PPM case ( REMAPPING_PPM_IH4 ) - call edge_values_implicit_h4( n0, h0, u0, ppoly_r_E, h_neglect_edge, answer_date=CS%answer_date ) + call edge_values_implicit_h4( n0, h0, u0, ppoly_r_E, h_neg_edge, answer_date=CS%answer_date ) call PPM_reconstruction( n0, h0, u0, ppoly_r_E, ppoly_r_coefs, h_neglect, answer_date=CS%answer_date ) if ( CS%boundary_extrapolation ) then call PPM_boundary_extrapolation( n0, h0, u0, ppoly_r_E, ppoly_r_coefs, h_neglect ) @@ -460,7 +465,7 @@ subroutine build_reconstructions_1d( CS, n0, h0, u0, ppoly_r_coefs, & call PPM_boundary_extrapolation( n0, h0, u0, ppoly_r_E, ppoly_r_coefs, h_neglect ) iMethod = INTEGRATION_PPM case ( REMAPPING_PQM_IH4IH3 ) - call edge_values_implicit_h4( n0, h0, u0, ppoly_r_E, h_neglect_edge, answer_date=CS%answer_date ) + call edge_values_implicit_h4( n0, h0, u0, ppoly_r_E, h_neg_edge, answer_date=CS%answer_date ) call edge_slopes_implicit_h3( n0, h0, u0, ppoly_r_S, h_neglect, answer_date=CS%answer_date ) call PQM_reconstruction( n0, h0, u0, ppoly_r_E, ppoly_r_S, ppoly_r_coefs, h_neglect, & answer_date=CS%answer_date ) @@ -470,7 +475,7 @@ subroutine build_reconstructions_1d( CS, n0, h0, u0, ppoly_r_coefs, & endif iMethod = INTEGRATION_PQM case ( REMAPPING_PQM_IH6IH5 ) - call edge_values_implicit_h6( n0, h0, u0, ppoly_r_E, h_neglect_edge, answer_date=CS%answer_date ) + call edge_values_implicit_h6( n0, h0, u0, ppoly_r_E, h_neg_edge, answer_date=CS%answer_date ) call edge_slopes_implicit_h5( n0, h0, u0, ppoly_r_S, h_neglect, answer_date=CS%answer_date ) call PQM_reconstruction( n0, h0, u0, ppoly_r_E, ppoly_r_S, ppoly_r_coefs, h_neglect, & answer_date=CS%answer_date ) @@ -1483,7 +1488,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 @@ -1494,12 +1500,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 @@ -1580,21 +1591,17 @@ logical function remapping_unit_tests(verbose) real :: u02_err ! Error in remaping [A] integer, allocatable, dimension(:) :: isrc_start, isrc_end, isrc_max, itgt_start, itgt_end, isub_src ! Indices integer :: answer_date ! The vintage of the expressions to test - real, parameter :: hNeglect_dflt = 1.0e-30 ! A thickness [H ~> m or kg m-2] that can be - ! added to thicknesses in a denominator without - ! changing the numerical result, except where - ! a division by zero would otherwise occur. real :: err ! Errors in the remapped thicknesses [H] or values [A] real :: h_neglect, h_neglect_edge ! Tiny thicknesses used in remapping [H] type(testing) :: test ! Unit testing convenience functions - integer :: om4 + integer :: i, om4 character(len=4) :: om4_tag call test%set( verbose=verbose ) ! Sets the verbosity flag in test answer_date = 20190101 ! 20181231 - h_neglect = hNeglect_dflt - h_neglect_edge = hNeglect_dflt ; if (answer_date < 20190101) h_neglect_edge = 1.0e-10 + h_neglect = 1.0e-30 + h_neglect_edge = h_neglect ; if (answer_date < 20190101) h_neglect_edge = 1.0e-10 if (verbose) write(test%stdout,*) ' ===== MOM_remapping: remapping_unit_tests =================' @@ -1603,7 +1610,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 @@ -1623,7 +1631,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)) @@ -2067,7 +2075,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 @@ -2079,27 +2087,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 @@ -2107,7 +2115,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/P1M_functions.F90 b/src/ALE/P1M_functions.F90 index 7889966135..d2051cc702 100644 --- a/src/ALE/P1M_functions.F90 +++ b/src/ALE/P1M_functions.F90 @@ -31,7 +31,7 @@ subroutine P1M_interpolation( N, h, u, edge_values, ppoly_coef, h_neglect, answe real, dimension(:,:), intent(inout) :: edge_values !< Potentially modified edge values [A] real, dimension(:,:), intent(inout) :: ppoly_coef !< Potentially modified !! piecewise polynomial coefficients, mainly [A] - real, optional, intent(in) :: h_neglect !< A negligibly small width [H] + real, intent(in) :: h_neglect !< A negligibly small width [H] integer, optional, intent(in) :: answer_date !< The vintage of the expressions to use ! Local variables diff --git a/src/ALE/P3M_functions.F90 b/src/ALE/P3M_functions.F90 index 6039b197fb..e9c234db32 100644 --- a/src/ALE/P3M_functions.F90 +++ b/src/ALE/P3M_functions.F90 @@ -10,9 +10,6 @@ module P3M_functions public P3M_interpolation public P3M_boundary_extrapolation -real, parameter :: hNeglect_dflt = 1.E-30 !< Default value of a negligible cell thickness -real, parameter :: hNeglect_edge_dflt = 1.E-10 !< Default value of a negligible edge thickness - contains !> Set up a piecewise cubic interpolation from cell averages and estimated @@ -32,7 +29,7 @@ subroutine P3M_interpolation( N, h, u, edge_values, ppoly_S, ppoly_coef, h_negle real, dimension(:,:), intent(inout) :: edge_values !< Edge value of polynomial [A] real, dimension(:,:), intent(inout) :: ppoly_S !< Edge slope of polynomial [A H-1]. real, dimension(:,:), intent(inout) :: ppoly_coef !< Coefficients of polynomial [A] - real, optional, intent(in) :: h_neglect !< A negligibly small width for the + real, intent(in) :: h_neglect !< A negligibly small width for the !! purpose of cell reconstructions [H] integer, optional, intent(in) :: answer_date !< The vintage of the expressions to use @@ -66,7 +63,7 @@ subroutine P3M_limiter( N, h, u, edge_values, ppoly_S, ppoly_coef, h_neglect, an real, dimension(:,:), intent(inout) :: edge_values !< Edge value of polynomial [A] real, dimension(:,:), intent(inout) :: ppoly_S !< Edge slope of polynomial [A H-1] real, dimension(:,:), intent(inout) :: ppoly_coef !< Coefficients of polynomial [A] - real, optional, intent(in) :: h_neglect !< A negligibly small width for + real, intent(in) :: h_neglect !< A negligibly small width for !! the purpose of cell reconstructions [H] integer, optional, intent(in) :: answer_date !< The vintage of the expressions to use @@ -79,15 +76,9 @@ subroutine P3M_limiter( N, h, u, edge_values, ppoly_S, ppoly_coef, h_neglect, an real :: h_l, h_c, h_r ! left, center and right cell widths [H] real :: sigma_l, sigma_c, sigma_r ! left, center and right van Leer slopes [A H-1] real :: slope ! retained PLM slope [A H-1] - real :: eps - real :: hNeglect ! A negligibly small thickness [H] - - hNeglect = hNeglect_dflt ; if (present(h_neglect)) hNeglect = h_neglect - - eps = 1e-10 ! 1. Bound edge values (boundary cells are assumed to be local extrema) - call bound_edge_values( N, h, u, edge_values, hNeglect, answer_date=answer_date ) + call bound_edge_values( N, h, u, edge_values, h_neglect, answer_date=answer_date ) ! 2. Systematically average discontinuous edge values call average_discontinuous_edge_values( N, edge_values ) @@ -127,9 +118,9 @@ subroutine P3M_limiter( N, h, u, edge_values, ppoly_S, ppoly_coef, h_neglect, an endif ! Compute limited slope - sigma_l = 2.0 * ( u_c - u_l ) / ( h_c + hNeglect ) - sigma_c = 2.0 * ( u_r - u_l ) / ( h_l + 2.0*h_c + h_r + hNeglect ) - sigma_r = 2.0 * ( u_r - u_c ) / ( h_c + hNeglect ) + sigma_l = 2.0 * ( u_c - u_l ) / ( h_c + h_neglect ) + sigma_c = 2.0 * ( u_r - u_l ) / ( h_l + 2.0*h_c + h_r + h_neglect ) + sigma_r = 2.0 * ( u_r - u_c ) / ( h_c + h_neglect ) if ( (sigma_l * sigma_r) > 0.0 ) then slope = sign( min(abs(sigma_l),abs(sigma_c),abs(sigma_r)), sigma_c ) @@ -197,10 +188,10 @@ subroutine P3M_boundary_extrapolation( N, h, u, edge_values, ppoly_S, ppoly_coef real, dimension(:,:), intent(inout) :: edge_values !< Edge value of polynomial [A] real, dimension(:,:), intent(inout) :: ppoly_S !< Edge slope of polynomial [A H-1] real, dimension(:,:), intent(inout) :: ppoly_coef !< Coefficients of polynomial [A] - real, optional, intent(in) :: h_neglect !< A negligibly small width for the + real, intent(in) :: h_neglect !< A negligibly small width for the !! purpose of cell reconstructions [H] - real, optional, intent(in) :: h_neglect_edge !< A negligibly small width - !! for the purpose of finding edge values [H] + real, optional, intent(in) :: h_neglect_edge !< A negligibly small width for the purpose + !! of finding edge values [H]. The default is h_neglect. ! Local variables integer :: i0, i1 logical :: monotonic ! boolean indicating whether the cubic is monotonic @@ -210,10 +201,9 @@ subroutine P3M_boundary_extrapolation( N, h, u, edge_values, ppoly_S, ppoly_coef real :: u0_l, u0_r ! Left and right edge values [A] real :: u1_l, u1_r ! Left and right edge slopes [A H-1] real :: slope ! The cell center slope [A H-1] - real :: hNeglect, hNeglect_edge ! Negligibly small thickness [H] + real :: hNeglect_edge ! Negligibly small thickness [H] - hNeglect = hNeglect_dflt ; if (present(h_neglect)) hNeglect = h_neglect - hNeglect_edge = hNeglect_edge_dflt ; if (present(h_neglect_edge)) hNeglect_edge = h_neglect_edge + hNeglect_edge = h_neglect ; if (present(h_neglect_edge)) hNeglect_edge = h_neglect_edge ! ----- Left boundary ----- i0 = 1 @@ -229,7 +219,7 @@ subroutine P3M_boundary_extrapolation( N, h, u, edge_values, ppoly_S, ppoly_coef u1_r = b / h1 ! derivative evaluated at xi = 0.0, expressed w.r.t. x ! Limit the right slope by the PLM limited slope - slope = 2.0 * ( u1 - u0 ) / ( h0 + hNeglect ) + slope = 2.0 * ( u1 - u0 ) / ( h0 + h_neglect ) if ( abs(u1_r) > abs(slope) ) then u1_r = slope endif @@ -242,7 +232,7 @@ subroutine P3M_boundary_extrapolation( N, h, u, edge_values, ppoly_S, ppoly_coef ! edge value and slope by computing the parabola as determined by ! the right edge value and slope and the boundary cell average u0_l = 3.0 * u0 + 0.5 * h0*u1_r - 2.0 * u0_r - u1_l = ( - 6.0 * u0 - 2.0 * h0*u1_r + 6.0 * u0_r) / ( h0 + hNeglect ) + u1_l = ( - 6.0 * u0 - 2.0 * h0*u1_r + 6.0 * u0_r) / ( h0 + h_neglect ) ! Check whether the edge values are monotonic. For example, if the left edge ! value is larger than the right edge value while the slope is positive, the @@ -286,10 +276,10 @@ subroutine P3M_boundary_extrapolation( N, h, u, edge_values, ppoly_S, ppoly_coef b = ppoly_coef(i0,2) c = ppoly_coef(i0,3) d = ppoly_coef(i0,4) - u1_l = (b + 2*c + 3*d) / ( h0 + hNeglect ) ! derivative evaluated at xi = 1.0 + u1_l = (b + 2*c + 3*d) / ( h0 + h_neglect ) ! derivative evaluated at xi = 1.0 ! Limit the left slope by the PLM limited slope - slope = 2.0 * ( u1 - u0 ) / ( h1 + hNeglect ) + slope = 2.0 * ( u1 - u0 ) / ( h1 + h_neglect ) if ( abs(u1_l) > abs(slope) ) then u1_l = slope endif @@ -302,7 +292,7 @@ subroutine P3M_boundary_extrapolation( N, h, u, edge_values, ppoly_S, ppoly_coef ! edge value and slope by computing the parabola as determined by ! the left edge value and slope and the boundary cell average u0_r = 3.0 * u1 - 0.5 * h1*u1_l - 2.0 * u0_l - u1_r = ( 6.0 * u1 - 2.0 * h1*u1_l - 6.0 * u0_l) / ( h1 + hNeglect ) + u1_r = ( 6.0 * u1 - 2.0 * h1*u1_l - 6.0 * u0_l) / ( h1 + h_neglect ) ! Check whether the edge values are monotonic. For example, if the right edge ! value is smaller than the left edge value while the slope is positive, the diff --git a/src/ALE/PLM_functions.F90 b/src/ALE/PLM_functions.F90 index c0c4516fe2..6d6afd3885 100644 --- a/src/ALE/PLM_functions.F90 +++ b/src/ALE/PLM_functions.F90 @@ -12,8 +12,6 @@ module PLM_functions public PLM_slope_wa public PLM_slope_cw -real, parameter :: hNeglect_dflt = 1.E-30 !< Default negligible cell thickness - contains !> Returns a limited PLM slope following White and Adcroft, 2008, in the same arbitrary @@ -195,7 +193,7 @@ subroutine PLM_reconstruction( N, h, u, edge_values, ppoly_coef, h_neglect ) !! with the same units as u [A]. real, dimension(:,:), intent(inout) :: ppoly_coef !< coefficients of piecewise polynomials, mainly !! with the same units as u [A]. - real, optional, intent(in) :: h_neglect !< A negligibly small width for + real, intent(in) :: h_neglect !< A negligibly small width for !! the purpose of cell reconstructions !! in the same units as h [H] @@ -208,15 +206,12 @@ subroutine PLM_reconstruction( N, h, u, edge_values, ppoly_coef, h_neglect ) real :: almost_one ! A value that is slightly smaller than 1 [nondim] real, dimension(N) :: slp ! The first guess at the normalized tracer slopes [A] real, dimension(N) :: mslp ! The monotonized normalized tracer slopes [A] - real :: hNeglect ! A negligibly small width used in cell reconstructions [H] - - hNeglect = hNeglect_dflt ; if (present(h_neglect)) hNeglect = h_neglect almost_one = 1. - epsilon(slope) ! Loop on interior cells do k = 2,N-1 - slp(k) = PLM_slope_wa(h(k-1), h(k), h(k+1), hNeglect, u(k-1), u(k), u(k+1)) + slp(k) = PLM_slope_wa(h(k-1), h(k), h(k+1), h_neglect, u(k-1), u(k), u(k+1)) enddo ! end loop on interior cells ! Boundary cells use PCM. Extrapolation is handled after monotonization. @@ -277,17 +272,14 @@ subroutine PLM_boundary_extrapolation( N, h, u, edge_values, ppoly_coef, h_negle !! with the same units as u [A]. real, dimension(:,:), intent(inout) :: ppoly_coef !< coefficients of piecewise polynomials, mainly !! with the same units as u [A]. - real, optional, intent(in) :: h_neglect !< A negligibly small width for + real, intent(in) :: h_neglect !< A negligibly small width for !! the purpose of cell reconstructions !! in the same units as h [H] ! Local variables real :: slope ! retained PLM slope for a normalized cell width [A] - real :: hNeglect ! A negligibly small width used in cell reconstructions [H] - - hNeglect = hNeglect_dflt ; if (present(h_neglect)) hNeglect = h_neglect ! Extrapolate from 2 to 1 to estimate slope - slope = - PLM_extrapolate_slope( h(2), h(1), hNeglect, u(2), u(1) ) + slope = - PLM_extrapolate_slope( h(2), h(1), h_neglect, u(2), u(1) ) edge_values(1,1) = u(1) - 0.5 * slope edge_values(1,2) = u(1) + 0.5 * slope @@ -296,7 +288,7 @@ subroutine PLM_boundary_extrapolation( N, h, u, edge_values, ppoly_coef, h_negle ppoly_coef(1,2) = edge_values(1,2) - edge_values(1,1) ! Extrapolate from N-1 to N to estimate slope - slope = PLM_extrapolate_slope( h(N-1), h(N), hNeglect, u(N-1), u(N) ) + slope = PLM_extrapolate_slope( h(N-1), h(N), h_neglect, u(N-1), u(N) ) edge_values(N,1) = u(N) - 0.5 * slope edge_values(N,2) = u(N) + 0.5 * slope diff --git a/src/ALE/PPM_functions.F90 b/src/ALE/PPM_functions.F90 index ef6841f635..c11ec6e741 100644 --- a/src/ALE/PPM_functions.F90 +++ b/src/ALE/PPM_functions.F90 @@ -15,13 +15,6 @@ module PPM_functions public PPM_reconstruction, PPM_boundary_extrapolation, PPM_monotonicity -!> A tiny width that is so small that adding it to cell widths does not -!! change the value due to a computational representation. It is used -!! to avoid division by zero. -!! @note This is a dimensional parameter and should really include a unit -!! conversion. -real, parameter :: hNeglect_dflt = 1.E-30 - contains !> Builds quadratic polynomials coefficients from cell mean and edge values. @@ -31,7 +24,7 @@ subroutine PPM_reconstruction( N, h, u, edge_values, ppoly_coef, h_neglect, answ real, dimension(N), intent(in) :: u !< Cell averages in arbitrary coordinates [A] real, dimension(N,2), intent(inout) :: edge_values !< Edge values [A] real, dimension(N,3), intent(inout) :: ppoly_coef !< Polynomial coefficients, mainly [A] - real, optional, intent(in) :: h_neglect !< A negligibly small width [H] + real, intent(in) :: h_neglect !< A negligibly small width [H] integer, optional, intent(in) :: answer_date !< The vintage of the expressions to use ! Local variables @@ -64,7 +57,7 @@ subroutine PPM_limiter_standard( N, h, u, edge_values, h_neglect, answer_date ) real, dimension(:), intent(in) :: h !< cell widths (size N) [H] real, dimension(:), intent(in) :: u !< cell average properties (size N) [A] real, dimension(:,:), intent(inout) :: edge_values !< Potentially modified edge values [A] - real, optional, intent(in) :: h_neglect !< A negligibly small width [H] + real, intent(in) :: h_neglect !< A negligibly small width [H] integer, optional, intent(in) :: answer_date !< The vintage of the expressions to use ! Local variables @@ -190,7 +183,7 @@ subroutine PPM_boundary_extrapolation( N, h, u, edge_values, ppoly_coef, h_negle real, dimension(:), intent(in) :: u !< cell averages (size N) [A] real, dimension(:,:), intent(inout) :: edge_values !< edge values of piecewise polynomials [A] real, dimension(:,:), intent(inout) :: ppoly_coef !< coefficients of piecewise polynomials, mainly [A] - real, optional, intent(in) :: h_neglect !< A negligibly small width for + real, intent(in) :: h_neglect !< A negligibly small width for !! the purpose of cell reconstructions [H] ! Local variables @@ -204,9 +197,6 @@ subroutine PPM_boundary_extrapolation( N, h, u, edge_values, ppoly_coef, h_negle ! the cell being worked on [A] real :: slope ! The normalized slope [A] real :: exp1, exp2 ! Temporary expressions [A2] - real :: hNeglect ! A negligibly small width used in cell reconstructions [H] - - hNeglect = hNeglect_dflt ; if (present(h_neglect)) hNeglect = h_neglect ! ----- Left boundary ----- i0 = 1 @@ -219,7 +209,7 @@ subroutine PPM_boundary_extrapolation( N, h, u, edge_values, ppoly_coef, h_negle ! Compute the left edge slope in neighboring cell and express it in ! the global coordinate system b = ppoly_coef(i1,2) - u1_r = b *((h0+hNeglect)/(h1+hNeglect)) ! derivative evaluated at xi = 0.0, + u1_r = b *((h0+h_neglect)/(h1+h_neglect)) ! derivative evaluated at xi = 0.0, ! expressed w.r.t. xi (local coord. system) ! Limit the right slope by the PLM limited slope @@ -273,7 +263,7 @@ subroutine PPM_boundary_extrapolation( N, h, u, edge_values, ppoly_coef, h_negle b = ppoly_coef(i0,2) c = ppoly_coef(i0,3) u1_l = (b + 2*c) ! derivative evaluated at xi = 1.0 - u1_l = u1_l * ((h1+hNeglect)/(h0+hNeglect)) + u1_l = u1_l * ((h1+h_neglect)/(h0+h_neglect)) ! Limit the left slope by the PLM limited slope slope = 2.0 * ( u1 - u0 ) diff --git a/src/ALE/PQM_functions.F90 b/src/ALE/PQM_functions.F90 index ef42fb9f01..418a4b47a2 100644 --- a/src/ALE/PQM_functions.F90 +++ b/src/ALE/PQM_functions.F90 @@ -9,8 +9,6 @@ module PQM_functions public PQM_reconstruction, PQM_boundary_extrapolation, PQM_boundary_extrapolation_v1 -real, parameter :: hNeglect_dflt = 1.E-30 !< Default negligible cell thickness - contains !> Reconstruction by quartic polynomials within each cell. @@ -24,7 +22,7 @@ subroutine PQM_reconstruction( N, h, u, edge_values, edge_slopes, ppoly_coef, h_ real, dimension(:,:), intent(inout) :: edge_values !< Edge value of polynomial [A] real, dimension(:,:), intent(inout) :: edge_slopes !< Edge slope of polynomial [A H-1] real, dimension(:,:), intent(inout) :: ppoly_coef !< Coefficients of polynomial, mainly [A] - real, optional, intent(in) :: h_neglect !< A negligibly small width for + real, intent(in) :: h_neglect !< A negligibly small width for !! the purpose of cell reconstructions [H] integer, optional, intent(in) :: answer_date !< The vintage of the expressions to use @@ -78,7 +76,7 @@ subroutine PQM_limiter( N, h, u, edge_values, edge_slopes, h_neglect, answer_dat real, dimension(:), intent(in) :: u !< cell average properties (size N) [A] real, dimension(:,:), intent(inout) :: edge_values !< Potentially modified edge values [A] real, dimension(:,:), intent(inout) :: edge_slopes !< Potentially modified edge slopes [A H-1] - real, optional, intent(in) :: h_neglect !< A negligibly small width for + real, intent(in) :: h_neglect !< A negligibly small width for !! the purpose of cell reconstructions [H] integer, optional, intent(in) :: answer_date !< The vintage of the expressions to use @@ -98,12 +96,9 @@ subroutine PQM_limiter( N, h, u, edge_values, edge_slopes, h_neglect, answer_dat real :: sqrt_rho ! The square root of rho [A] real :: gradient1, gradient2 ! Normalized gradients [A] real :: x1, x2 ! Fractional inflection point positions in a cell [nondim] - real :: hNeglect ! A negligibly small width for the purpose of cell reconstructions [H] - - hNeglect = hNeglect_dflt ; if (present(h_neglect)) hNeglect = h_neglect ! Bound edge values - call bound_edge_values( N, h, u, edge_values, hNeglect, answer_date=answer_date ) + call bound_edge_values( N, h, u, edge_values, h_neglect, answer_date=answer_date ) ! Make discontinuous edge values monotonic (thru averaging) call check_discontinuous_edge_values( N, u, edge_values ) @@ -132,9 +127,9 @@ subroutine PQM_limiter( N, h, u, edge_values, edge_slopes, h_neglect, answer_dat u_r = u(k+1) ! Compute limited slope - sigma_l = 2.0 * ( u_c - u_l ) / ( h_c + hNeglect ) - sigma_c = 2.0 * ( u_r - u_l ) / ( h_l + 2.0*h_c + h_r + hNeglect ) - sigma_r = 2.0 * ( u_r - u_c ) / ( h_c + hNeglect ) + sigma_l = 2.0 * ( u_c - u_l ) / ( h_c + h_neglect ) + sigma_c = 2.0 * ( u_r - u_l ) / ( h_l + 2.0*h_c + h_r + h_neglect ) + sigma_r = 2.0 * ( u_r - u_c ) / ( h_c + h_neglect ) if ( (sigma_l * sigma_r) > 0.0 ) then slope = sign( min(abs(sigma_l),abs(sigma_c),abs(sigma_r)), sigma_c ) @@ -272,8 +267,8 @@ subroutine PQM_limiter( N, h, u, edge_values, edge_slopes, h_neglect, answer_dat ! We modify the edge slopes so that both inflexion points ! collapse onto the left edge - u1_l = ( 10.0 * u_c - 2.0 * u0_r - 8.0 * u0_l ) / (3.0*h_c + hNeglect ) - u1_r = ( -10.0 * u_c + 6.0 * u0_r + 4.0 * u0_l ) / ( h_c + hNeglect ) + u1_l = ( 10.0 * u_c - 2.0 * u0_r - 8.0 * u0_l ) / (3.0*h_c + h_neglect ) + u1_r = ( -10.0 * u_c + 6.0 * u0_r + 4.0 * u0_l ) / ( h_c + h_neglect ) ! One of the modified slopes might be inconsistent. When that happens, ! the inconsistent slope is set equal to zero and the opposite edge value @@ -283,13 +278,13 @@ subroutine PQM_limiter( N, h, u, edge_values, edge_slopes, h_neglect, answer_dat u1_l = 0.0 u0_r = 5.0 * u_c - 4.0 * u0_l - u1_r = 20.0 * (u_c - u0_l) / ( h_c + hNeglect ) + u1_r = 20.0 * (u_c - u0_l) / ( h_c + h_neglect ) elseif ( u1_r * slope < 0.0 ) then u1_r = 0.0 u0_l = (5.0*u_c - 3.0*u0_r) / 2.0 - u1_l = 10.0 * (-u_c + u0_r) / (3.0 * h_c + hNeglect) + u1_l = 10.0 * (-u_c + u0_r) / (3.0 * h_c + h_neglect) endif @@ -297,8 +292,8 @@ subroutine PQM_limiter( N, h, u, edge_values, edge_slopes, h_neglect, answer_dat ! We modify the edge slopes so that both inflexion points ! collapse onto the right edge - u1_r = ( -10.0 * u_c + 8.0 * u0_r + 2.0 * u0_l ) / (3.0 * h_c + hNeglect) - u1_l = ( 10.0 * u_c - 4.0 * u0_r - 6.0 * u0_l ) / (h_c + hNeglect) + u1_r = ( -10.0 * u_c + 8.0 * u0_r + 2.0 * u0_l ) / (3.0 * h_c + h_neglect) + u1_l = ( 10.0 * u_c - 4.0 * u0_r - 6.0 * u0_l ) / (h_c + h_neglect) ! One of the modified slopes might be inconsistent. When that happens, ! the inconsistent slope is set equal to zero and the opposite edge value @@ -308,13 +303,13 @@ subroutine PQM_limiter( N, h, u, edge_values, edge_slopes, h_neglect, answer_dat u1_l = 0.0 u0_r = ( 5.0 * u_c - 3.0 * u0_l ) / 2.0 - u1_r = 10.0 * (u_c - u0_l) / (3.0 * h_c + hNeglect) + u1_r = 10.0 * (u_c - u0_l) / (3.0 * h_c + h_neglect) elseif ( u1_r * slope < 0.0 ) then u1_r = 0.0 u0_l = 5.0 * u_c - 4.0 * u0_r - u1_l = 20.0 * ( -u_c + u0_r ) / (h_c + hNeglect) + u1_l = 20.0 * ( -u_c + u0_r ) / (h_c + h_neglect) endif @@ -506,7 +501,7 @@ subroutine PQM_boundary_extrapolation_v1( N, h, u, edge_values, edge_slopes, ppo real, dimension(:,:), intent(inout) :: edge_values !< Edge value of polynomial [A] real, dimension(:,:), intent(inout) :: edge_slopes !< Edge slope of polynomial [A H-1] real, dimension(:,:), intent(inout) :: ppoly_coef !< Coefficients of polynomial, mainly [A] - real, optional, intent(in) :: h_neglect !< A negligibly small width for + real, intent(in) :: h_neglect !< A negligibly small width for !! the purpose of cell reconstructions [H] ! Local variables integer :: i0, i1 @@ -526,9 +521,6 @@ subroutine PQM_boundary_extrapolation_v1( N, h, u, edge_values, edge_slopes, ppo real :: sqrt_rho ! The square root of rho [A] real :: gradient1, gradient2 ! Normalized gradients [A] real :: x1, x2 ! Fractional inflection point positions in a cell [nondim] - real :: hNeglect ! A negligibly small width for the purpose of cell reconstructions [H] - - hNeglect = hNeglect_dflt ; if (present(h_neglect)) hNeglect = h_neglect ! ----- Left boundary (TOP) ----- i0 = 1 @@ -541,7 +533,7 @@ subroutine PQM_boundary_extrapolation_v1( N, h, u, edge_values, edge_slopes, ppo ! Compute real slope and express it w.r.t. local coordinate system ! within boundary cell - slope = 2.0 * ( u1 - u0 ) / ( ( h0 + h1 ) + hNeglect ) + slope = 2.0 * ( u1 - u0 ) / ( ( h0 + h1 ) + h_neglect ) slope = slope * h0 ! The right edge value and slope of the boundary cell are taken to be the @@ -550,12 +542,12 @@ subroutine PQM_boundary_extrapolation_v1( N, h, u, edge_values, edge_slopes, ppo b = ppoly_coef(i1,2) u0_r = a ! edge value - u1_r = b / (h1 + hNeglect) ! edge slope (w.r.t. global coord.) + u1_r = b / (h1 + h_neglect) ! edge slope (w.r.t. global coord.) ! Compute coefficient for rational function based on mean and right ! edge value and slope if (u1_r /= 0.) then ! HACK by AJA - beta = 2.0 * ( u0_r - um ) / ( (h0 + hNeglect)*u1_r) - 1.0 + beta = 2.0 * ( u0_r - um ) / ( (h0 + h_neglect)*u1_r) - 1.0 else beta = 0. endif ! HACK by AJA @@ -574,10 +566,10 @@ subroutine PQM_boundary_extrapolation_v1( N, h, u, edge_values, edge_slopes, ppo ! compute corresponding slope. if ( abs(um-u0_l) < abs(um-u_plm) ) then u1_l = 2.0 * ( br - ar*beta) - u1_l = u1_l / (h0 + hNeglect) + u1_l = u1_l / (h0 + h_neglect) else u0_l = u_plm - u1_l = slope / (h0 + hNeglect) + u1_l = slope / (h0 + h_neglect) endif ! Monotonize quartic @@ -635,8 +627,8 @@ subroutine PQM_boundary_extrapolation_v1( N, h, u, edge_values, edge_slopes, ppo ! We modify the edge slopes so that both inflexion points ! collapse onto the left edge - u1_l = ( 10.0 * um - 2.0 * u0_r - 8.0 * u0_l ) / (3.0*h0 + hNeglect) - u1_r = ( -10.0 * um + 6.0 * u0_r + 4.0 * u0_l ) / (h0 + hNeglect) + u1_l = ( 10.0 * um - 2.0 * u0_r - 8.0 * u0_l ) / (3.0*h0 + h_neglect) + u1_r = ( -10.0 * um + 6.0 * u0_r + 4.0 * u0_l ) / (h0 + h_neglect) ! One of the modified slopes might be inconsistent. When that happens, ! the inconsistent slope is set equal to zero and the opposite edge value @@ -646,13 +638,13 @@ subroutine PQM_boundary_extrapolation_v1( N, h, u, edge_values, edge_slopes, ppo u1_l = 0.0 u0_r = 5.0 * um - 4.0 * u0_l - u1_r = 20.0 * (um - u0_l) / ( h0 + hNeglect ) + u1_r = 20.0 * (um - u0_l) / ( h0 + h_neglect ) elseif ( u1_r * slope < 0.0 ) then u1_r = 0.0 u0_l = (5.0*um - 3.0*u0_r) / 2.0 - u1_l = 10.0 * (-um + u0_r) / (3.0 * h0 + hNeglect ) + u1_l = 10.0 * (-um + u0_r) / (3.0 * h0 + h_neglect ) endif diff --git a/src/ALE/coord_hycom.F90 b/src/ALE/coord_hycom.F90 index ddc569e45e..1e5474770a 100644 --- a/src/ALE/coord_hycom.F90 +++ b/src/ALE/coord_hycom.F90 @@ -119,7 +119,7 @@ subroutine build_hycom1_column(CS, remapCS, eqn_of_state, nz, depth, h, T, S, p_ real, optional, intent(in) :: zScale !< Scaling factor from the input coordinate thicknesses in [Z ~> m] !! to desired units for zInterface, perhaps GV%Z_to_H in which !! case this has units of [H Z-1 ~> nondim or kg m-3] - real, optional, intent(in) :: h_neglect !< A negligibly small width for the purpose of + real, intent(in) :: h_neglect !< A negligibly small width for the purpose of !! cell reconstruction [H ~> m or kg m-2] real, optional, intent(in) :: h_neglect_edge !< A negligibly small width for the purpose of !! edge value calculation [H ~> m or kg m-2] @@ -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 @@ -225,7 +225,7 @@ subroutine build_hycom1_target_anomaly(CS, remapCS, eqn_of_state, nz, depth, h, real, dimension(nz+1), intent(out) :: RiAnom !< The interface density anomaly !! w.r.t. the interface target !! densities [R ~> kg m-3] - real, optional, intent(in) :: h_neglect !< A negligibly small width for the purpose of + real, intent(in) :: h_neglect !< A negligibly small width for the purpose of !! cell reconstruction [H ~> m or kg m-2] real, optional, intent(in) :: h_neglect_edge !< A negligibly small width for the purpose of !! edge value calculation [H ~> m or kg m-2] diff --git a/src/ALE/coord_rho.F90 b/src/ALE/coord_rho.F90 index 3ed769f4e4..c967687dc8 100644 --- a/src/ALE/coord_rho.F90 +++ b/src/ALE/coord_rho.F90 @@ -105,7 +105,7 @@ subroutine build_rho_column(CS, nz, depth, h, T, S, eqn_of_state, z_interface, & !! units as depth) [H ~> m or kg m-2] real, optional, intent(in) :: eta_orig !< The actual original height of the top in the same !! units as depth) [H ~> m or kg m-2] - real, optional, intent(in) :: h_neglect !< A negligibly small width for the purpose + real, intent(in) :: h_neglect !< A negligibly small width for the purpose !! of cell reconstructions [H ~> m or kg m-2] real, optional, intent(in) :: h_neglect_edge !< A negligibly small width for the purpose !! of edge value calculations [H ~> m or kg m-2] @@ -201,7 +201,7 @@ subroutine build_rho_column_iteratively(CS, remapCS, nz, depth, h, T, S, eqn_of_ real, dimension(nz), intent(in) :: S !< S for column [S ~> ppt] type(EOS_type), intent(in) :: eqn_of_state !< Equation of state structure real, dimension(nz+1), intent(inout) :: zInterface !< Absolute positions of interfaces [Z ~> m] - real, optional, intent(in) :: h_neglect !< A negligibly small width for the + real, intent(in) :: h_neglect !< A negligibly small width for the !! purpose of cell reconstructions !! in the same units as h [Z ~> m] real, optional, intent(in) :: h_neglect_edge !< A negligibly small width @@ -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/ALE/regrid_edge_values.F90 b/src/ALE/regrid_edge_values.F90 index 0814c6a907..fa3bb55a43 100644 --- a/src/ALE/regrid_edge_values.F90 +++ b/src/ALE/regrid_edge_values.F90 @@ -18,15 +18,10 @@ module regrid_edge_values public edge_values_implicit_h4, edge_values_implicit_h6 public edge_slopes_implicit_h3, edge_slopes_implicit_h5 -! The following parameters are used to avoid singular matrices for boundary -! extrapolation. The are needed only in the case where thicknesses vanish +! The following parameter is used to avoid singular matrices for boundary +! extrapolation. It is needed only in the case where thicknesses vanish ! to a small enough values such that the eigenvalues of the matrix can not ! be separated. -! Specifying a dimensional parameter value, as is done here, is a terrible idea. -real, parameter :: hNeglect_edge_dflt = 1.e-10 !< The default value for cut-off minimum - !! thickness for sum(h) in edge value inversions -real, parameter :: hNeglect_dflt = 1.e-30 !< The default value for cut-off minimum - !! thickness for sum(h) in other calculations real, parameter :: hMinFrac = 1.e-5 !< A minimum fraction for min(h)/sum(h) [nondim] contains @@ -47,20 +42,16 @@ subroutine bound_edge_values( N, h, u, edge_val, h_neglect, answer_date ) real, dimension(N), intent(in) :: u !< cell average properties in arbitrary units [A] real, dimension(N,2), intent(inout) :: edge_val !< Potentially modified edge values [A]; the !! second index is for the two edges of each cell. - real, optional, intent(in) :: h_neglect !< A negligibly small width [H] + real, intent(in) :: h_neglect !< A negligibly small width [H] integer, optional, intent(in) :: answer_date !< The vintage of the expressions to use ! Local variables real :: sigma_l, sigma_c, sigma_r ! left, center and right van Leer slopes [A H-1] or [A] real :: slope_x_h ! retained PLM slope times half grid step [A] - real :: hNeglect ! A negligible thickness [H]. logical :: use_2018_answers ! If true use older, less accurate expressions. integer :: k, km1, kp1 ! Loop index and the values to either side. use_2018_answers = .true. ; if (present(answer_date)) use_2018_answers = (answer_date < 20190101) - if (use_2018_answers) then - hNeglect = hNeglect_dflt ; if (present(h_neglect)) hNeglect = h_neglect - endif ! Loop on cells to bound edge value do k = 1,N @@ -73,9 +64,9 @@ subroutine bound_edge_values( N, h, u, edge_val, h_neglect, answer_date ) slope_x_h = 0.0 if (use_2018_answers) then - sigma_l = 2.0 * ( u(k) - u(km1) ) / ( h(k) + hNeglect ) - sigma_c = 2.0 * ( u(kp1) - u(km1) ) / ( h(km1) + 2.0*h(k) + h(kp1) + hNeglect ) - sigma_r = 2.0 * ( u(kp1) - u(k) ) / ( h(k) + hNeglect ) + sigma_l = 2.0 * ( u(k) - u(km1) ) / ( h(k) + h_neglect ) + sigma_c = 2.0 * ( u(kp1) - u(km1) ) / ( h(km1) + 2.0*h(k) + h(kp1) + h_neglect ) + sigma_r = 2.0 * ( u(kp1) - u(k) ) / ( h(k) + h_neglect ) ! The limiter is used in the local coordinate system to each cell, so for convenience store ! the slope times a half grid spacing. (See White and Adcroft JCP 2008 Eqs 19 and 20) @@ -225,7 +216,7 @@ subroutine edge_values_explicit_h4( N, h, u, edge_val, h_neglect, answer_date ) real, dimension(N), intent(in) :: u !< cell average properties in arbitrary units [A] real, dimension(N,2), intent(inout) :: edge_val !< Returned edge values [A]; the second index !! is for the two edges of each cell. - real, optional, intent(in) :: h_neglect !< A negligibly small width [H] + real, intent(in) :: h_neglect !< A negligibly small width [H] integer, optional, intent(in) :: answer_date !< The vintage of the expressions to use ! Local variables @@ -248,16 +239,10 @@ subroutine edge_values_explicit_h4( N, h, u, edge_val, h_neglect, answer_date ) real, dimension(4) :: B ! The right hand side of the system to solve for C [A H] real, dimension(4) :: C ! The coefficients of a fit polynomial in units that vary ! with the index (j) as [A H^(j-1)] - real :: hNeglect ! A negligible thickness in the same units as h [H]. integer :: i, j logical :: use_2018_answers ! If true use older, less accurate expressions. use_2018_answers = .true. ; if (present(answer_date)) use_2018_answers = (answer_date < 20190101) - if (use_2018_answers) then - hNeglect = hNeglect_edge_dflt ; if (present(h_neglect)) hNeglect = h_neglect - else - hNeglect = hNeglect_dflt ; if (present(h_neglect)) hNeglect = h_neglect - endif ! Loop on interior cells do i = 3,N-1 @@ -270,9 +255,9 @@ subroutine edge_values_explicit_h4( N, h, u, edge_val, h_neglect, answer_date ) ! Avoid singularities when consecutive pairs of h vanish if (h0+h1==0.0 .or. h1+h2==0.0 .or. h2+h3==0.0) then if (use_2018_answers) then - h_min = hMinFrac*max( hNeglect, h0+h1+h2+h3 ) + h_min = hMinFrac*max( h_neglect, h0+h1+h2+h3 ) else - h_min = hMinFrac*max( hNeglect, (h0+h1)+(h2+h3) ) + h_min = hMinFrac*max( h_neglect, (h0+h1)+(h2+h3) ) endif h0 = max( h_min, h(i-2) ) h1 = max( h_min, h(i-1) ) @@ -307,7 +292,7 @@ subroutine edge_values_explicit_h4( N, h, u, edge_val, h_neglect, answer_date ) ! Determine first two edge values if (use_2018_answers) then - h_min = max( hNeglect, hMinFrac*sum(h(1:4)) ) + h_min = max( h_neglect, hMinFrac*sum(h(1:4)) ) x(1) = 0.0 do i = 1,4 dx = max(h_min, h(i) ) @@ -322,7 +307,7 @@ subroutine edge_values_explicit_h4( N, h, u, edge_val, h_neglect, answer_date ) edge_val(1,1) = evaluation_polynomial( C, 4, x(1) ) edge_val(1,2) = evaluation_polynomial( C, 4, x(2) ) else ! Use expressions with less sensitivity to roundoff - do i=1,4 ; dz(i) = max(hNeglect, h(i) ) ; u_tmp(i) = u(i) ; enddo + do i=1,4 ; dz(i) = max(h_neglect, h(i) ) ; u_tmp(i) = u(i) ; enddo call end_value_h4(dz, u_tmp, C) ! Set the edge values of the first cell @@ -333,7 +318,7 @@ subroutine edge_values_explicit_h4( N, h, u, edge_val, h_neglect, answer_date ) ! Determine two edge values of the last cell if (use_2018_answers) then - h_min = max( hNeglect, hMinFrac*sum(h(N-3:N)) ) + h_min = max( h_neglect, hMinFrac*sum(h(N-3:N)) ) x(1) = 0.0 do i = 1,4 @@ -351,7 +336,7 @@ subroutine edge_values_explicit_h4( N, h, u, edge_val, h_neglect, answer_date ) else ! Use expressions with less sensitivity to roundoff, including using a coordinate ! system that sets the origin at the last interface in the domain. - do i=1,4 ; dz(i) = max(hNeglect, h(N+1-i) ) ; u_tmp(i) = u(N+1-i) ; enddo + do i=1,4 ; dz(i) = max(h_neglect, h(N+1-i) ) ; u_tmp(i) = u(N+1-i) ; enddo call end_value_h4(dz, u_tmp, C) ! Set the last and second to last edge values @@ -384,11 +369,10 @@ subroutine edge_values_explicit_h4cw( N, h, u, edge_val, h_neglect ) real, dimension(N), intent(in) :: u !< cell average properties in arbitrary units [A] real, dimension(N,2), intent(inout) :: edge_val !< Returned edge values [A]; the second index !! is for the two edges of each cell. - real, optional, intent(in) :: h_neglect !< A negligibly small width [H] + real, intent(in) :: h_neglect !< A negligibly small width [H] ! Local variables real :: dp(N) ! Input grid layer thicknesses, but with a minimum thickness [H ~> m or kg m-2] - real :: hNeglect ! A negligible thickness in the same units as h [H] real :: da ! Difference between the unlimited scalar edge value estimates [A] real :: a6 ! Scalar field differences that are proportional to the curvature [A] real :: slk, srk ! Differences between adjacent cell averages of scalars [A] @@ -403,10 +387,8 @@ subroutine edge_values_explicit_h4cw( N, h, u, edge_val, h_neglect ) real :: h23_h122(N+1) ! A ratio of sums of adjacent thicknesses [nondim], 2/3 in the limit of uniform thicknesses. integer :: k - hNeglect = hNeglect_dflt ; if (present(h_neglect)) hNeglect = h_neglect - ! Set the thicknesses for very thin layers to some minimum value. - do k=1,N ; dp(k) = max(h(k), hNeglect) ; enddo + do k=1,N ; dp(k) = max(h(k), h_neglect) ; enddo !compute grid metrics do k=2,N @@ -494,7 +476,7 @@ subroutine edge_values_implicit_h4( N, h, u, edge_val, h_neglect, answer_date ) real, dimension(N), intent(in) :: u !< cell average properties in arbitrary units [A] real, dimension(N,2), intent(inout) :: edge_val !< Returned edge values [A]; the second index !! is for the two edges of each cell. - real, optional, intent(in) :: h_neglect !< A negligibly small width [H] + real, intent(in) :: h_neglect !< A negligibly small width [H] integer, optional, intent(in) :: answer_date !< The vintage of the expressions to use ! Local variables @@ -524,15 +506,9 @@ subroutine edge_values_implicit_h4( N, h, u, edge_val, h_neglect, answer_date ) tri_u, & ! tridiagonal system (upper diagonal) [nondim] tri_b, & ! tridiagonal system (right hand side) [A] tri_x ! tridiagonal system (solution vector) [A] - real :: hNeglect ! A negligible thickness [H] logical :: use_2018_answers ! If true use older, less accurate expressions. use_2018_answers = .true. ; if (present(answer_date)) use_2018_answers = (answer_date < 20190101) - if (use_2018_answers) then - hNeglect = hNeglect_edge_dflt ; if (present(h_neglect)) hNeglect = h_neglect - else - hNeglect = hNeglect_dflt ; if (present(h_neglect)) hNeglect = h_neglect - endif ! Loop on cells (except last one) do i = 1,N-1 @@ -542,8 +518,8 @@ subroutine edge_values_implicit_h4( N, h, u, edge_val, h_neglect, answer_date ) h1 = h(i+1) ! Avoid singularities when h0+h1=0 if (h0+h1==0.) then - h0 = hNeglect - h1 = hNeglect + h0 = h_neglect + h1 = h_neglect endif ! Auxiliary calculations @@ -562,8 +538,8 @@ subroutine edge_values_implicit_h4( N, h, u, edge_val, h_neglect, answer_date ) tri_d(i+1) = 1.0 else ! Use expressions with less sensitivity to roundoff ! Get cell widths - h0 = max(h(i), hNeglect) - h1 = max(h(i+1), hNeglect) + h0 = max(h(i), h_neglect) + h1 = max(h(i+1), h_neglect) ! The 1e-12 here attempts to balance truncation errors from the differences of ! large numbers against errors from approximating thin layers as non-vanishing. if (abs(h0) < 1.0e-12*abs(h1)) h0 = 1.0e-12*h1 @@ -587,7 +563,7 @@ subroutine edge_values_implicit_h4( N, h, u, edge_val, h_neglect, answer_date ) ! Boundary conditions: set the first boundary value if (use_2018_answers) then - h_min = max( hNeglect, hMinFrac*sum(h(1:4)) ) + h_min = max( h_neglect, hMinFrac*sum(h(1:4)) ) x(1) = 0.0 do i = 1,4 dx = max(h_min, h(i) ) @@ -601,7 +577,7 @@ subroutine edge_values_implicit_h4( N, h, u, edge_val, h_neglect, answer_date ) tri_b(1) = evaluation_polynomial( Csys, 4, x(1) ) ! Set the first edge value tri_d(1) = 1.0 else ! Use expressions with less sensitivity to roundoff - do i=1,4 ; dz(i) = max(hNeglect, h(i) ) ; u_tmp(i) = u(i) ; enddo + do i=1,4 ; dz(i) = max(h_neglect, h(i) ) ; u_tmp(i) = u(i) ; enddo call end_value_h4(dz, u_tmp, Csys) tri_b(1) = Csys(1) ! Set the first edge value. @@ -611,7 +587,7 @@ subroutine edge_values_implicit_h4( N, h, u, edge_val, h_neglect, answer_date ) ! Boundary conditions: set the last boundary value if (use_2018_answers) then - h_min = max( hNeglect, hMinFrac*sum(h(N-3:N)) ) + h_min = max( h_neglect, hMinFrac*sum(h(N-3:N)) ) x(1) = 0.0 do i=1,4 dx = max(h_min, h(N-4+i) ) @@ -629,7 +605,7 @@ subroutine edge_values_implicit_h4( N, h, u, edge_val, h_neglect, answer_date ) else ! Use expressions with less sensitivity to roundoff, including using a coordinate ! system that sets the origin at the last interface in the domain. - do i=1,4 ; dz(i) = max(hNeglect, h(N+1-i) ) ; u_tmp(i) = u(N+1-i) ; enddo + do i=1,4 ; dz(i) = max(h_neglect, h(N+1-i) ) ; u_tmp(i) = u(N+1-i) ; enddo call end_value_h4(dz, u_tmp, Csys) tri_b(N+1) = Csys(1) ! Set the last edge value @@ -806,7 +782,7 @@ subroutine edge_slopes_implicit_h3( N, h, u, edge_slopes, h_neglect, answer_date real, dimension(N), intent(in) :: u !< cell average properties in arbitrary units [A] real, dimension(N,2), intent(inout) :: edge_slopes !< Returned edge slopes [A H-1]; the !! second index is for the two edges of each cell. - real, optional, intent(in) :: h_neglect !< A negligibly small width [H] + real, intent(in) :: h_neglect !< A negligibly small width [H] integer, optional, intent(in) :: answer_date !< The vintage of the expressions to use ! Local variables @@ -837,12 +813,10 @@ subroutine edge_slopes_implicit_h3( N, h, u, edge_slopes, h_neglect, answer_date tri_u, & ! tridiagonal system (upper diagonal) [nondim] tri_b, & ! tridiagonal system (right hand side) [A H-1] tri_x ! tridiagonal system (solution vector) [A H-1] - real :: hNeglect ! A negligible thickness [H]. real :: hNeglect3 ! hNeglect^3 [H3]. logical :: use_2018_answers ! If true use older, less accurate expressions. - hNeglect = hNeglect_dflt ; if (present(h_neglect)) hNeglect = h_neglect - hNeglect3 = hNeglect**3 + hNeglect3 = h_neglect**3 use_2018_answers = .true. ; if (present(answer_date)) use_2018_answers = (answer_date < 20190101) ! Loop on cells (except last one) @@ -875,8 +849,8 @@ subroutine edge_slopes_implicit_h3( N, h, u, edge_slopes, h_neglect, answer_date tri_b(i+1) = a * u(i) + b * u(i+1) else ! Get cell widths - h0 = max(h(i), hNeglect) - h1 = max(h(i+1), hNeglect) + h0 = max(h(i), h_neglect) + h1 = max(h(i+1), h_neglect) I_h = 1.0 / (h0 + h1) h0 = h0 * I_h ; h1 = h1 * I_h @@ -917,7 +891,7 @@ subroutine edge_slopes_implicit_h3( N, h, u, edge_slopes, h_neglect, answer_date tri_b(1) = evaluation_polynomial( Dsys, 3, x(1) ) ! Set the first edge slope tri_d(1) = 1.0 else ! Use expressions with less sensitivity to roundoff - do i=1,4 ; dz(i) = max(hNeglect, h(i) ) ; u_tmp(i) = u(i) ; enddo + do i=1,4 ; dz(i) = max(h_neglect, h(i) ) ; u_tmp(i) = u(i) ; enddo call end_value_h4(dz, u_tmp, Csys) ! Set the first edge slope @@ -945,7 +919,7 @@ subroutine edge_slopes_implicit_h3( N, h, u, edge_slopes, h_neglect, answer_date else ! Use expressions with less sensitivity to roundoff, including using a coordinate ! system that sets the origin at the last interface in the domain. - do i=1,4 ; dz(i) = max(hNeglect, h(N+1-i) ) ; u_tmp(i) = u(N+1-i) ; enddo + do i=1,4 ; dz(i) = max(h_neglect, h(N+1-i) ) ; u_tmp(i) = u(N+1-i) ; enddo call end_value_h4(dz, u_tmp, Csys) @@ -980,7 +954,7 @@ subroutine edge_slopes_implicit_h5( N, h, u, edge_slopes, h_neglect, answer_date real, dimension(N), intent(in) :: u !< cell average properties in arbitrary units [A] real, dimension(N,2), intent(inout) :: edge_slopes !< Returned edge slopes [A H-1]; the !! second index is for the two edges of each cell. - real, optional, intent(in) :: h_neglect !< A negligibly small width [H] + real, intent(in) :: h_neglect !< A negligibly small width [H] integer, optional, intent(in) :: answer_date !< The vintage of the expressions to use ! ----------------------------------------------------------------------------- @@ -1021,7 +995,6 @@ subroutine edge_slopes_implicit_h5( N, h, u, edge_slopes, h_neglect, answer_date real :: hMin ! The minimum thickness used in these calculations [H] real :: h01, h01_2 ! Summed thicknesses to various powers [H^n ~> m^n or kg^n m-2n] real :: h23, h23_2 ! Summed thicknesses to various powers [H^n ~> m^n or kg^n m-2n] - real :: hNeglect ! A negligible thickness [H]. real :: h1_2, h2_2 ! Squares of thicknesses [H2] real :: h1_3, h2_3 ! Cubes of thicknesses [H3] real :: h1_4, h2_4 ! Fourth powers of thicknesses [H4] @@ -1045,12 +1018,10 @@ subroutine edge_slopes_implicit_h5( N, h, u, edge_slopes, h_neglect, answer_date real :: h_Min_Frac = 1.0e-4 ! A minimum fractional thickness [nondim] integer :: i, k ! loop indexes - hNeglect = hNeglect_dflt ; if (present(h_neglect)) hNeglect = h_neglect - ! Loop on cells (except the first and last ones) do k = 2,N-2 ! Store temporary cell widths, avoiding singularities from zero thicknesses or extreme changes. - hMin = max(hNeglect, h_Min_Frac*((h(k-1) + h(k)) + (h(k+1) + h(k+2)))) + hMin = max(h_neglect, h_Min_Frac*((h(k-1) + h(k)) + (h(k+1) + h(k+2)))) h0 = max(h(k-1), hMin) ; h1 = max(h(k), hMin) h2 = max(h(k+1), hMin) ; h3 = max(h(k+2), hMin) @@ -1091,7 +1062,7 @@ subroutine edge_slopes_implicit_h5( N, h, u, edge_slopes, h_neglect, answer_date ! Use a right-biased stencil for the second row, as described in Eq. (53) of White and Adcroft (2009). ! Store temporary cell widths, avoiding singularities from zero thicknesses or extreme changes. - hMin = max(hNeglect, h_Min_Frac*((h(1) + h(2)) + (h(3) + h(4)))) + hMin = max(h_neglect, h_Min_Frac*((h(1) + h(2)) + (h(3) + h(4)))) h0 = max(h(1), hMin) ; h1 = max(h(2), hMin) h2 = max(h(3), hMin) ; h3 = max(h(4), hMin) @@ -1147,7 +1118,7 @@ subroutine edge_slopes_implicit_h5( N, h, u, edge_slopes, h_neglect, answer_date ! Use a left-biased stencil for the second to last row, as described in Eq. (54) of White and Adcroft (2009). ! Store temporary cell widths, avoiding singularities from zero thicknesses or extreme changes. - hMin = max(hNeglect, h_Min_Frac*((h(N-3) + h(N-2)) + (h(N-1) + h(N)))) + hMin = max(h_neglect, h_Min_Frac*((h(N-3) + h(N-2)) + (h(N-1) + h(N)))) h0 = max(h(N-3), hMin) ; h1 = max(h(N-2), hMin) h2 = max(h(N-1), hMin) ; h3 = max(h(N), hMin) @@ -1255,7 +1226,7 @@ subroutine edge_values_implicit_h6( N, h, u, edge_val, h_neglect, answer_date ) real, dimension(N), intent(in) :: u !< cell average properties (size N) in arbitrary units [A] real, dimension(N,2), intent(inout) :: edge_val !< Returned edge values [A]; the second index !! is for the two edges of each cell. - real, optional, intent(in) :: h_neglect !< A negligibly small width [H] + real, intent(in) :: h_neglect !< A negligibly small width [H] integer, optional, intent(in) :: answer_date !< The vintage of the expressions to use ! Local variables @@ -1263,7 +1234,6 @@ subroutine edge_values_implicit_h6( N, h, u, edge_val, h_neglect, answer_date ) real :: hMin ! The minimum thickness used in these calculations [H] real :: h01, h01_2, h01_3 ! Summed thicknesses to various powers [H^n ~> m^n or kg^n m-2n] real :: h23, h23_2, h23_3 ! Summed thicknesses to various powers [H^n ~> m^n or kg^n m-2n] - real :: hNeglect ! A negligible thickness [H]. real :: h1_2, h2_2, h1_3, h2_3 ! Cell widths raised to the 2nd and 3rd powers [H2] or [H3] real :: h1_4, h2_4, h1_5, h2_5 ! Cell widths raised to the 4th and 5th powers [H4] or [H5] real :: alpha, beta ! stencil coefficients [nondim] @@ -1286,12 +1256,10 @@ subroutine edge_values_implicit_h6( N, h, u, edge_val, h_neglect, answer_date ) tri_x ! trid. system (unknowns vector) [A] integer :: i, k ! loop indexes - hNeglect = hNeglect_edge_dflt ; if (present(h_neglect)) hNeglect = h_neglect - ! Loop on interior cells do k = 2,N-2 ! Store temporary cell widths, avoiding singularities from zero thicknesses or extreme changes. - hMin = max(hNeglect, hMinFrac*((h(k-1) + h(k)) + (h(k+1) + h(k+2)))) + hMin = max(h_neglect, hMinFrac*((h(k-1) + h(k)) + (h(k+1) + h(k+2)))) h0 = max(h(k-1), hMin) ; h1 = max(h(k), hMin) h2 = max(h(k+1), hMin) ; h3 = max(h(k+2), hMin) @@ -1329,7 +1297,7 @@ subroutine edge_values_implicit_h6( N, h, u, edge_val, h_neglect, answer_date ) ! Use a right-biased stencil for the second row, as described in Eq. (49) of White and Adcroft (2009). ! Store temporary cell widths, avoiding singularities from zero thicknesses or extreme changes. - hMin = max(hNeglect, hMinFrac*((h(1) + h(2)) + (h(3) + h(4)))) + hMin = max(h_neglect, hMinFrac*((h(1) + h(2)) + (h(3) + h(4)))) h0 = max(h(1), hMin) ; h1 = max(h(2), hMin) h2 = max(h(3), hMin) ; h3 = max(h(4), hMin) @@ -1364,7 +1332,7 @@ subroutine edge_values_implicit_h6( N, h, u, edge_val, h_neglect, answer_date ) tri_b(2) = Csys(3) * u(1) + Csys(4) * u(2) + Csys(5) * u(3) + Csys(6) * u(4) ! Boundary conditions: left boundary - hMin = max( hNeglect, hMinFrac*((h(1)+h(2)) + (h(5)+h(6)) + (h(3)+h(4))) ) + hMin = max( h_neglect, hMinFrac*((h(1)+h(2)) + (h(5)+h(6)) + (h(3)+h(4))) ) x(1) = 0.0 do i = 1,6 dx = max( hMin, h(i) ) @@ -1386,7 +1354,7 @@ subroutine edge_values_implicit_h6( N, h, u, edge_val, h_neglect, answer_date ) ! Use a left-biased stencil for the second to last row, as described in Eq. (50) of White and Adcroft (2009). ! Store temporary cell widths, avoiding singularities from zero thicknesses or extreme changes. - hMin = max(hNeglect, hMinFrac*((h(N-3) + h(N-2)) + (h(N-1) + h(N)))) + hMin = max(h_neglect, hMinFrac*((h(N-3) + h(N-2)) + (h(N-1) + h(N)))) h0 = max(h(N-3), hMin) ; h1 = max(h(N-2), hMin) h2 = max(h(N-1), hMin) ; h3 = max(h(N), hMin) @@ -1421,7 +1389,7 @@ subroutine edge_values_implicit_h6( N, h, u, edge_val, h_neglect, answer_date ) tri_b(N) = Csys(3) * u(N-3) + Csys(4) * u(N-2) + Csys(5) * u(N-1) + Csys(6) * u(N) ! Boundary conditions: right boundary - hMin = max( hNeglect, hMinFrac*(h(N-3) + h(N-2)) + ((h(N-1) + h(N)) + (h(N-5) + h(N-4))) ) + hMin = max( h_neglect, hMinFrac*(h(N-3) + h(N-2)) + ((h(N-1) + h(N)) + (h(N-5) + h(N-4))) ) x(1) = 0.0 do i = 1,6 dx = max( hMin, h(N+1-i) ) diff --git a/src/ALE/regrid_interp.F90 b/src/ALE/regrid_interp.F90 index b3100fe8ae..6e0be9ebba 100644 --- a/src/ALE/regrid_interp.F90 +++ b/src/ALE/regrid_interp.F90 @@ -87,15 +87,19 @@ subroutine regridding_set_ppolys(CS, densities, n0, h0, ppoly0_E, ppoly0_S, & real, dimension(n0,2), intent(inout) :: ppoly0_S !< Edge slope of polynomial [A H-1] real, dimension(n0,DEGREE_MAX+1), intent(inout) :: ppoly0_coefs !< Coefficients of polynomial [A] integer, intent(inout) :: degree !< The degree of the polynomials - real, optional, intent(in) :: h_neglect !< A negligibly small width for the + real, intent(in) :: h_neglect !< A negligibly small width for the !! purpose of cell reconstructions [H] !! in the same units as h0. real, optional, intent(in) :: h_neglect_edge !< A negligibly small width !! for the purpose of edge value calculations [H] !! in the same units as h0. ! Local variables + real :: h_neg_edge ! A negligibly small width for the purpose of edge value + ! calculations in the same units as h0 [H] logical :: extrapolate + h_neg_edge = h_neglect ; if (present(h_neglect_edge)) h_neg_edge = h_neglect_edge + ! Reset piecewise polynomials ppoly0_E(:,:) = 0.0 ppoly0_S(:,:) = 0.0 @@ -117,7 +121,7 @@ subroutine regridding_set_ppolys(CS, densities, n0, h0, ppoly0_E, ppoly0_S, & case ( INTERPOLATION_P1M_H4 ) degree = DEGREE_1 if ( n0 >= 4 ) then - call edge_values_explicit_h4( n0, h0, densities, ppoly0_E, h_neglect_edge, answer_date=CS%answer_date ) + call edge_values_explicit_h4( n0, h0, densities, ppoly0_E, h_neg_edge, answer_date=CS%answer_date ) else call edge_values_explicit_h2( n0, h0, densities, ppoly0_E ) endif @@ -129,7 +133,7 @@ subroutine regridding_set_ppolys(CS, densities, n0, h0, ppoly0_E, ppoly0_S, & case ( INTERPOLATION_P1M_IH4 ) degree = DEGREE_1 if ( n0 >= 4 ) then - call edge_values_implicit_h4( n0, h0, densities, ppoly0_E, h_neglect_edge, answer_date=CS%answer_date ) + call edge_values_implicit_h4( n0, h0, densities, ppoly0_E, h_neg_edge, answer_date=CS%answer_date ) else call edge_values_explicit_h2( n0, h0, densities, ppoly0_E ) endif @@ -148,7 +152,7 @@ subroutine regridding_set_ppolys(CS, densities, n0, h0, ppoly0_E, ppoly0_S, & case ( INTERPOLATION_PPM_CW ) if ( n0 >= 4 ) then degree = DEGREE_2 - call edge_values_explicit_h4cw( n0, h0, densities, ppoly0_E, h_neglect_edge ) + call edge_values_explicit_h4cw( n0, h0, densities, ppoly0_E, h_neg_edge ) call PPM_monotonicity( n0, densities, ppoly0_E ) call PPM_reconstruction( n0, h0, densities, ppoly0_E, ppoly0_coefs, h_neglect, answer_date=CS%answer_date ) if (extrapolate) then @@ -167,7 +171,7 @@ subroutine regridding_set_ppolys(CS, densities, n0, h0, ppoly0_E, ppoly0_S, & case ( INTERPOLATION_PPM_H4 ) if ( n0 >= 4 ) then degree = DEGREE_2 - call edge_values_explicit_h4( n0, h0, densities, ppoly0_E, h_neglect_edge, answer_date=CS%answer_date ) + call edge_values_explicit_h4( n0, h0, densities, ppoly0_E, h_neg_edge, answer_date=CS%answer_date ) call PPM_reconstruction( n0, h0, densities, ppoly0_E, ppoly0_coefs, h_neglect, answer_date=CS%answer_date ) if (extrapolate) then call PPM_boundary_extrapolation( n0, h0, densities, ppoly0_E, & @@ -185,7 +189,7 @@ subroutine regridding_set_ppolys(CS, densities, n0, h0, ppoly0_E, ppoly0_S, & case ( INTERPOLATION_PPM_IH4 ) if ( n0 >= 4 ) then degree = DEGREE_2 - call edge_values_implicit_h4( n0, h0, densities, ppoly0_E, h_neglect_edge, answer_date=CS%answer_date ) + call edge_values_implicit_h4( n0, h0, densities, ppoly0_E, h_neg_edge, answer_date=CS%answer_date ) call PPM_reconstruction( n0, h0, densities, ppoly0_E, ppoly0_coefs, h_neglect, answer_date=CS%answer_date ) if (extrapolate) then call PPM_boundary_extrapolation( n0, h0, densities, ppoly0_E, & @@ -203,13 +207,13 @@ subroutine regridding_set_ppolys(CS, densities, n0, h0, ppoly0_E, ppoly0_S, & case ( INTERPOLATION_P3M_IH4IH3 ) if ( n0 >= 4 ) then degree = DEGREE_3 - call edge_values_implicit_h4( n0, h0, densities, ppoly0_E, h_neglect_edge, answer_date=CS%answer_date ) + call edge_values_implicit_h4( n0, h0, densities, ppoly0_E, h_neg_edge, answer_date=CS%answer_date ) call edge_slopes_implicit_h3( n0, h0, densities, ppoly0_S, h_neglect, answer_date=CS%answer_date ) call P3M_interpolation( n0, h0, densities, ppoly0_E, ppoly0_S, & ppoly0_coefs, h_neglect, answer_date=CS%answer_date ) if (extrapolate) then call P3M_boundary_extrapolation( n0, h0, densities, ppoly0_E, ppoly0_S, & - ppoly0_coefs, h_neglect, h_neglect_edge ) + ppoly0_coefs, h_neglect, h_neg_edge ) endif else degree = DEGREE_1 @@ -223,7 +227,7 @@ subroutine regridding_set_ppolys(CS, densities, n0, h0, ppoly0_E, ppoly0_S, & case ( INTERPOLATION_P3M_IH6IH5 ) if ( n0 >= 6 ) then degree = DEGREE_3 - call edge_values_implicit_h6( n0, h0, densities, ppoly0_E, h_neglect_edge, answer_date=CS%answer_date ) + call edge_values_implicit_h6( n0, h0, densities, ppoly0_E, h_neg_edge, answer_date=CS%answer_date ) call edge_slopes_implicit_h5( n0, h0, densities, ppoly0_S, h_neglect, answer_date=CS%answer_date ) call P3M_interpolation( n0, h0, densities, ppoly0_E, ppoly0_S, & ppoly0_coefs, h_neglect, answer_date=CS%answer_date ) @@ -243,7 +247,7 @@ subroutine regridding_set_ppolys(CS, densities, n0, h0, ppoly0_E, ppoly0_S, & case ( INTERPOLATION_PQM_IH4IH3 ) if ( n0 >= 4 ) then degree = DEGREE_4 - call edge_values_implicit_h4( n0, h0, densities, ppoly0_E, h_neglect_edge, answer_date=CS%answer_date ) + call edge_values_implicit_h4( n0, h0, densities, ppoly0_E, h_neg_edge, answer_date=CS%answer_date ) call edge_slopes_implicit_h3( n0, h0, densities, ppoly0_S, h_neglect, answer_date=CS%answer_date ) call PQM_reconstruction( n0, h0, densities, ppoly0_E, ppoly0_S, & ppoly0_coefs, h_neglect, answer_date=CS%answer_date ) @@ -263,7 +267,7 @@ subroutine regridding_set_ppolys(CS, densities, n0, h0, ppoly0_E, ppoly0_S, & case ( INTERPOLATION_PQM_IH6IH5 ) if ( n0 >= 6 ) then degree = DEGREE_4 - call edge_values_implicit_h6( n0, h0, densities, ppoly0_E, h_neglect_edge, answer_date=CS%answer_date ) + call edge_values_implicit_h6( n0, h0, densities, ppoly0_E, h_neg_edge, answer_date=CS%answer_date ) call edge_slopes_implicit_h5( n0, h0, densities, ppoly0_S, h_neglect, answer_date=CS%answer_date ) call PQM_reconstruction( n0, h0, densities, ppoly0_E, ppoly0_S, & ppoly0_coefs, h_neglect, answer_date=CS%answer_date ) @@ -335,7 +339,7 @@ subroutine build_and_interpolate_grid(CS, densities, n0, h0, x0, target_values, real, dimension(n0+1), intent(in) :: x0 !< Source interface positions [H] real, dimension(n1), intent(inout) :: h1 !< Output cell widths [H] real, dimension(n1+1), intent(inout) :: x1 !< Target interface positions [H] - real, optional, intent(in) :: h_neglect !< A negligibly small width for the + real, intent(in) :: h_neglect !< A negligibly small width for the !! purpose of cell reconstructions [H] !! in the same units as h0. real, optional, intent(in) :: h_neglect_edge !< A negligibly small width diff --git a/src/ALE/remapping_attic.F90 b/src/ALE/remapping_attic.F90 index be20a27466..ab345dc53e 100644 --- a/src/ALE/remapping_attic.F90 +++ b/src/ALE/remapping_attic.F90 @@ -28,11 +28,6 @@ module remapping_attic ! outside of the range 0 to 1. #define __USE_ROUNDOFF_SAFE_ADJUSTMENTS__ -real, parameter :: hNeglect_dflt = 1.E-30 !< A thickness [H ~> m or kg m-2] that can be - !! added to thicknesses in a denominator without - !! changing the numerical result, except where - !! a division by zero would otherwise occur. - contains !> Compare two summation estimates of positive data and judge if due to more @@ -83,7 +78,7 @@ subroutine remapByProjection( n0, h0, u0, ppoly0_E, ppoly0_coefs, & real, intent(in) :: h1(:) !< Target grid widths (size n1) [H] integer, intent(in) :: method !< Remapping scheme to use real, intent(out) :: u1(:) !< Target cell averages (size n1) [A] - real, optional, intent(in) :: h_neglect !< A negligibly small width for the + real, intent(in) :: h_neglect !< A negligibly small width for the !! purpose of cell reconstructions !! in the same units as h [H]. ! Local variables @@ -132,7 +127,7 @@ subroutine remapByDeltaZ( n0, h0, u0, ppoly0_E, ppoly0_coefs, n1, dx1, & real, dimension(:), intent(out) :: u1 !< Target cell averages (size n1) [A] real, dimension(:), & optional, intent(out) :: h1 !< Target grid widths (size n1) [H] - real, optional, intent(in) :: h_neglect !< A negligibly small width for the + real, intent(in) :: h_neglect !< A negligibly small width for the !! purpose of cell reconstructions !! in the same units as h [H]. ! Local variables @@ -181,7 +176,7 @@ subroutine remapByDeltaZ( n0, h0, u0, ppoly0_E, ppoly0_coefs, n1, dx1, & ! hFlux is the positive width of the remapped volume hFlux = abs(dx1(iTarget+1)) call integrateReconOnInterval( n0, h0, u0, ppoly0_E, ppoly0_coefs, method, & - xL, xR, hFlux, uAve, jStart, xStart ) + xL, xR, hFlux, uAve, jStart, xStart, h_neglect ) ! uAve is the average value of u, independent of sign of dx1 fluxR = dx1(iTarget+1)*uAve ! Includes sign of dx1 @@ -218,7 +213,7 @@ subroutine integrateReconOnInterval( n0, h0, u0, ppoly0_E, ppoly0_coefs, method, !< On exit, contains index of last cell used real, intent(inout) :: xStart !< The left edge position of cell jStart [H] !< On first entry should be 0. - real, optional, intent(in) :: h_neglect !< A negligibly small width for the + real, intent(in) :: h_neglect !< A negligibly small width for the !! purpose of cell reconstructions !! in the same units as h [H] ! Local variables @@ -232,11 +227,8 @@ subroutine integrateReconOnInterval( n0, h0, u0, ppoly0_E, ppoly0_coefs, method, ! (notionally xR - xL) which differs due to roundoff [H]. real :: x0_2, x1_2 ! Squares of normalized positions used to evaluate polynomials [nondim] real :: x0px1, x02px12 ! Sums of normalized positions and their squares [nondim] - real :: hNeglect ! A negligible thickness in the same units as h [H] real, parameter :: r_3 = 1.0/3.0 ! Used in evaluation of integrated polynomials [nondim] - hNeglect = hNeglect_dflt ; if (present(h_neglect)) hNeglect = h_neglect - q = -1.E30 x0jLl = -1.E30 x0jRl = -1.E30 @@ -288,7 +280,7 @@ subroutine integrateReconOnInterval( n0, h0, u0, ppoly0_E, ppoly0_coefs, method, uAve = 0.5 * ( ppoly0_E(jL,1) + ppoly0_E(jL,2) ) else ! WHY IS THIS NOT WRITTEN AS xi0 = ( xL - x0jLl ) / h0(jL) ---AJA - xi0 = xL / ( h0(jL) + hNeglect ) - x0jLl / ( h0(jL) + hNeglect ) + xi0 = xL / ( h0(jL) + h_neglect ) - x0jLl / ( h0(jL) + h_neglect ) select case ( method ) case ( INTEGRATION_PCM ) @@ -347,11 +339,11 @@ subroutine integrateReconOnInterval( n0, h0, u0, ppoly0_E, ppoly0_coefs, method, ! ! Determine normalized coordinates #ifdef __USE_ROUNDOFF_SAFE_ADJUSTMENTS__ - xi0 = max( 0., min( 1., ( xL - x0jLl ) / ( h0(jL) + hNeglect ) ) ) - xi1 = max( 0., min( 1., ( xR - x0jLl ) / ( h0(jL) + hNeglect ) ) ) + xi0 = max( 0., min( 1., ( xL - x0jLl ) / ( h0(jL) + h_neglect ) ) ) + xi1 = max( 0., min( 1., ( xR - x0jLl ) / ( h0(jL) + h_neglect ) ) ) #else - xi0 = xL / h0(jL) - x0jLl / ( h0(jL) + hNeglect ) - xi1 = xR / h0(jL) - x0jLl / ( h0(jL) + hNeglect ) + xi0 = xL / h0(jL) - x0jLl / ( h0(jL) + h_neglect ) + xi1 = xR / h0(jL) - x0jLl / ( h0(jL) + h_neglect ) #endif hAct = h0(jL) * ( xi1 - xi0 ) @@ -403,9 +395,9 @@ subroutine integrateReconOnInterval( n0, h0, u0, ppoly0_E, ppoly0_coefs, method, ! Integrate from xL up to right boundary of cell jL #ifdef __USE_ROUNDOFF_SAFE_ADJUSTMENTS__ - xi0 = max( 0., min( 1., ( xL - x0jLl ) / ( h0(jL) + hNeglect ) ) ) + xi0 = max( 0., min( 1., ( xL - x0jLl ) / ( h0(jL) + h_neglect ) ) ) #else - xi0 = (xL - x0jLl) / ( h0(jL) + hNeglect ) + xi0 = (xL - x0jLl) / ( h0(jL) + h_neglect ) #endif xi1 = 1.0 @@ -449,9 +441,9 @@ subroutine integrateReconOnInterval( n0, h0, u0, ppoly0_E, ppoly0_coefs, method, ! Integrate from left boundary of cell jR up to xR xi0 = 0.0 #ifdef __USE_ROUNDOFF_SAFE_ADJUSTMENTS__ - xi1 = max( 0., min( 1., ( xR - x0jRl ) / ( h0(jR) + hNeglect ) ) ) + xi1 = max( 0., min( 1., ( xR - x0jRl ) / ( h0(jR) + h_neglect ) ) ) #else - xi1 = (xR - x0jRl) / ( h0(jR) + hNeglect ) + xi1 = (xR - x0jRl) / ( h0(jR) + h_neglect ) #endif hAct = hAct + h0(jR) * ( xi1 - xi0 ) @@ -568,8 +560,8 @@ logical function remapping_attic_unit_tests(verbose) v = verbose answer_date = 20190101 ! 20181231 - h_neglect = hNeglect_dflt - h_neglect_edge = hNeglect_dflt ; if (answer_date < 20190101) h_neglect_edge = 1.0e-10 + h_neglect = 1.0E-30 + h_neglect_edge = h_neglect ; if (answer_date < 20190101) h_neglect_edge = 1.0e-10 write(stdout,*) '==== remapping_attic: remapping_attic_unit_tests =================' remapping_attic_unit_tests = .false. ! Normally return false 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..c8a14de773 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 @@ -346,7 +347,10 @@ module MOM_open_boundary real :: rx_max !< The maximum magnitude of the baroclinic radiation velocity (or speed of !! characteristics) in units of grid points per timestep [nondim]. logical :: OBC_pe !< Is there an open boundary on this tile? - type(remapping_CS), pointer :: remap_CS=> NULL() !< ALE remapping control structure for segments only + type(remapping_CS), pointer :: remap_z_CS=> NULL() !< ALE remapping control structure for + !! z-space data on segments + type(remapping_CS), pointer :: remap_h_CS=> NULL() !< ALE remapping control structure for + !! thickness-based fields on segments type(OBC_registry_type), pointer :: OBC_Reg => NULL() !< Registry type for boundaries real, allocatable :: rx_normal(:,:,:) !< Array storage for normal phase speed for EW radiation OBCs in units of !! grid points per timestep [nondim] @@ -382,6 +386,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 +434,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 @@ -667,23 +675,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 +704,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 +731,41 @@ subroutine open_boundary_config(G, US, param_file, OBC) end subroutine open_boundary_config +!> Setup vertical remapping 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 in vertical height units [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_z_CS) + call initialize_remapping(OBC%remap_z_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) + allocate(OBC%remap_h_CS) + call initialize_remapping(OBC%remap_h_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=GV%H_subroundoff, h_neglect_edge=GV%H_subroundoff) + 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) @@ -1973,6 +2010,8 @@ subroutine open_boundary_dealloc(OBC) if (allocated(OBC%cff_normal_v)) deallocate(OBC%cff_normal_v) if (allocated(OBC%tres_x)) deallocate(OBC%tres_x) if (allocated(OBC%tres_y)) deallocate(OBC%tres_y) + if (associated(OBC%remap_z_CS)) deallocate(OBC%remap_z_CS) + if (associated(OBC%remap_h_CS)) deallocate(OBC%remap_h_CS) deallocate(OBC) end subroutine open_boundary_dealloc @@ -3867,7 +3906,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 +3918,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) @@ -4176,25 +4206,22 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) segment%field(m)%buffer_dst(I,J,:) = 0.0 ! initialize remap destination buffer if (G%mask2dCu(I,j)>0. .and. G%mask2dCu(I,j+1)>0.) then dz_stack(:) = 0.5*(dz(i+ishift,j,:) + dz(i+ishift,j+1,:)) - call remapping_core_h(OBC%remap_CS, & + call remapping_core_h(OBC%remap_z_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, & + call remapping_core_h(OBC%remap_z_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, & + call remapping_core_h(OBC%remap_z_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 @@ -4206,11 +4233,10 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) net_dz_src = sum( segment%field(m)%dz_src(I,j,:) ) net_dz_int = sum( dz(i+ishift,j,:) ) scl_fac = net_dz_int / net_dz_src - call remapping_core_h(OBC%remap_CS, & + call remapping_core_h(OBC%remap_z_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 @@ -4226,25 +4252,22 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) ! Using the h remapping approach ! Pretty sure we need to check for source/target grid consistency here dz_stack(:) = 0.5*(dz(i,j+jshift,:) + dz(i+1,j+jshift,:)) - call remapping_core_h(OBC%remap_CS, & + call remapping_core_h(OBC%remap_z_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, & + call remapping_core_h(OBC%remap_z_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, & + call remapping_core_h(OBC%remap_z_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 @@ -4256,11 +4279,10 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) net_dz_src = sum( segment%field(m)%dz_src(i,J,:) ) net_dz_int = sum( dz(i,j+jshift,:) ) scl_fac = net_dz_int / net_dz_src - call remapping_core_h(OBC%remap_CS, & + call remapping_core_h(OBC%remap_z_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 @@ -5528,7 +5550,6 @@ subroutine remap_OBC_fields(G, GV, h_old, h_new, OBC, PCM_cell) real :: h2(GV%ke) ! A column of target grid layer thicknesses [H ~> m or kg m-2] real :: I_scale ! The inverse of the scaling factor for the tracers. ! For salinity the units would be [ppt S-1 ~> 1]. - real :: h_neglect ! Tiny thickness used in remapping [H ~> m or kg m-2] logical :: PCM(GV%ke) ! If true, do PCM remapping from a cell. integer :: i, j, k, m, n, ntr, nz @@ -5536,7 +5557,6 @@ subroutine remap_OBC_fields(G, GV, h_old, h_new, OBC, PCM_cell) nz = GV%ke ntr = OBC%ntr - h_neglect = GV%H_subroundoff if (.not.present(PCM_cell)) PCM(:) = .false. @@ -5566,11 +5586,10 @@ subroutine remap_OBC_fields(G, GV, h_old, h_new, OBC, PCM_cell) I_scale = 1.0 ; if (segment%tr_Reg%Tr(m)%scale /= 0.0) I_scale = 1.0 / segment%tr_Reg%Tr(m)%scale 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) + call remapping_core_h(OBC%remap_h_CS, nz, h1, segment%tr_Reg%Tr(m)%tres(I,j,:), nz, h2, tr_column, & + 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_h_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? @@ -5584,8 +5603,8 @@ subroutine remap_OBC_fields(G, GV, h_old, h_new, OBC, PCM_cell) endif ; enddo 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) + call remapping_core_h(OBC%remap_h_CS, nz, h1, segment%rx_norm_rad(I,j,:), nz, h2, r_norm_col, & + PCM_cell=PCM) do k=1,nz segment%rx_norm_rad(I,j,k) = r_norm_col(k) @@ -5594,14 +5613,14 @@ subroutine remap_OBC_fields(G, GV, h_old, h_new, OBC, PCM_cell) endif 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) + call remapping_core_h(OBC%remap_h_CS, nz, h1, segment%rx_norm_obl(I,j,:), nz, h2, rxy_col, & + 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) + call remapping_core_h(OBC%remap_h_CS, nz, h1, segment%ry_norm_obl(I,j,:), nz, h2, rxy_col, & + 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) + call remapping_core_h(OBC%remap_h_CS, nz, h1, segment%cff_normal(I,j,:), nz, h2, rxy_col, & + PCM_cell=PCM) segment%cff_normal(I,j,:) = rxy_col(:) do k=1,nz @@ -5634,11 +5653,10 @@ subroutine remap_OBC_fields(G, GV, h_old, h_new, OBC, PCM_cell) I_scale = 1.0 ; if (segment%tr_Reg%Tr(m)%scale /= 0.0) I_scale = 1.0 / segment%tr_Reg%Tr(m)%scale 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) + call remapping_core_h(OBC%remap_h_CS, nz, h1, segment%tr_Reg%Tr(m)%tres(i,J,:), nz, h2, tr_column, & + 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_h_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? @@ -5652,8 +5670,8 @@ subroutine remap_OBC_fields(G, GV, h_old, h_new, OBC, PCM_cell) endif ; enddo 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) + call remapping_core_h(OBC%remap_h_CS, nz, h1, segment%ry_norm_rad(i,J,:), nz, h2, r_norm_col, & + PCM_cell=PCM) do k=1,nz segment%ry_norm_rad(i,J,k) = r_norm_col(k) @@ -5662,14 +5680,14 @@ subroutine remap_OBC_fields(G, GV, h_old, h_new, OBC, PCM_cell) endif 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) + call remapping_core_h(OBC%remap_h_CS, nz, h1, segment%rx_norm_obl(i,J,:), nz, h2, rxy_col, & + 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) + call remapping_core_h(OBC%remap_h_CS, nz, h1, segment%ry_norm_obl(i,J,:), nz, h2, rxy_col, & + 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) + call remapping_core_h(OBC%remap_h_CS, nz, h1, segment%cff_normal(i,J,:), nz, h2, rxy_col, & + PCM_cell=PCM) segment%cff_normal(i,J,:) = rxy_col(:) do k=1,nz @@ -5861,10 +5879,14 @@ subroutine rotate_OBC_config(OBC_in, G_in, OBC, G, turns) OBC%rx_max = OBC_in%rx_max OBC%OBC_pe = OBC_in%OBC_pe - ! remap_CS is set up by initialize_segment_data, so we copy the fields here. - if (ASSOCIATED(OBC_in%remap_CS)) then - allocate(OBC%remap_CS) - OBC%remap_CS = OBC_in%remap_CS + ! remap_z_CS and remap_h_CS are set up by initialize_segment_data, so we copy the fields here. + if (ASSOCIATED(OBC_in%remap_z_CS)) then + allocate(OBC%remap_z_CS) + OBC%remap_z_CS = OBC_in%remap_z_CS + endif + if (ASSOCIATED(OBC_in%remap_h_CS)) then + allocate(OBC%remap_h_CS) + OBC%remap_h_CS = OBC_in%remap_h_CS endif ! TODO: The OBC registry seems to be a list of "registered" OBC types. 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..1c508ec490 100644 --- a/src/diagnostics/MOM_wave_speed.F90 +++ b/src/diagnostics/MOM_wave_speed.F90 @@ -51,7 +51,9 @@ module MOM_wave_speed !! are simply reported as 0 [L T-1 ~> m s-1]. A non-negative !! value must be specified via a call to wave_speed_init for !! the subroutine wave_speeds to be used (but not wave_speed). - type(remapping_CS) :: remapping_CS !< Used for vertical remapping when calculating equivalent barotropic + type(remapping_CS) :: remap_2018_CS !< Used for vertical remapping when calculating equivalent barotropic + !! mode structure for answer dates below 20190101. + type(remapping_CS) :: remap_CS !< Used for vertical remapping when calculating equivalent barotropic !! mode structure. integer :: remap_answer_date = 99991231 !< The vintage of the order of arithmetic and expressions to use !! for remapping. Values below 20190101 recover the remapping @@ -674,13 +676,11 @@ subroutine wave_speed(h, tv, G, GV, US, cg1, CS, halo_size, use_ebt_mode, mono_N 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) + call remapping_core_h(CS%remap_2018_CS, kc, Hc(:), mode_struct, & + nz, h(i,j,:), modal_structure(i,j,:)) 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) + call remapping_core_h(CS%remap_CS, kc, Hc(:), mode_struct, & + nz, h(i,j,:), modal_structure(i,j,:)) endif endif else @@ -1357,9 +1357,8 @@ subroutine wave_speeds(h, tv, G, GV, US, nmodes, cn, CS, w_struct, u_struct, u_s nz, h(i,j,:), modal_structure(:), .false.) ! 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) + call remapping_core_h(CS%remap_CS, kc, Hc(1:kc), mode_struct_fder(1:kc), & + nz, h(i,j,:), modal_structure_fder(:)) ! write the wave structure do k=1,nz+1 @@ -1533,9 +1532,8 @@ subroutine wave_speeds(h, tv, G, GV, US, nmodes, cn, CS, w_struct, u_struct, u_s nz, h(i,j,:), modal_structure(:), .false.) ! 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) + call remapping_core_h(CS%remap_CS, kc, Hc(1:kc), mode_struct_fder(1:kc), & + nz, h(i,j,:), modal_structure_fder(:)) ! write the wave structure ! note that m=1 solves for 2nd mode,... @@ -1610,10 +1608,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 @@ -1657,10 +1656,18 @@ subroutine wave_speed_init(CS, use_ebt_mode, mono_N2_column_fraction, mono_N2_de remap_answers_2018=remap_answers_2018, remap_answer_date=remap_answer_date, & 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., & + ! The following remapping is only used for wave_speed with pre-2019 answers. + if (CS%remap_answer_date < 20190101) & + call initialize_remapping(CS%remap_2018_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) + + ! This is used in wave_speeds in all cases, and in wave_speed with newer answers. + call initialize_remapping(CS%remap_CS, 'PLM', boundary_extrapolation=.false., & om4_remap_via_sub_cells=om4_remap_via_sub_cells, & - answer_date=CS%remap_answer_date) + answer_date=CS%remap_answer_date, & + h_neglect=GV%H_subroundoff, h_neglect_edge=GV%H_subroundoff) 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_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index 1467cdaaad..51996f896c 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -92,7 +92,7 @@ module MOM_state_initialization use MOM_ALE, only : ALE_remap_scalar, ALE_regrid_accelerated, TS_PLM_edge_values use MOM_regridding, only : regridding_CS, set_regrid_params, getCoordinateResolution use MOM_regridding, only : regridding_main, regridding_preadjust_reqs, convective_adjustment -use MOM_regridding, only : set_dz_neglect +use MOM_regridding, only : set_dz_neglect, set_h_neglect use MOM_remapping, only : remapping_CS, initialize_remapping, remapping_core_h use MOM_horizontal_regridding, only : horiz_interp_and_extrap_tracer, homogenize_field use MOM_oda_incupd, only: oda_incupd_CS, initialize_oda_incupd_fixed, initialize_oda_incupd @@ -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 @@ -1204,10 +1210,9 @@ subroutine trim_for_ice(PF, G, GV, US, ALE_CSp, tv, h, just_read) endif do j=G%jsc,G%jec ; do i=G%isc,G%iec - call cut_off_column_top(GV%ke, tv, GV, US, GV%g_Earth, G%bathyT(i,j)+G%Z_ref, & - min_thickness, tv%T(i,j,:), T_t(i,j,:), T_b(i,j,:), & - tv%S(i,j,:), S_t(i,j,:), S_b(i,j,:), p_surf(i,j), h(i,j,:), remap_CS, & - z_tol=z_tolerance, remap_answer_date=remap_answer_date) + call cut_off_column_top(GV%ke, tv, GV, US, GV%g_Earth, G%bathyT(i,j)+G%Z_ref, min_thickness, & + tv%T(i,j,:), T_t(i,j,:), T_b(i,j,:), tv%S(i,j,:), S_t(i,j,:), S_b(i,j,:), & + p_surf(i,j), h(i,j,:), remap_CS, z_tol=z_tolerance) enddo ; enddo end subroutine trim_for_ice @@ -1298,7 +1303,7 @@ end subroutine calc_sfc_displacement !> Adjust the layer thicknesses by removing the top of the water column above the !! depth where the hydrostatic pressure matches p_surf subroutine cut_off_column_top(nk, tv, GV, US, G_earth, depth, min_thickness, T, T_t, T_b, & - S, S_t, S_b, p_surf, h, remap_CS, z_tol, remap_answer_date) + S, S_t, S_b, p_surf, h, remap_CS, z_tol) integer, intent(in) :: nk !< Number of layers type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. @@ -1318,10 +1323,6 @@ subroutine cut_off_column_top(nk, tv, GV, US, G_earth, depth, min_thickness, T, !! if associated real, intent(in) :: z_tol !< The tolerance with which to find the depth !! matching the specified pressure [Z ~> m]. - integer, optional, intent(in) :: remap_answer_date !< The vintage of the order of arithmetic and - !! expressions to use for remapping. Values below 20190101 - !! recover the remapping answers from 2018, while higher - !! values use more robust forms of the same remapping expressions. ! Local variables real, dimension(nk+1) :: e ! Top and bottom edge positions for reconstructions [Z ~> m] @@ -1332,11 +1333,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 +1405,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 +2751,14 @@ 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 + if (remap_general) then + dz_neglect = set_h_neglect(GV, remap_answer_date, dz_neglect_edge) + else + dz_neglect = set_dz_neglect(GV, US, remap_answer_date, dz_neglect_edge) + endif 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 @@ -2774,9 +2773,9 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just deallocate( dz_interface ) call ALE_remap_scalar(remapCS, G, GV, nkd, h1, tmpT1dIn, h, tv%T, all_cells=remap_full_column, & - old_remap=remap_old_alg, answer_date=remap_answer_date ) + old_remap=remap_old_alg ) call ALE_remap_scalar(remapCS, G, GV, nkd, h1, tmpS1dIn, h, tv%S, all_cells=remap_full_column, & - old_remap=remap_old_alg, answer_date=remap_answer_date ) + old_remap=remap_old_alg ) else ! This is the old way of initializing to z* coordinates only allocate( hTarget(nz) ) @@ -2799,11 +2798,9 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just 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) 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) if (GV%Boussinesq .or. GV%semi_Boussinesq) then ! This is a simple conversion of the target grid to thickness units that is not @@ -3106,9 +3103,17 @@ subroutine MOM_state_init_tests(G, GV, US, tv) write(0,*) ' ==================================================================== ' write(0,*) '' write(0,*) GV%H_to_m*h(:) + + ! For consistency with the usual call, add the following: + ! if (use_remapping) then + ! allocate(remap_CS) + ! call initialize_remapping(remap_CS, 'PLM', boundary_extrapolation=.true., & + ! h_neglect=GV%H_subroundoff, h_neglect_edge=GV%H_subroundoff) + ! endif call cut_off_column_top(nk, tv, GV, US, GV%g_Earth, -e(nk+1), GV%Angstrom_H, & T, T_t, T_b, S, S_t, S_b, 0.5*P_tot, h, remap_CS, z_tol=z_tol) write(0,*) GV%H_to_m*h(:) + if (associated(remap_CS)) deallocate(remap_CS) end subroutine MOM_state_init_tests diff --git a/src/initialization/MOM_tracer_initialization_from_Z.F90 b/src/initialization/MOM_tracer_initialization_from_Z.F90 index bafa5d8c36..6e3da385ce 100644 --- a/src/initialization/MOM_tracer_initialization_from_Z.F90 +++ b/src/initialization/MOM_tracer_initialization_from_Z.F90 @@ -13,7 +13,7 @@ module MOM_tracer_initialization_from_Z use MOM_grid, only : ocean_grid_type use MOM_horizontal_regridding, only : myStats, horiz_interp_and_extrap_tracer use MOM_interface_heights, only : dz_to_thickness_simple -use MOM_regridding, only : set_dz_neglect +use MOM_regridding, only : set_dz_neglect, set_h_neglect use MOM_remapping, only : remapping_CS, initialize_remapping use MOM_unit_scaling, only : unit_scale_type use MOM_verticalGrid, only : verticalGrid_type @@ -93,9 +93,9 @@ subroutine MOM_initialize_tracer_from_Z(h, tr, G, GV, US, PF, src_file, src_var_ real :: missing_value ! A value indicating that there is no valid input data at this point [CU ~> conc] real :: dz_neglect ! A negligibly small vertical layer extent used in - ! remapping cell reconstructions [Z ~> m] + ! remapping cell reconstructions [Z ~> m] or [H ~> m or kg m-2] real :: dz_neglect_edge ! A negligibly small vertical layer extent used in - ! remapping edge value calculations [Z ~> m] + ! remapping edge value calculations [Z ~> m] or [H ~> m or kg m-2] logical :: om4_remap_via_sub_cells ! If true, use the OM4 remapping algorithm integer :: nPoints ! The number of valid input data points in a column integer :: id_clock_routine, id_clock_ALE @@ -117,7 +117,7 @@ subroutine MOM_initialize_tracer_from_Z(h, tr, G, GV, US, PF, src_file, src_var_ is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed - call callTree_enter(trim(mdl)//"(), MOM_state_initialization.F90") + call callTree_enter(trim(mdl)//"(), MOM_tracer_initialization_from_Z.F90") call get_param(PF, mdl, "Z_INIT_HOMOGENIZE", homog, & "If True, then horizontally homogenize the interpolated "//& @@ -178,9 +178,15 @@ subroutine MOM_initialize_tracer_from_Z(h, tr, G, GV, US, PF, src_file, src_var_ allocate( h1(kd) ) allocate( dzSrc(isd:ied,jsd:jed,kd) ) allocate( hSrc(isd:ied,jsd:jed,kd) ) - ! Set parameters for reconstructions + ! Set parameters for reconstructions in the right units + if (h_is_in_Z_units) then + dz_neglect = set_dz_neglect(GV, US, remap_answer_date, dz_neglect_edge) + else + dz_neglect = set_h_neglect(GV, remap_answer_date, dz_neglect_edge) + endif 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 @@ -206,18 +212,15 @@ subroutine MOM_initialize_tracer_from_Z(h, tr, G, GV, US, PF, src_file, src_var_ enddo ; enddo 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) + ! Because h is in units of [Z ~> m], dzSrc is already in the right units. + call ALE_remap_scalar(remapCS, G, GV, kd, dzSrc, tr_z, h, tr, all_cells=.false.) 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. GV_loc = GV ; GV_loc%ke = kd call dz_to_thickness_simple(dzSrc, hSrc, G, GV_loc, US) - call ALE_remap_scalar(remapCS, G, GV, kd, hSrc, tr_z, h, tr, all_cells=.false., answer_date=remap_answer_date ) + call ALE_remap_scalar(remapCS, G, GV, kd, hSrc, tr_z, h, tr, all_cells=.false.) endif deallocate( hSrc ) diff --git a/src/ocean_data_assim/MOM_oda_driver.F90 b/src/ocean_data_assim/MOM_oda_driver.F90 index 9275555afc..9d7c5cf05a 100644 --- a/src/ocean_data_assim/MOM_oda_driver.F90 +++ b/src/ocean_data_assim/MOM_oda_driver.F90 @@ -55,7 +55,7 @@ module MOM_oda_driver_mod use MOM_domains, only : MOM_domains_init, MOM_domain_type, clone_MOM_domain use MOM_remapping, only : remapping_CS, initialize_remapping, remapping_core_h use MOM_regridding, only : regridding_CS, initialize_regridding -use MOM_regridding, only : regridding_main, set_regrid_params +use MOM_regridding, only : regridding_main, set_regrid_params, set_h_neglect use MOM_unit_scaling, only : unit_scale_type, unit_scaling_init use MOM_variables, only : thermo_var_ptrs use MOM_verticalGrid, only : verticalGrid_type, verticalGridInit @@ -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,8 +324,12 @@ 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.) + 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) + + h_neglect = set_h_neglect(GV, CS%answer_date, h_neglect_edge) + call initialize_remapping(CS%remapCS, remap_scheme, om4_remap_via_sub_cells=om4_remap_via_sub_cells, & + h_neglect=h_neglect, h_neglect_edge=h_neglect_edge) call set_regrid_params(CS%regridCS, min_thickness=0.) isd = G%isd; ied = G%ied; jsd = G%jsd; jed = G%jed @@ -415,7 +420,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 +431,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 +439,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 @@ -683,7 +679,6 @@ subroutine apply_oda_tracer_increments(dt, Time_end, G, GV, tv, h, CS) !! DA [C T-1 ~> degC s-1] real, dimension(SZI_(G),SZJ_(G),SZK_(CS%Grid)) :: S_tend !< The salinity tendency adjustment from DA !! [S T-1 ~> ppt s-1] - real :: h_neglect, h_neglect_edge ! small thicknesses [H ~> m or kg m-2] if (.not. associated(CS)) return if (CS%assim_method == NO_ASSIM .and. (.not. CS%do_bias_adjustment)) return @@ -700,20 +695,12 @@ subroutine apply_oda_tracer_increments(dt, Time_end, G, GV, tv, h, CS) S_tend = S_tend + CS%S_bc_tend endif - 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 - 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..ad3423cde5 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) @@ -554,7 +550,6 @@ subroutine apply_oda_incupd(h, tv, u, v, dt, G, GV, US, CS) integer :: isB, ieB, jsB, jeB ! integer :: ncount ! time step counter real :: inc_wt ! weight of the update for this time-step [nondim] - 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] character(len=256) :: mesg @@ -578,12 +573,6 @@ subroutine apply_oda_incupd(h, tv, u, v, dt, G, GV, US, CS) write(mesg,'(f10.8)') inc_wt if (is_root_pe()) call MOM_error(NOTE,"updating fields with weight inc_wt:"//trim(mesg)) - 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) @@ -621,7 +610,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) @@ -633,8 +622,8 @@ subroutine apply_oda_incupd(h, tv, u, v, dt, G, GV, US, CS) tmp_val2(k) = CS%Inc(2)%p(i,j,k) 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) + call remapping_core_h(CS%remap_cs, nz_data, tmp_h(1:nz_data), tmp_val2, & + 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 +669,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 +707,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..5b57103078 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..6d8fe881d1 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) + call extract_member_remapping_CS(CS%remap_CS, degree=CS%deg) 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.