From 0f2a69d52558d5a44486acdbc83272bd9b451b0c Mon Sep 17 00:00:00 2001 From: Kate Hedstrom Date: Wed, 8 Nov 2023 09:03:09 -0900 Subject: [PATCH] Obc tracer fix (#507) * +Fix for issue #506, tracer OBC bug - it only happens in the advection for certain flow directions, inflow from OBC plus along-boundary flow. * Tracer OBCs need more of an h halo update. - This one should finally fix the processor count issues with OBCs. * Correct the "if" statement. * +Adding more halo points to an exchange - This will change answers if you start with a non-zero velocity. You need three halo points (or maybe cont_stencil) for the continuity solver. - Also trying to put in some initial DEBUG_OBC code. * Fix some DEBUG_OBC logic * Writing to temporary arrays for tres_xy - Way to trick some compiler. * Another fix for DEBUG_OBC * Fixing the whalo troubles for grids that are 2 wide/long. * Exchange all the h_new points. - without this, because we're remapping all the tres points, it dies in debug mode on bad h_new values. * Trying a different exchange - as it was, it was passing my tests, failing the auto-tests. * Fixing the DEBUG_OBC logging * Putting the logging statement back. - Making an error more verbose too. --- src/core/MOM.F90 | 12 ++-- src/core/MOM_open_boundary.F90 | 57 ++++++++++++++----- src/framework/MOM_restart.F90 | 4 +- .../lateral/MOM_hor_visc.F90 | 5 +- 4 files changed, 54 insertions(+), 24 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index e1fa004fb8..ae794e02e2 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -1627,9 +1627,11 @@ subroutine step_MOM_thermo(CS, G, GV, US, u, v, h, tv, fluxes, dtdia, & if (CS%split) & call remap_dyn_split_RK2_aux_vars(G, GV, CS%dyn_split_RK2_CSp, h, h_new, CS%ALE_CSp, CS%OBC, dzRegrid) - if (associated(CS%OBC)) & - call pass_var(h_new, G%Domain) + if (associated(CS%OBC)) then + call pass_var(h, G%Domain, complete=.false.) + call pass_var(h_new, G%Domain, complete=.true.) call remap_OBC_fields(G, GV, h, h_new, CS%OBC, PCM_cell=PCM_cell) + endif call remap_vertvisc_aux_vars(G, GV, CS%visc, h, h_new, CS%ALE_CSp, CS%OBC) if (associated(CS%visc%Kv_shear)) & @@ -3016,10 +3018,10 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, & call cpu_clock_begin(id_clock_pass_init) call create_group_pass(tmp_pass_uv_T_S_h, CS%u, CS%v, G%Domain) if (use_temperature) then - call create_group_pass(tmp_pass_uv_T_S_h, CS%tv%T, G%Domain, halo=1) - call create_group_pass(tmp_pass_uv_T_S_h, CS%tv%S, G%Domain, halo=1) + call create_group_pass(tmp_pass_uv_T_S_h, CS%tv%T, G%Domain) + call create_group_pass(tmp_pass_uv_T_S_h, CS%tv%S, G%Domain) endif - call create_group_pass(tmp_pass_uv_T_S_h, CS%h, G%Domain, halo=1) + call create_group_pass(tmp_pass_uv_T_S_h, CS%h, G%Domain) call do_group_pass(tmp_pass_uv_T_S_h, G%Domain) call cpu_clock_end(id_clock_pass_init) diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index 13ce524006..5cf7c92fe9 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -7,33 +7,33 @@ module MOM_open_boundary use MOM_array_transform, only : allocate_rotated_array use MOM_coms, only : sum_across_PEs, Set_PElist, Get_PElist, PE_here, num_PEs use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, CLOCK_ROUTINE +use MOM_debugging, only : hchksum, uvchksum use MOM_diag_mediator, only : diag_ctrl, time_type use MOM_domains, only : pass_var, pass_vector use MOM_domains, only : create_group_pass, do_group_pass, group_pass_type use MOM_domains, only : To_All, EAST_FACE, NORTH_FACE, SCALAR_PAIR, CGRID_NE, CORNER +use MOM_dyn_horgrid, only : dyn_horgrid_type use MOM_error_handler, only : MOM_mesg, MOM_error, FATAL, WARNING, NOTE, is_root_pe use MOM_file_parser, only : get_param, log_version, param_file_type, log_param use MOM_grid, only : ocean_grid_type, hor_index_type -use MOM_dyn_horgrid, only : dyn_horgrid_type use MOM_interface_heights, only : thickness_to_dz +use MOM_interpolate, only : init_external_field, time_interp_external, time_interp_external_init +use MOM_interpolate, only : external_field use MOM_io, only : slasher, field_size, SINGLE_FILE use MOM_io, only : vardesc, query_vardesc, var_desc +use MOM_obsolete_params, only : obsolete_logical, obsolete_int, obsolete_real, obsolete_char +use MOM_regridding, only : regridding_CS +use MOM_remapping, only : remappingSchemesDoc, remappingDefaultScheme, remapping_CS +use MOM_remapping, only : initialize_remapping, remapping_core_h, end_remapping use MOM_restart, only : register_restart_field, register_restart_pair use MOM_restart, only : query_initialized, MOM_restart_CS -use MOM_obsolete_params, only : obsolete_logical, obsolete_int, obsolete_real, obsolete_char -use MOM_string_functions, only : extract_word, remove_spaces, uppercase +use MOM_string_functions, only : extract_word, remove_spaces, uppercase, lowercase use MOM_tidal_forcing, only : astro_longitudes, astro_longitudes_init, eq_phase, nodal_fu, tidal_frequency use MOM_time_manager, only : set_date, time_type, time_type_to_real, operator(-) use MOM_tracer_registry, only : tracer_type, tracer_registry_type, tracer_name_lookup -use MOM_interpolate, only : init_external_field, time_interp_external, time_interp_external_init -use MOM_interpolate, only : external_field -use MOM_remapping, only : remappingSchemesDoc, remappingDefaultScheme, remapping_CS -use MOM_remapping, only : initialize_remapping, remapping_core_h, end_remapping -use MOM_regridding, only : regridding_CS use MOM_unit_scaling, only : unit_scale_type use MOM_variables, only : thermo_var_ptrs use MOM_verticalGrid, only : verticalGrid_type -use MOM_string_functions, only : lowercase implicit none ; private @@ -528,13 +528,16 @@ subroutine open_boundary_config(G, US, param_file, OBC) OBC%add_tide_constituents = .false. endif - call get_param(param_file, mdl, "DEBUG", OBC%debug, default=.false.) - call get_param(param_file, mdl, "DEBUG_OBC", debug_OBC, default=.false.) - if (debug_OBC .or. OBC%debug) & + call get_param(param_file, mdl, "DEBUG", debug_OBC, default=.false.) + call get_param(param_file, mdl, "DEBUG_OBC", debug_OBC, default=debug_OBC, & + do_not_log=.not.debug_OBC) + if (debug_OBC) then call log_param(param_file, mdl, "DEBUG_OBC", debug_OBC, & "If true, do additional calls to help debug the performance "//& "of the open boundary condition code.", default=.false., & debuggingParam=.true.) + OBC%debug = debug_OBC + endif call get_param(param_file, mdl, "OBC_SILLY_THICK", OBC%silly_h, & "A silly value of thicknesses used outside of open boundary "//& @@ -854,6 +857,8 @@ subroutine initialize_segment_data(G, GV, US, OBC, PF) ! if (siz(4) == 1) segment%values_needed = .false. if (segment%on_pe) then if (OBC%brushcutter_mode .and. (modulo(siz(1),2) == 0 .or. modulo(siz(2),2) == 0)) then + write(mesg,'("Brushcutter mode sizes ", I6, I6))') siz(1), siz(2) + call MOM_error(WARNING, mesg // " " // trim(filename) // " " // trim(fieldname)) call MOM_error(FATAL,'segment data are not on the supergrid') endif siz2(1)=1 @@ -2224,6 +2229,8 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, GV, US, type(OBC_segment_type), pointer :: segment => NULL() integer :: i, j, k, is, ie, js, je, m, nz, n integer :: is_obc, ie_obc, js_obc, je_obc + logical :: sym + character(len=3) :: var_num is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke @@ -3298,6 +3305,29 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, GV, US, call pass_vector(u_new, v_new, G%Domain, clock=id_clock_pass) + if (OBC%debug) then + sym = G%Domain%symmetric + if (OBC%radiation_BCs_exist_globally) then + call uvchksum("radiation_OBCs: OBC%r[xy]_normal", OBC%rx_normal, OBC%ry_normal, G%HI, & + haloshift=0, symmetric=sym, scale=1.0) + endif + if (OBC%oblique_BCs_exist_globally) then + call uvchksum("radiation_OBCs: OBC%r[xy]_oblique_[uv]", OBC%rx_oblique_u, OBC%ry_oblique_v, G%HI, & + haloshift=0, symmetric=sym, scale=1.0/US%L_T_to_m_s**2) + call uvchksum("radiation_OBCs: OBC%r[yx]_oblique_[uv]", OBC%ry_oblique_u, OBC%rx_oblique_v, G%HI, & + haloshift=0, symmetric=sym, scale=1.0/US%L_T_to_m_s**2) + call uvchksum("radiation_OBCs: OBC%cff_normal_[uv]", OBC%cff_normal_u, OBC%cff_normal_v, G%HI, & + haloshift=0, symmetric=sym, scale=1.0/US%L_T_to_m_s**2) + endif + if (OBC%ntr == 0) return + if (.not. allocated (OBC%tres_x) .or. .not. allocated (OBC%tres_y)) return + do m=1,OBC%ntr + write(var_num,'(I3.3)') m + call uvchksum("radiation_OBCs: OBC%tres_[xy]_"//var_num, OBC%tres_x(:,:,:,m), OBC%tres_y(:,:,:,m), G%HI, & + haloshift=0, symmetric=sym, scale=1.0) + enddo + endif + end subroutine radiation_open_bdry_conds !> Applies OBC values stored in segments to 3d u,v fields @@ -5638,9 +5668,6 @@ subroutine remap_OBC_fields(G, GV, h_old, h_new, OBC, PCM_cell) To_All+Scalar_Pair) if (OBC%oblique_BCs_exist_globally) then call do_group_pass(OBC%pass_oblique, G%Domain) -! call pass_vector(OBC%rx_oblique_u, OBC%ry_oblique_v, G%Domain, To_All+Scalar_Pair) -! call pass_vector(OBC%ry_oblique_u, OBC%rx_oblique_v, G%Domain, To_All+Scalar_Pair) -! call pass_vector(OBC%cff_normal_u, OBC%cff_normal_v, G%Domain, To_All+Scalar_Pair) endif end subroutine remap_OBC_fields diff --git a/src/framework/MOM_restart.F90 b/src/framework/MOM_restart.F90 index 252f14bfac..188cfbb2ec 100644 --- a/src/framework/MOM_restart.F90 +++ b/src/framework/MOM_restart.F90 @@ -660,7 +660,7 @@ end subroutine register_restart_field_0d !> query_initialized_name determines whether a named field has been successfully -!! read from a restart file or has otherwise been recored as being initialzed. +!! read from a restart file or has otherwise been recorded as being initialized. function query_initialized_name(name, CS) result(query_initialized) character(len=*), intent(in) :: name !< The name of the field that is being queried type(MOM_restart_CS), intent(in) :: CS !< MOM restart control struct @@ -1271,7 +1271,7 @@ subroutine only_read_restart_pair_3d(a_ptr, b_ptr, a_name, b_name, G, CS, & end subroutine only_read_restart_pair_3d -!> Return an indicationof whether the named variable is the restart files, and provie the full path +!> Return an indication of whether the named variable is in the restart files, and provide the full path !! to the restart file in which a variable is found. function find_var_in_restart_files(varname, G, CS, file_path, filename, directory, is_global) result (found) character(len=*), intent(in) :: varname !< The variable name to be used in the restart file diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index 02b4ec66a6..e3249afb73 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -414,12 +414,13 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB h_neglect = GV%H_subroundoff - h_neglect3 = h_neglect**3 + !h_neglect3 = h_neglect**3 + h_neglect3 = h_neglect*h_neglect*h_neglect inv_PI3 = 1.0/((4.0*atan(1.0))**3) inv_PI2 = 1.0/((4.0*atan(1.0))**2) inv_PI6 = inv_PI3 * inv_PI3 - m_leithy(:,:) = 0. ! Initialize + m_leithy(:,:) = 0.0 ! Initialize if (present(OBC)) then ; if (associated(OBC)) then ; if (OBC%OBC_pe) then apply_OBC = OBC%Flather_u_BCs_exist_globally .or. OBC%Flather_v_BCs_exist_globally