From 57ea17b16c58151f4d78c8a2df0d7b25c4692351 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Fri, 29 Aug 2025 16:37:07 -0400 Subject: [PATCH 1/6] +Restore open_boundary_setup_vert Returned the functionality of open_boundary_setup_vert() to a separate routine, and no longer have initialize_segment_data() carry out these initialization steps. These routines had previously been combined in the expectation that they would always occur at about the same point in the code but this turns out not to be the case for generic tracers with open boundary conditions. This renewed routine is called directly form initialize_MOM. All answers are bitwise identical but there is one added public OBC interface. --- src/core/MOM.F90 | 7 ++-- src/core/MOM_open_boundary.F90 | 66 ++++++++++++++++++++-------------- 2 files changed, 44 insertions(+), 29 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 984ba8f620..26c533f21d 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -116,7 +116,7 @@ module MOM use MOM_open_boundary, only : register_temp_salt_segments, update_segment_tracer_reservoirs use MOM_open_boundary, only : setup_OBC_tracer_reservoirs use MOM_open_boundary, only : open_boundary_register_restarts, remap_OBC_fields -use MOM_open_boundary, only : initialize_segment_data, rotate_OBC_config +use MOM_open_boundary, only : open_boundary_setup_vert, initialize_segment_data, rotate_OBC_config use MOM_open_boundary, only : update_OBC_segment_data, open_boundary_halo_update use MOM_open_boundary, only : write_OBC_info, chksum_OBC_segments use MOM_porous_barriers, only : porous_widths_layer, porous_widths_interface, porous_barriers_init @@ -2873,9 +2873,10 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, & CS%HFrz = (US%Z_to_m * GV%m_to_H) * HFrz_z if (associated(OBC_in)) then - ! This call allocates the arrays on the segments for open boundary data and initializes the - ! relevant vertical remapping structures. It can only occur after the vertical grid has been + ! These calls allocate the arrays on the segments for open boundary data and initialize the + ! relevant vertical remapping structures. They can only occur after the vertical grid has been ! initialized. + call open_boundary_setup_vert(GV, US, OBC_in) call initialize_segment_data(G_in, GV, US, OBC_in, param_file) endif diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index 68366e91f4..7df8f11c6e 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -40,6 +40,7 @@ module MOM_open_boundary public open_boundary_apply_normal_flow public open_boundary_config +public open_boundary_setup_vert public open_boundary_halo_update public open_boundary_query public open_boundary_end @@ -801,9 +802,46 @@ 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 + ! Set up vertical remapping for open boundaries. Remapping happens independently on each PE, + ! so this block could be skipped for PEs without open boundary conditions that use remapping. + 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 -!> Set up vertical remapping and 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. +end subroutine open_boundary_setup_vert + + +!> Get and store properties about the fields on the OBC segments and allocate space for reading +!! OBC data from files. In the process, it does funky stuff with the MPI processes. subroutine initialize_segment_data(G, GV, US, OBC, PF) type(ocean_grid_type), intent(in) :: G !< Ocean grid structure type(verticalGrid_type), intent(in) :: GV !< Container for vertical grid information @@ -835,30 +873,6 @@ subroutine initialize_segment_data(G, GV, US, OBC, PF) is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec - if (OBC%number_of_segments > 0) then - ! Set up vertical remapping for open boundaries. Remapping happens independently on each PE, - ! so this block could be skipped for PEs without open boundary conditions that use remapping. - 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 - ! There is a problem with the order of the OBC initialization ! with respect to ALE_init. Currently handling this by copying the ! param file so that I can use it later in step_MOM in order to finish From 70b37225bff8eb4d88b1a2723a5510c1e26cd6a8 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Fri, 22 Aug 2025 13:38:17 -0400 Subject: [PATCH 2/6] +Separate OBC_for_remap and OBC_for_bug arguments Add separate OBC_for_remap and OBC_for_bug arguments to MOM_initialize_state() to clarify the distinct roles of the two arguments, to make it clear that one can not be used without rotation, and to prepare for the eventual obsoleting of OBC_RESERVOIR_INIT_BUG. This commit also simplifies the logic setting ntstep when do_thermo is false. All answers are bitwise identical, but there is a new optional argument to MOM_initialize_state and another argument has been renamed. --- src/core/MOM.F90 | 15 ++++++------ .../MOM_state_initialization.F90 | 24 ++++++++++++------- 2 files changed, 23 insertions(+), 16 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 26c533f21d..001be33812 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -673,7 +673,7 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS if (do_dyn) then n_max = 1 if (time_interval > CS%dt) n_max = ceiling(time_interval/CS%dt - 0.001) - ntstep = 1 ! initialization + dt = time_interval / real(n_max) thermo_does_span_coupling = (CS%thermo_spans_coupling .and. & (CS%dt_therm > 1.5*cycle_time)) @@ -686,7 +686,8 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS elseif (.not.do_thermo) then dt_therm = CS%dt_therm if (present(cycle_length)) dt_therm = min(CS%dt_therm, cycle_length) - ! ntstep is not used. + ntstep = 1 ! ntstep is initialized to avoid an error in a secondary logical test, + ! but the nonzero value of ntstep does not matter when do_thermo is false. else ntstep = MAX(1, MIN(n_max, floor(CS%dt_therm/dt + 0.001))) dt_therm = dt*ntstep @@ -698,7 +699,7 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS elseif (.not.do_thermo) then dt_tr_adv = CS%dt_tr_adv if (present(cycle_length)) dt_tr_adv = min(CS%dt_tr_adv, cycle_length) - ! ntstep is not used. + ! ntastep is not used. else ntastep = MAX(1, MIN(n_max, floor(CS%dt_tr_adv/dt + 0.001))) dt_tr_adv = dt*ntastep @@ -3179,7 +3180,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, & CS%tv%T => T_in CS%tv%S => S_in - if (associated(OBC_in)) then + if (associated(CS%OBC)) then ! Log this parameter in MOM_initialize_state call get_param(param_file, "MOM", "OBC_RESERVOIR_INIT_BUG", OBC_reservoir_init_bug, & "If true, set the OBC tracer reservoirs at the startup of a new run from the "//& @@ -3260,12 +3261,12 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, & call ice_shelf_query(ice_shelf_CSp,G,CS%frac_shelf_h, CS%mass_shelf) call MOM_initialize_state(CS%u, CS%v, CS%h, CS%tv, Time, G, GV, US, & param_file, dirs, restart_CSp, CS%ALE_CSp, CS%tracer_Reg, & - CS%sponge_CSp, CS%ALE_sponge_CSp,CS%oda_incupd_CSp, CS%OBC, Time_in, & - frac_shelf_h=CS%frac_shelf_h, mass_shelf=CS%mass_shelf) + CS%sponge_CSp, CS%ALE_sponge_CSp, CS%oda_incupd_CSp, CS%OBC, Time_in, & + frac_shelf_h=CS%frac_shelf_h, mass_shelf=CS%mass_shelf, OBC_for_bug=CS%OBC) else call MOM_initialize_state(CS%u, CS%v, CS%h, CS%tv, Time, G, GV, US, & param_file, dirs, restart_CSp, CS%ALE_CSp, CS%tracer_Reg, & - CS%sponge_CSp, CS%ALE_sponge_CSp, CS%oda_incupd_CSp, CS%OBC, Time_in) + CS%sponge_CSp, CS%ALE_sponge_CSp, CS%oda_incupd_CSp, CS%OBC, Time_in, OBC_for_bug=CS%OBC) endif ! Reset the first direction if it was found in a restart file. diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index 2238f95e74..bcefb0022e 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -113,7 +113,8 @@ module MOM_state_initialization !! conditions or by reading them from a restart (or saves) file. subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & restart_CS, ALE_CSp, tracer_Reg, sponge_CSp, & - ALE_sponge_CSp, oda_incupd_CSp, OBC, Time_in, frac_shelf_h, mass_shelf) + ALE_sponge_CSp, oda_incupd_CSp, OBC_for_remap, & + Time_in, frac_shelf_h, mass_shelf, OBC_for_bug) type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type @@ -137,8 +138,10 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & type(tracer_registry_type), pointer :: tracer_Reg !< A pointer to the tracer registry type(sponge_CS), pointer :: sponge_CSp !< The layerwise sponge control structure. type(ALE_sponge_CS), pointer :: ALE_sponge_CSp !< The ALE sponge control structure. - type(ocean_OBC_type), pointer :: OBC !< The open boundary condition control structure. - ! OBC is only used in MOM_initialize_state if OBC_RESERVOIR_INIT_BUG is true. + type(ocean_OBC_type), pointer :: OBC_for_remap !< The open boundary condition control + !! structure that may be used for remapping velocities. + !! This must be on the unrotated grid, but only the + !! position and directions of the OBC faces are used. type(oda_incupd_CS), pointer :: oda_incupd_CSp !< The oda_incupd control structure. type(time_type), optional, intent(in) :: Time_in !< Time at the start of the run segment. real, dimension(SZI_(G),SZJ_(G)), & @@ -147,6 +150,9 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & real, dimension(SZI_(G),SZJ_(G)), & optional, intent(in) :: mass_shelf !< The mass per unit area of the overlying !! ice shelf [ R Z ~> kg m-2 ] + type(ocean_OBC_type), optional, pointer :: OBC_for_bug !< An open boundary condition control structure + !! that might be used to store OBC temperatures and + !! salinities if OBC_RESERVOIR_INIT_BUG is true. ! Local variables real :: depth_tot(SZI_(G),SZJ_(G)) ! The nominal total depth of the ocean [Z ~> m] real :: dz(SZI_(G),SZJ_(G),SZK_(GV)) ! The layer thicknesses in geopotential (z) units [Z ~> m] @@ -432,7 +438,7 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & endif endif ! not from_Z_file. - if (use_temperature .and. associated(OBC)) then + if (present(OBC_for_bug)) then ; if (use_temperature .and. associated(OBC_for_bug)) then call get_param(PF, mdl, "ENABLE_BUGS_BY_DEFAULT", enable_bugs, & default=.true., do_not_log=.true.) ! This is logged from MOM.F90. ! Log this parameter later with the other OBC parameters. @@ -445,9 +451,9 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & ! the temperatures and salinities can change due to the remapping and reading from the restarts. call pass_var(tv%T, G%Domain, complete=.false.) call pass_var(tv%S, G%Domain, complete=.true.) - call fill_temp_salt_segments(G, GV, US, OBC, tv) + call fill_temp_salt_segments(G, GV, US, OBC_for_bug, tv) endif - endif + endif ; endif ! Convert thicknesses from geometric distances in depth units to thickness units or mass-per-unit-area. if (new_sim .and. convert) call dz_to_thickness(dz, tv, h, G, GV, US) @@ -496,10 +502,10 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & if (new_sim .and. debug) & call hchksum(h, "Pre-ALE_regrid: h ", G%HI, haloshift=1, unscale=GV%H_to_MKS) - ! In this call, OBC is only used for the directions of OBCs when setting thicknesses at + ! In this call, OBC_for_remap is only used for the directions of OBCs when setting thicknesses at ! velocity points. - call ALE_regrid_accelerated(ALE_CSp, G, GV, US, h, tv, regrid_iterations, u, v, OBC, tracer_Reg, & - dt=dt, initial=.true.) + call ALE_regrid_accelerated(ALE_CSp, G, GV, US, h, tv, regrid_iterations, u, v, OBC_for_remap, & + tracer_Reg, dt=dt, initial=.true.) endif endif From 3f8de40f6cdba791a6f627cf426011ca0bca289f Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Sun, 24 Aug 2025 09:20:18 -0400 Subject: [PATCH 3/6] +Add chksum_OBC_segment_data Added the new subroutine chksum_OBC_segment_data() from the contents of the loop over the segments in chksum_OBC_segments, as this can be useful when called separately for debugging some of the constructs that will be coming in subsequent commits. Also added rotate_OBC_segment_values_needed() to copy over the values_needed and _index variables between segments, and to avoid the need for duplicate code blocks in rotate_OBC_segment_data() and rotate_OBC_config(). There is also some minor revisions to rename segnam to segname for greater clarity. This commit also makes greater use of the I0 format that was introduced with Fortran 95 to simplify or shorten some error messages. All answers are bitwise identical but there are two new internal subroutines in the MOM_open_boundary module. --- src/core/MOM_open_boundary.F90 | 219 +++++++++++++++++---------------- 1 file changed, 113 insertions(+), 106 deletions(-) diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index 7df8f11c6e..c3f02a0e9b 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -751,18 +751,18 @@ subroutine open_boundary_config(G, US, param_file, OBC) 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.) + "detected then a FATAL error is issued.", default=.false., do_not_log=.true.) 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.) + "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", 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.) + "round off.", default=.false., do_not_log=.true.) call get_param(param_file, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, & "This sets the default value for the various _ANSWER_DATE parameters.", & default=99991231) @@ -852,7 +852,7 @@ subroutine initialize_segment_data(G, GV, US, OBC, PF) integer :: n, n_seg, m, num_fields, mm character(len=1024) :: segstr character(len=256) :: filename - character(len=20) :: segnam, suffix + character(len=20) :: segname, suffix character(len=32) :: fieldname real :: value ! A value that is parsed from the segment data string [various units] real :: dz_neglect, dz_neglect_edge ! Small thicknesses in vertical height units [Z ~> m] @@ -873,6 +873,7 @@ subroutine initialize_segment_data(G, GV, US, OBC, PF) is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec + ! ###IS THIS STILL TRUE? ! There is a problem with the order of the OBC initialization ! with respect to ALE_init. Currently handling this by copying the ! param file so that I can use it later in step_MOM in order to finish @@ -885,8 +886,8 @@ subroutine initialize_segment_data(G, GV, US, OBC, PF) ! Try this here just for the documentation. It is repeated below. do n=1,OBC%number_of_segments - write(segnam,"('OBC_SEGMENT_',i3.3,'_DATA')") n - call get_param(PF, mdl, segnam, segstr, 'OBC segment docs') + write(segname, "('OBC_SEGMENT_',i3.3,'_DATA')") n + call get_param(PF, mdl, segname, segstr, 'OBC segment docs') enddo !< temporarily disable communication in order to read segment data independently @@ -901,18 +902,18 @@ subroutine initialize_segment_data(G, GV, US, OBC, PF) n_seg = n ; if (OBC%reverse_segment_order) n_seg = OBC%number_of_segments + 1 - n segment => OBC%segment(n_seg) ! segment%values_needed is only true if this segment is on the local PE and some values need to be read. - if (.not. segment%values_needed) cycle + if (.not. OBC%segment(n_seg)%values_needed) cycle - write(segnam,"('OBC_SEGMENT_',i3.3,'_DATA')") n - write(suffix,"('_segment_',i3.3)") n + write(segname, "('OBC_SEGMENT_',i3.3,'_DATA')") n + write(suffix, "('_segment_',i3.3)") n ! needs documentation !! Yet, unsafe for now, causes grief for ! MOM_parameter_docs in circle_obcs on two processes. -! call get_param(PF, mdl, segnam, segstr, 'xyz') +! call get_param(PF, mdl, segname, segstr, 'xyz') ! Clear out any old values segstr = '' - call get_param(PF, mdl, segnam, segstr) + call get_param(PF, mdl, segname, segstr) if (segstr == '') then - write(mesg,'("No OBC_SEGMENT_XXX_DATA string for OBC segment ",I3)') n + write(mesg,'("No OBC_SEGMENT_XXX_DATA string for OBC segment ",I0)') n call MOM_error(FATAL, mesg) endif @@ -983,7 +984,7 @@ subroutine initialize_segment_data(G, GV, US, OBC, PF) call MOM_error(FATAL," Unable to open OBC file " // trim(filename)) 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) + write(mesg,'("Brushcutter mode sizes ", I0, I0)') 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 @@ -1074,7 +1075,7 @@ subroutine initialize_segment_data(G, GV, US, OBC, PF) segment%v_values_needed .or. segment%vamp_values_needed .or. segment%vphase_values_needed .or. & segment%t_values_needed .or. segment%s_values_needed .or. segment%g_values_needed .or. & segment%z_values_needed .or. segment%zamp_values_needed .or. segment%zphase_values_needed ) then - write(mesg,'("Values needed for OBC segment ",I3)') n + write(mesg,'("Values needed for OBC segment ",I0)') n call MOM_error(FATAL, mesg) endif enddo @@ -1812,10 +1813,10 @@ subroutine parse_segment_str(ni_global, nj_global, segment_str, l, m, n, action_ ! checking if the number of provided OBC types is less than or equal to 8 if (extract_word(segment_str,',',3+size(action_str))/="") then - write(max_words, '(I3)') size(action_str) + write(max_words, '(I0)') size(action_str) call MOM_error(FATAL, "MOM_open_boundary.F90, parse_segment_str: "// & "Number of OBC descriptor words in '" // trim(segment_str) // "' is too large. " // & - "There can be at most " // trim(adjustl(max_words)) // " descriptor words.") + "There can be at most " // trim(max_words) // " descriptor words.") endif ! Type of open boundary condition @@ -1938,7 +1939,7 @@ subroutine parse_for_tracer_reservoirs(OBC, PF, use_temperature) integer :: na character(len=1024) :: segstr character(len=256) :: filename - character(len=20) :: segnam, suffix + character(len=20) :: segname, suffix character(len=32) :: fieldname real :: value ! A value that is parsed from the segment data string [various units] character(len=32), dimension(MAX_OBC_FIELDS) :: fields ! segment field names @@ -1947,11 +1948,11 @@ subroutine parse_for_tracer_reservoirs(OBC, PF, use_temperature) do n=1,OBC%number_of_segments n_seg = n ; if (OBC%reverse_segment_order) n_seg = OBC%number_of_segments + 1 - n segment => OBC%segment(n_seg) - write(segnam,"('OBC_SEGMENT_',i3.3,'_DATA')") n - write(suffix,"('_segment_',i3.3)") n + write(segname, "('OBC_SEGMENT_',i3.3,'_DATA')") n + write(suffix, "('_segment_',i3.3)") n ! Clear out any old values segstr = '' - call get_param(PF, mdl, segnam, segstr) + call get_param(PF, mdl, segname, segstr) if (segstr == '') cycle call parse_segment_manifest_str(trim(segstr), num_fields, fields) @@ -4610,8 +4611,7 @@ subroutine update_OBC_ramp(Time, OBC, US, activate) OBC%ramp_value = wghtA endif write(msg(1:12),'(es12.3)') OBC%ramp_value - call MOM_error(NOTE, "MOM_open_boundary: update_OBC_ramp set OBC"// & - " ramp to "//trim(msg)) + call MOM_error(NOTE, "MOM_open_boundary: update_OBC_ramp set OBC ramp to "//trim(msg)) end subroutine update_OBC_ramp !> register open boundary objects for boundary updates. @@ -4625,7 +4625,7 @@ subroutine register_OBC(name, param_file, Reg) if (.not. associated(Reg)) call OBC_registry_init(param_file, Reg) if (Reg%nobc>=MAX_FIELDS_) then - write(mesg,'("Increase MAX_FIELDS_ in MOM_memory.h to at least ",I3," to allow for & + write(mesg, '("Increase MAX_FIELDS_ in MOM_memory.h to at least ",I0," to allow for & &all the open boundaries being registered via register_OBC.")') Reg%nobc+1 call MOM_error(FATAL,"MOM register_OBC: "//mesg) endif @@ -4658,9 +4658,8 @@ subroutine OBC_registry_init(param_file, Reg) init_calls = init_calls + 1 if (init_calls > 1) then - write(mesg,'("OBC_registry_init called ",I3, & - &" times with different registry pointers.")') init_calls - if (is_root_pe()) call MOM_error(WARNING,"MOM_open_boundary"//mesg) + write(mesg,'("OBC_registry_init called ",I0," times with different registry pointers.")') init_calls + if (is_root_pe()) call MOM_error(WARNING,"MOM_open_boundary: "//trim(mesg)) endif end subroutine OBC_registry_init @@ -4719,13 +4718,6 @@ subroutine segment_tracer_registry_init(param_file, segment) ! Read all relevant parameters and write them to the model log. if (init_calls == 1) call log_version(param_file, mdl, version, "") -! Need to call once per segment with tracers... -! if (init_calls > 1) then -! write(mesg,'("segment_tracer_registry_init called ",I3, & -! &" times with different registry pointers.")') init_calls -! if (is_root_pe()) call MOM_error(WARNING,"MOM_tracer"//mesg) -! endif - end subroutine segment_tracer_registry_init !> Register a tracer array that is active on an OBC segment, potentially also specifying how the @@ -4764,7 +4756,7 @@ subroutine register_segment_tracer(tr_ptr, ntr_index, param_file, GV, segment, & call segment_tracer_registry_init(param_file, segment) if (segment%tr_Reg%ntseg>=MAX_FIELDS_) then - write(mesg,'("Increase MAX_FIELDS_ in MOM_memory.h to at least ",I3," to allow for & + write(mesg,'("Increase MAX_FIELDS_ in MOM_memory.h to at least ",I0," to allow for & &all the tracers being registered via register_segment_tracer.")') segment%tr_Reg%ntseg+1 call MOM_error(FATAL,"MOM register_segment_tracer: "//mesg) endif @@ -5183,7 +5175,7 @@ subroutine mask_outside_OBCs(G, US, param_file, OBC) do j=G%jsd,G%jed ; do i=G%isd,G%ied if (color(i,j) /= color2(i,j)) then fatal_error = .True. - write(mesg,'("MOM_open_boundary: problem with OBC segments specification at ",I5,",",I5," during\n", & + write(mesg,'("MOM_open_boundary: problem with OBC segments specification at ",I0,",",I0," during\n", & &"the masking of the outside grid points.")') i, j call MOM_error(WARNING,"MOM mask_outside_OBCs: "//mesg, all_print=.true.) endif @@ -5391,7 +5383,7 @@ subroutine open_boundary_register_restarts(HI, GV, US, OBC, Reg, param_file, res ! This would be coming from user code such as DOME. if (OBC%ntr /= Reg%ntr) then ! call MOM_error(FATAL, "open_boundary_register_restarts: Inconsistent value for ntr") - write(mesg,'("Inconsistent values for ntr ", I8," and ",I8,".")') OBC%ntr, Reg%ntr + write(mesg,'("Inconsistent values for ntr ", I0," and ",I0,".")') OBC%ntr, Reg%ntr call MOM_error(WARNING, 'open_boundary_register_restarts: '//mesg) endif endif @@ -5868,13 +5860,13 @@ subroutine adjustSegmentEtaToFitBathymetry(G, GV, US, segment, fld, at_node) ! call sum_across_PEs(contractions) ! if ((contractions > 0) .and. (is_root_pe())) then ! write(mesg,'("Thickness OBCs were contracted ",'// & - ! '"to fit topography in ",I8," places.")') contractions + ! '"to fit topography in ",I0," places.")') contractions ! call MOM_error(WARNING, 'adjustEtaToFitBathymetry: '//mesg) ! endif ! call sum_across_PEs(dilations) ! if ((dilations > 0) .and. (is_root_pe())) then ! write(mesg,'("Thickness OBCs were dilated ",'// & - ! '"to fit topography in ",I8," places.")') dilations + ! '"to fit topography in ",I0," places.")') dilations ! call MOM_error(WARNING, 'adjustEtaToFitBathymetry: '//mesg) ! endif @@ -6072,39 +6064,7 @@ subroutine rotate_OBC_segment_config(segment_in, G_in, segment, G, turns) segment%open = segment_in%open segment%gradient = segment_in%gradient - if ((qturns == 0) .or. (qturns == 2)) then - segment%u_values_needed = segment_in%u_values_needed - segment%v_values_needed = segment_in%v_values_needed - segment%uamp_values_needed = segment_in%uamp_values_needed - segment%vamp_values_needed = segment_in%vamp_values_needed - segment%uphase_values_needed = segment_in%uphase_values_needed - segment%vphase_values_needed = segment_in%vphase_values_needed - segment%uamp_index = segment_in%uamp_index ! ### Perhaps this should not be set here. - segment%vamp_index = segment_in%vamp_index - segment%uphase_index = segment_in%uphase_index - segment%vphase_index = segment_in%vphase_index - else ! NOTE: [uv]_values_needed are swapped - segment%u_values_needed = segment_in%v_values_needed - segment%v_values_needed = segment_in%u_values_needed - segment%uamp_values_needed = segment_in%vamp_values_needed - segment%vamp_values_needed = segment_in%uamp_values_needed - segment%uphase_values_needed = segment_in%vphase_values_needed - segment%vphase_values_needed = segment_in%uphase_values_needed - segment%uamp_index = segment_in%vamp_index ! ### Perhaps this should not be set here. - segment%vamp_index = segment_in%uamp_index - segment%uphase_index = segment_in%vphase_index - segment%vphase_index = segment_in%uphase_index - endif - segment%z_values_needed = segment_in%z_values_needed - segment%g_values_needed = segment_in%g_values_needed - segment%t_values_needed = segment_in%t_values_needed - segment%s_values_needed = segment_in%s_values_needed - segment%zamp_values_needed = segment_in%zamp_values_needed - segment%zphase_values_needed = segment_in%zphase_values_needed - segment%zamp_index = segment_in%zamp_index ! ### Perhaps this should not be set here. - segment%zphase_index = segment_in%zphase_index - - segment%values_needed = segment_in%values_needed + call rotate_OBC_segment_values_needed(segment_in, segment, qturns) ! These are conditionally set if nudged segment%Velocity_nudging_timescale_in = segment_in%Velocity_nudging_timescale_in @@ -6229,6 +6189,51 @@ function rotate_OBC_segment_direction(direction, turns) result(rotated_dir) end function rotate_OBC_segment_direction +!> Copies which values are needed and field indices from one OBC segment type to another, +!! taking the difference in the number of turns into account. +subroutine rotate_OBC_segment_values_needed(segment_in, segment, turns) + type(OBC_segment_type), intent(in) :: segment_in !< The unrotated segment to use as a source + type(OBC_segment_type), intent(inout) :: segment !< The rotated segment to initialize + integer, intent(in) :: turns !< The number of quarter turns of the grid to apply + + integer :: qturns ! The number of quarter turns in the range of 0 to 3 + + qturns = modulo(turns, 4) + + if ((qturns == 0) .or. (qturns == 2)) then + segment%u_values_needed = segment_in%u_values_needed + segment%v_values_needed = segment_in%v_values_needed + segment%uamp_values_needed = segment_in%uamp_values_needed + segment%vamp_values_needed = segment_in%vamp_values_needed + segment%uphase_values_needed = segment_in%uphase_values_needed + segment%vphase_values_needed = segment_in%vphase_values_needed + segment%uamp_index = segment_in%uamp_index + segment%vamp_index = segment_in%vamp_index + segment%uphase_index = segment_in%uphase_index + segment%vphase_index = segment_in%vphase_index + else ! NOTE: [uv]_values_needed are swapped + segment%u_values_needed = segment_in%v_values_needed + segment%v_values_needed = segment_in%u_values_needed + segment%uamp_values_needed = segment_in%vamp_values_needed + segment%vamp_values_needed = segment_in%uamp_values_needed + segment%uphase_values_needed = segment_in%vphase_values_needed + segment%vphase_values_needed = segment_in%uphase_values_needed + segment%uamp_index = segment_in%vamp_index + segment%vamp_index = segment_in%uamp_index + segment%uphase_index = segment_in%vphase_index + segment%vphase_index = segment_in%uphase_index + endif + segment%z_values_needed = segment_in%z_values_needed + segment%g_values_needed = segment_in%g_values_needed + segment%t_values_needed = segment_in%t_values_needed + segment%s_values_needed = segment_in%s_values_needed + segment%zamp_values_needed = segment_in%zamp_values_needed + segment%zphase_values_needed = segment_in%zphase_values_needed + segment%zamp_index = segment_in%zamp_index + segment%zphase_index = segment_in%zphase_index + segment%values_needed = segment_in%values_needed + +end subroutine rotate_OBC_segment_values_needed !> Rotate an OBC segment's fields from the input to the model index map. subroutine rotate_OBC_segment_data(segment_in, segment, turns) @@ -6250,19 +6255,7 @@ subroutine rotate_OBC_segment_data(segment_in, segment, turns) IsdB = segment%HI%IsdB ; IedB = segment%HI%IedB JsdB = segment%HI%JsdB ; JedB = segment%HI%JedB - if ((turns == 0) .or. (turns == 2)) then - segment%uamp_index = segment_in%uamp_index - segment%vamp_index = segment_in%vamp_index - segment%uphase_index = segment_in%uphase_index - segment%vphase_index = segment_in%vphase_index - else ! NOTE: [uv]_values_needed are swapped - segment%uamp_index = segment_in%vamp_index - segment%vamp_index = segment_in%uamp_index - segment%uphase_index = segment_in%vphase_index - segment%vphase_index = segment_in%uphase_index - endif - segment%zamp_index = segment_in%zamp_index - segment%zphase_index = segment_in%zphase_index + call rotate_OBC_segment_values_needed(segment_in, segment, turns) segment%num_fields = segment_in%num_fields do n = 1, num_fields @@ -6466,7 +6459,7 @@ subroutine write_OBC_info(OBC, G, GV, US) turns = modulo(G%HI%turns, 4) - write(mesg, '("OBC has ", I3, " segments.")') OBC%number_of_segments + write(mesg, '("OBC has ", I0, " segments.")') OBC%number_of_segments call MOM_mesg(mesg, verb=1) ! call MOM_error(WARNING, mesg) @@ -6516,7 +6509,7 @@ subroutine write_OBC_info(OBC, G, GV, US) if (OBC%debug) call MOM_mesg("debug", verb=1) if (OBC%ramp) call MOM_mesg("ramp", verb=1) if (OBC%ramping_is_activated) call MOM_mesg("ramping_is_activated", verb=1) - write(mesg, '("n_tide_constituents ", I3)') OBC%n_tide_constituents + write(mesg, '("n_tide_constituents ", I0)') OBC%n_tide_constituents call MOM_mesg(mesg, verb=1) if (OBC%n_tide_constituents > 0) then do c=1,OBC%n_tide_constituents @@ -6542,18 +6535,18 @@ subroutine write_OBC_info(OBC, G, GV, US) dir = segment%direction unrot_dir = rotate_OBC_segment_direction(dir, -turns) - write(mesg, '(" Segment ", I3, " has direction ", I3)') n, unrot_dir - if (unrot_dir == OBC_DIRECTION_N) write(mesg, '(" Segment ", I3, " is Northern")') n - if (unrot_dir == OBC_DIRECTION_S) write(mesg, '(" Segment ", I3, " is Southern")') n - if (unrot_dir == OBC_DIRECTION_E) write(mesg, '(" Segment ", I3, " is Eastern")') n - if (unrot_dir == OBC_DIRECTION_W) write(mesg, '(" Segment ", I3, " is Western")') n + write(mesg, '(" Segment ", I0, " has direction ", I0)') n, unrot_dir + if (unrot_dir == OBC_DIRECTION_N) write(mesg, '(" Segment ", I0, " is Northern")') n + if (unrot_dir == OBC_DIRECTION_S) write(mesg, '(" Segment ", I0, " is Southern")') n + if (unrot_dir == OBC_DIRECTION_E) write(mesg, '(" Segment ", I0, " is Eastern")') n + if (unrot_dir == OBC_DIRECTION_W) write(mesg, '(" Segment ", I0, " is Western")') n call MOM_mesg(mesg, verb=1) - ! write(mesg, '(" range: ", 4I3)') segment%Is_obc, segment%Ie_obc, segment%Js_obc, segment%Je_obc + ! write(mesg, '(" range:", 4(x,I0))') segment%Is_obc, segment%Ie_obc, segment%Js_obc, segment%Je_obc if (modulo(turns, 2) == 0) then - write(mesg, '(" size: ", 4I3)') 1+abs(segment%Ie_obc-segment%Is_obc), 1+abs(segment%Je_obc-segment%Js_obc) + write(mesg, '(" size:", 2(x,I0))') 1+abs(segment%Ie_obc-segment%Is_obc), 1+abs(segment%Je_obc-segment%Js_obc) else - write(mesg, '(" size: ", 4I3)') 1+abs(segment%Je_obc-segment%Js_obc), 1+abs(segment%Ie_obc-segment%Is_obc) + write(mesg, '(" size:", 2(x,I0))') 1+abs(segment%Je_obc-segment%Js_obc), 1+abs(segment%Ie_obc-segment%Is_obc) endif call MOM_mesg(mesg, verb=1) @@ -6615,31 +6608,46 @@ subroutine write_OBC_info(OBC, G, GV, US) end subroutine write_OBC_info -!> Write checksums and perhaps the values of all the allocated arrays on an OBC segments. +!> Write checksums and perhaps some or all of the values of all the allocated arrays on the OBC segments. subroutine chksum_OBC_segments(OBC, G, GV, US, nk) - type(ocean_OBC_type), pointer :: OBC !< OBC on input map + type(ocean_OBC_type), intent(in) :: OBC !< An open boundary condition control structure type(ocean_grid_type), intent(in) :: G !< Rotated grid metric type(verticalGrid_type), intent(in) :: GV !< Vertical grid type(unit_scale_type), intent(in) :: US !< Unit scaling integer, intent(in) :: nk !< The number of layers to print ! Local variables - type(OBC_segment_type), pointer :: segment => NULL() ! pointer to segment type list + integer :: n ! The segment number reported in output + integer :: n_seg ! The internal segment number + + do n=1,OBC%number_of_segments + n_seg = n ; if (OBC%reverse_segment_order) n_seg = OBC%number_of_segments + 1 - n + + call chksum_OBC_segment_data(OBC%segment(n_seg), GV, US, nk, n) + enddo + +end subroutine chksum_OBC_segments + + +!> Write checksums and perhaps some or all of the values of all the allocated arrays on a single OBC segment. +subroutine chksum_OBC_segment_data(segment, GV, US, nk, nseg_out) + type(OBC_segment_type), intent(in) :: segment !< Segment type to checksum + type(verticalGrid_type), intent(in) :: GV !< Vertical grid + type(unit_scale_type), intent(in) :: US !< Unit scaling + integer, intent(in) :: nk !< The number of layers to print + integer, intent(in) :: nseg_out !< The segment number reported in output + + ! Local variables real :: norm ! A sign change used when rotating a normal component [nondim] real :: tang ! A sign change used when rotating a tangential component [nondim] character(len=8) :: sn, segno character(len=1024) :: mesg - integer :: n ! The segment number reported in output - integer :: n_seg ! The internal segment number integer :: dir ! This indicates the internal logical orientation of a segment - do n=1,OBC%number_of_segments - n_seg = n ; if (OBC%reverse_segment_order) n_seg = OBC%number_of_segments + 1 - n - segment => OBC%segment(n_seg) dir = segment%direction - write(segno, '(I3)') n - sn = '('//trim(adjustl(segno))//')' + write(segno, '(I0)') nseg_out + sn = '('//trim(segno)//')' ! Turn each segment and write it as though it is an eastern face. norm = 0.0 ; tang = 0.0 @@ -6704,7 +6712,6 @@ subroutine chksum_OBC_segments(OBC, G, GV, US, nk) if (allocated(segment%nudged_tangential_grad)) & call write_3d_array_vals("nudged_tangential_grad"//trim(sn), segment%nudged_tangential_grad, dir, nk, & unscale=tang*norm*US%s_to_T) - enddo contains @@ -6809,7 +6816,7 @@ subroutine write_3d_array_vals(name, Array, seg_dir, nkp, unscale) end subroutine write_3d_array_vals -end subroutine chksum_OBC_segments +end subroutine chksum_OBC_segment_data !> \namespace mom_open_boundary !! This module implements some aspects of internal open boundary From 04ab2bd4f7663dcf209a5a61caf077bdbcd9716d Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Sun, 24 Aug 2025 11:12:22 -0400 Subject: [PATCH 4/6] +Call initialize_segment_data with rotated OBCs Call initialize_segment_data() with the rotated OBC type. This includes the addition of a new turns arguments to initialize_segment_data() and the elimination of the ocean_grid_type argument to this routine. This includes the addition of the new function rotated_field_name() to swap the names of u- and v- velocity fields when the grid is rotated, avoiding the duplication of code. Internally in initialize_segment_data(), num_fields was renamed to reflect the fact that it is actually the number of tracer or other fields that are initialized via the OBC manifest string (which does not include generic OBGC tracers), and not the total number of fields on a segment. Code was also added to allow for rotate_OBC_config() to be called either before or after initialize_segment_data(). With these changes, the call to initialize_segment_data() can be moved much later in the order of the initialization calls to address issues with the initialization of OBGC tracers with open boundary conditions. All answers are bitwise identical, but there are changes to the arguments of the publicly visible routine initialize_segment_data(). --- src/core/MOM.F90 | 25 ++++-- src/core/MOM_open_boundary.F90 | 136 ++++++++++++++++++--------------- 2 files changed, 91 insertions(+), 70 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 001be33812..a7a41e85f4 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -2803,6 +2803,8 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, & call get_param(param_file, "MOM", "INDEX_TURNS", turns, & "Number of counterclockwise quarter-turn index rotations.", & default=1, debuggingParam=.true.) + else + turns = 0 endif ! Set up the model domain and grids. @@ -2873,13 +2875,13 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, & endif CS%HFrz = (US%Z_to_m * GV%m_to_H) * HFrz_z - if (associated(OBC_in)) then - ! These calls allocate the arrays on the segments for open boundary data and initialize the - ! relevant vertical remapping structures. They can only occur after the vertical grid has been - ! initialized. - call open_boundary_setup_vert(GV, US, OBC_in) - call initialize_segment_data(G_in, GV, US, OBC_in, param_file) - endif +! if (associated(OBC_in)) then +! ! These calls allocate the arrays on the segments for open boundary data and initialize the +! ! relevant vertical remapping structures. They can only occur after the vertical grid has been +! ! initialized. +! call open_boundary_setup_vert(GV, US, OBC_in) +! call initialize_segment_data(GV, US, OBC_in, param_file, turns=0) +! endif ! 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. @@ -2909,6 +2911,13 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, & endif ! dG_in is retained for now so that it can be used with write_ocean_geometry_file() below. + if (associated(CS%OBC)) then + ! This call allocates the arrays on the segments for open boundary data and initializes the + ! relevant vertical remapping structures. + call open_boundary_setup_vert(GV, US, CS%OBC) + call initialize_segment_data(GV, US, CS%OBC, param_file, turns) + endif + if (is_root_PE()) call check_MOM6_scaling_factors(CS%GV, US) call callTree_waypoint("grids initialized (initialize_MOM)") @@ -3087,7 +3096,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, & ! could occur with the call to update_OBC_data or after the main initialization. if (use_temperature) & call register_temp_salt_segments(GV, US, CS%OBC, CS%tracer_Reg, param_file) - !This is the equivalent call to register_temp_salt_segments for external tracers with OBC + ! This is the equivalent call to register_temp_salt_segments for external tracers with OBC call call_tracer_register_obc_segments(GV, param_file, CS%tracer_flow_CSp, CS%tracer_Reg, CS%OBC) ! This needs the number of tracers and to have called any code that sets whether diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index c3f02a0e9b..3d42b74ac0 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -314,6 +314,7 @@ module MOM_open_boundary logical :: zero_biharmonic = .false. !< If True, zeros the Laplacian of flow on open boundaries for !! use in the biharmonic viscosity term. logical :: brushcutter_mode = .false. !< If True, read data on supergrid. + logical :: initialize_segment_data_called = .false. !< This is set to true when initialize_segment_data is called. logical, allocatable :: tracer_x_reservoirs_used(:) !< Dimensioned by the number of tracers, set globally, !! true for those with x reservoirs (needed for restarts). logical, allocatable :: tracer_y_reservoirs_used(:) !< Dimensioned by the number of tracers, set globally, @@ -839,31 +840,29 @@ subroutine open_boundary_setup_vert(GV, US, OBC) end subroutine open_boundary_setup_vert - !> Get and store properties about the fields on the OBC segments and allocate space for reading !! OBC data from files. In the process, it does funky stuff with the MPI processes. -subroutine initialize_segment_data(G, GV, US, OBC, PF) - type(ocean_grid_type), intent(in) :: G !< Ocean grid structure +subroutine initialize_segment_data(GV, US, OBC, PF, turns) 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), target, intent(inout) :: OBC !< Open boundary control structure type(param_file_type), intent(in) :: PF !< Parameter file handle + integer, intent(in) :: turns !< Number of quarter turns of the grid - integer :: n, n_seg, m, num_fields, mm + integer :: n, n_seg, m, num_manifest_fields, mm character(len=1024) :: segstr character(len=256) :: filename character(len=20) :: segname, suffix character(len=32) :: fieldname real :: value ! A value that is parsed from the segment data string [various units] - real :: dz_neglect, dz_neglect_edge ! Small thicknesses in vertical height units [Z ~> m] character(len=32), dimension(MAX_OBC_FIELDS) :: fields ! segment field names character(len=128) :: inputdir type(OBC_segment_type), pointer :: segment => NULL() ! pointer to segment type list character(len=256) :: mesg ! Message for error messages. integer, dimension(4) :: siz - integer :: is, ie, js, je integer :: isd, ied, jsd, jed integer :: IsdB, IedB, JsdB, JedB + integer :: qturns ! The number of quarter turns in the range of 0 to 3 integer, dimension(:), allocatable :: saved_pelist integer :: current_pe integer, dimension(1) :: single_pelist @@ -871,7 +870,9 @@ subroutine initialize_segment_data(G, GV, US, OBC, PF) !will be able to dynamically switch between sub-sampling refined grid data or model grid integer :: IO_needs(2) ! Sums to determine global OBC data use and update patterns. - is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec + qturns = modulo(turns, 4) + + OBC%initialize_segment_data_called = .true. ! ###IS THIS STILL TRUE? ! There is a problem with the order of the OBC initialization @@ -917,14 +918,15 @@ subroutine initialize_segment_data(G, GV, US, OBC, PF) call MOM_error(FATAL, mesg) endif - call parse_segment_manifest_str(trim(segstr), num_fields, fields) - if (num_fields == 0) then + call parse_segment_manifest_str(trim(segstr), num_manifest_fields, fields) + !There are OBC%num_obgc_tracers obgc tracers that are not listed in param file + segment%num_fields = num_manifest_fields + OBC%num_obgc_tracers + + if (segment%num_fields == 0) then call MOM_mesg('initialize_segment_data: num_fields = 0') cycle ! cycle to next segment endif - !There are OBC%num_obgc_tracers obgc tracers are there that are not listed in param file - segment%num_fields = num_fields + OBC%num_obgc_tracers allocate(segment%field(segment%num_fields)) segment%temp_segment_data_exists = .false. @@ -933,15 +935,13 @@ subroutine initialize_segment_data(G, GV, US, OBC, PF) ! CODE HERE FOR OTHER OPTIONS (CLAMPED, NUDGED,..) !! - isd = segment%HI%isd ; ied = segment%HI%ied - jsd = segment%HI%jsd ; jed = segment%HI%jed - IsdB = segment%HI%IsdB ; IedB = segment%HI%IedB - JsdB = segment%HI%JsdB ; JedB = segment%HI%JedB + isd = segment%HI%isd ; ied = segment%HI%ied ; IsdB = segment%HI%IsdB ; IedB = segment%HI%IedB + jsd = segment%HI%jsd ; jed = segment%HI%jed ; JsdB = segment%HI%JsdB ; JedB = segment%HI%JedB obgc_segments_props_list => OBC%obgc_segments_props !pointer to the head node do m=1,segment%num_fields - if (m <= num_fields) then + if (m <= num_manifest_fields) then ! These are tracers with segments specified in MOM6 style override files call parse_segment_data_str(trim(segstr), m, trim(fields(m)), value, filename, fieldname) segment%field(m)%genre = '' @@ -954,7 +954,7 @@ subroutine initialize_segment_data(G, GV, US, OBC, PF) call get_obgc_segments_props(obgc_segments_props_list, fields(m), filename, fieldname, & segment%field(m)%resrv_lfac_in, segment%field(m)%resrv_lfac_out) ! Make sure the obgc tracer is not specified in the MOM6 param file too. - do mm=1,num_fields + do mm=1,num_manifest_fields if (trim(fields(m)) == trim(fields(mm))) then if (is_root_pe()) & call MOM_error(FATAL,"MOM_open_boundary:initialize_segment_data(): obgc tracer " //trim(fields(m))// & @@ -963,11 +963,12 @@ subroutine initialize_segment_data(G, GV, US, OBC, PF) enddo endif - segment%field(m)%name = trim(fields(m)) + segment%field(m)%name = rotated_field_name(trim(fields(m)), turns) + ! The scale factor for tracers may also be set in register_segment_tracer, and a constant input ! value is rescaled there. segment%field(m)%scale = scale_factor_from_name(fields(m), GV, US, segment%tr_Reg) - segment%field(m)%on_face = field_is_on_face(fields(m), segment%is_E_or_W) + segment%field(m)%on_face = field_is_on_face(segment%field(m)%name, segment%is_E_or_W) if (trim(filename) /= 'none') then OBC%update_OBC = .true. ! Data is assumed to be time-dependent if we are reading from file @@ -1033,6 +1034,13 @@ subroutine initialize_segment_data(G, GV, US, OBC, PF) else ! This data is not being read from a file. segment%field(m)%value = segment%field(m)%scale * value + ! Change the sign of the specified velocities, depending on the number of quarter turns of the grid. + if ( ( ((segment%field(m)%name == 'U') .or. (segment%field(m)%name == 'Uamp')) .and. & + ((qturns == 1) .or. (qturns == 2)) ) .or. & + ( ((segment%field(m)%name == 'V') .or. (segment%field(m)%name == 'Vamp')) .and. & + ((qturns == 3) .or. (qturns == 2)) ) ) & + segment%field(m)%value = -segment%field(m)%value + segment%field(m)%use_IO = .false. ! Check if this is a tidal field. If so, the number @@ -1078,6 +1086,10 @@ subroutine initialize_segment_data(G, GV, US, OBC, PF) write(mesg,'("Values needed for OBC segment ",I0)') n call MOM_error(FATAL, mesg) endif + + ! write(stderr, '(A)') trim(suffix)//" segment checksum" + if (OBC%debug) call chksum_OBC_segment_data(OBC%segment(n_seg), GV, US, OBC%nk_OBC_debug, n) + enddo call Set_PElist(saved_pelist) @@ -3872,7 +3884,6 @@ subroutine deallocate_OBC_segment_data(segment) if (associated(segment%tr_Reg)) call segment_tracer_registry_end(segment%tr_Reg) - end subroutine deallocate_OBC_segment_data !> Set tangential velocities outside of open boundaries to silly values @@ -5911,7 +5922,8 @@ subroutine rotate_OBC_config(OBC_in, G_in, OBC, G, turns) ! Data stored in setup_[uv]_point_obc is needed for allocate_obc_segment_data call allocate_OBC_segment_data(OBC, OBC%segment(l_seg)) ! Initialize the field-related data of a rotated segment. - call rotate_OBC_segment_data(OBC_in%segment(l_seg), OBC%segment(l_seg), turns) + if (OBC%initialize_segment_data_called) & + call rotate_OBC_segment_data(OBC_in%segment(l_seg), OBC%segment(l_seg), turns) enddo ! The horizontal segment map @@ -6013,14 +6025,18 @@ subroutine rotate_OBC_config(OBC_in, G_in, OBC, G, turns) OBC%tidal_longitudes = OBC_in%tidal_longitudes endif - ! 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 + ! remap_z_CS and remap_h_CS are set up by initialize_segment_data, so we copy the types here + ! if they are already associated. + if (OBC%initialize_segment_data_called) then + ! This block of code is unnecessary if initialize_segment_data is called after rotate_OBC_config. + 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 endif end subroutine rotate_OBC_config @@ -6245,16 +6261,10 @@ subroutine rotate_OBC_segment_data(segment_in, segment, turns) logical :: flip_normal_vel_sign, flip_tang_vel_sign integer :: n integer :: num_fields - integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, ke num_fields = segment_in%num_fields allocate(segment%field(num_fields)) - isd = segment%HI%isd ; ied = segment%HI%ied - jsd = segment%HI%jsd ; jed = segment%HI%jed - IsdB = segment%HI%IsdB ; IedB = segment%HI%IedB - JsdB = segment%HI%JsdB ; JedB = segment%HI%JedB - call rotate_OBC_segment_values_needed(segment_in, segment, turns) segment%num_fields = segment_in%num_fields @@ -6274,30 +6284,7 @@ subroutine rotate_OBC_segment_data(segment_in, segment, turns) call rotate_array(segment_in%field(n)%buffer_dst, turns, segment%field(n)%buffer_dst) endif - if (modulo(turns, 2) /= 0) then - select case (segment_in%field(n)%name) - case ('U') - segment%field(n)%name = 'V' - case ('Uamp') - segment%field(n)%name = 'Vamp' - case ('Uphase') - segment%field(n)%name = 'Vphase' - case ('V') - segment%field(n)%name = 'U' - case ('Vamp') - segment%field(n)%name = 'Uamp' - case ('Vphase') - segment%field(n)%name = 'Uphase' - case ('DVDX') - segment%field(n)%name = 'DUDY' - case ('DUDY') - segment%field(n)%name = 'DVDX' - case default - segment%field(n)%name = segment_in%field(n)%name - end select - else - segment%field(n)%name = segment_in%field(n)%name - endif + segment%field(n)%name = rotated_field_name(segment_in%field(n)%name, turns) if (allocated(segment_in%field(n)%buffer_src)) then call allocate_rotated_seg_data(segment_in%field(n)%buffer_src, segment_in%HI, & @@ -6407,9 +6394,34 @@ subroutine rotate_OBC_segment_data(segment_in, segment, turns) segment%temp_segment_data_exists = segment_in%temp_segment_data_exists segment%salt_segment_data_exists = segment_in%salt_segment_data_exists + end subroutine rotate_OBC_segment_data +!> Return the that the field would have after being rotated by the given number of quarter turns +function rotated_field_name(input_name, turns) + character(len=*), intent(in) :: input_name !< The unrotated field name + integer, intent(in) :: turns !< Number of quarter turns of the grid + character(len=len(input_name)) :: rotated_field_name !< The rotated field name + + if (modulo(turns, 2) /= 0) then + select case (input_name) + case ('U') ; rotated_field_name = 'V' + case ('Uamp') ; rotated_field_name = 'Vamp' + case ('Uphase') ; rotated_field_name = 'Vphase' + case ('V') ; rotated_field_name = 'U' + case ('Vamp') ; rotated_field_name = 'Uamp' + case ('Vphase') ; rotated_field_name = 'Uphase' + case ('DVDX') ; rotated_field_name = 'DUDY' + case ('DUDY') ; rotated_field_name = 'DVDX' + case default ; rotated_field_name = input_name + end select + else + rotated_field_name = input_name + endif + +end function rotated_field_name + !> Allocate an array of data for a field on a segment based on the size of a potentially rotated source array subroutine allocate_rotated_seg_data(src_array, HI_in, tgt_array, segment) real, dimension(:,:,:), intent(in) :: src_array !< The segment data on the unrotated source grid @@ -6542,11 +6554,11 @@ subroutine write_OBC_info(OBC, G, GV, US) if (unrot_dir == OBC_DIRECTION_W) write(mesg, '(" Segment ", I0, " is Western")') n call MOM_mesg(mesg, verb=1) - ! write(mesg, '(" range:", 4(x,I0))') segment%Is_obc, segment%Ie_obc, segment%Js_obc, segment%Je_obc + ! write(mesg, '(" range:", 4(1x,I0))') segment%Is_obc, segment%Ie_obc, segment%Js_obc, segment%Je_obc if (modulo(turns, 2) == 0) then - write(mesg, '(" size:", 2(x,I0))') 1+abs(segment%Ie_obc-segment%Is_obc), 1+abs(segment%Je_obc-segment%Js_obc) + write(mesg, '(" size: ", I0," ",I0)') 1+abs(segment%Ie_obc-segment%Is_obc), 1+abs(segment%Je_obc-segment%Js_obc) else - write(mesg, '(" size:", 2(x,I0))') 1+abs(segment%Je_obc-segment%Js_obc), 1+abs(segment%Ie_obc-segment%Is_obc) + write(mesg, '(" size: ", I0," ",I0)') 1+abs(segment%Je_obc-segment%Js_obc), 1+abs(segment%Ie_obc-segment%Is_obc) endif call MOM_mesg(mesg, verb=1) From badeacf3daf9232aca95046f8ad1b8a39ed5acff Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Sun, 31 Aug 2025 06:18:53 -0400 Subject: [PATCH 5/6] *Call initialize_segment_data after call_tracer_register_obc_segments Move call to initialize_segment_data() after the call to call_tracer_register_obc_segments() to correct an answer changing bug in models with OBGC tracers and open boundaries that was causing the dev/CEFI branch answers not to reproduce with dev/gfdl. This commit will change answers (and restore previous answers) in cases that use open boundary conditions with biogeochemical tracers, such as those in COBALT. --- src/core/MOM.F90 | 31 +++++++++++-------------------- 1 file changed, 11 insertions(+), 20 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index a7a41e85f4..fc874c8cbd 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -116,9 +116,9 @@ module MOM use MOM_open_boundary, only : register_temp_salt_segments, update_segment_tracer_reservoirs use MOM_open_boundary, only : setup_OBC_tracer_reservoirs use MOM_open_boundary, only : open_boundary_register_restarts, remap_OBC_fields -use MOM_open_boundary, only : open_boundary_setup_vert, initialize_segment_data, rotate_OBC_config -use MOM_open_boundary, only : update_OBC_segment_data, open_boundary_halo_update -use MOM_open_boundary, only : write_OBC_info, chksum_OBC_segments +use MOM_open_boundary, only : open_boundary_setup_vert, initialize_segment_data +use MOM_open_boundary, only : update_OBC_segment_data, rotate_OBC_config +use MOM_open_boundary, only : open_boundary_halo_update, write_OBC_info, chksum_OBC_segments use MOM_porous_barriers, only : porous_widths_layer, porous_widths_interface, porous_barriers_init use MOM_porous_barriers, only : porous_barrier_CS use MOM_set_visc, only : set_viscous_BBL, set_viscous_ML, set_visc_CS @@ -2614,7 +2614,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, & "The time between OBC segment data updates for OBGC tracers. "//& "This must be an integer multiple of DT and DT_THERM. "//& "The default is set to DT.", & - units="s", default=US%T_to_s*CS%dt, scale=US%s_to_T, do_not_log=.not.associated(CS%OBC)) + units="s", default=US%T_to_s*CS%dt, scale=US%s_to_T, do_not_log=.not.associated(OBC_in)) ! This is here in case these values are used inappropriately. use_frazil = .false. ; bound_salinity = .false. ; use_p_surf_in_EOS = .false. @@ -2875,14 +2875,6 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, & endif CS%HFrz = (US%Z_to_m * GV%m_to_H) * HFrz_z -! if (associated(OBC_in)) then -! ! These calls allocate the arrays on the segments for open boundary data and initialize the -! ! relevant vertical remapping structures. They can only occur after the vertical grid has been -! ! initialized. -! call open_boundary_setup_vert(GV, US, OBC_in) -! call initialize_segment_data(GV, US, OBC_in, param_file, turns=0) -! endif - ! 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 @@ -2911,13 +2903,6 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, & endif ! dG_in is retained for now so that it can be used with write_ocean_geometry_file() below. - if (associated(CS%OBC)) then - ! This call allocates the arrays on the segments for open boundary data and initializes the - ! relevant vertical remapping structures. - call open_boundary_setup_vert(GV, US, CS%OBC) - call initialize_segment_data(GV, US, CS%OBC, param_file, turns) - endif - if (is_root_PE()) call check_MOM6_scaling_factors(CS%GV, US) call callTree_waypoint("grids initialized (initialize_MOM)") @@ -3088,6 +3073,9 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, & CS%mixedlayer_restrat_CSp, restart_CSp) if (associated(CS%OBC)) then + ! This call initializes the relevant vertical remapping structures. + call open_boundary_setup_vert(GV, US, CS%OBC) + ! Set up remaining information about open boundary conditions that is needed for OBCs. ! Package specific changes to OBCs occur here. call call_OBC_register(G, GV, US, param_file, CS%update_OBC_CSp, CS%OBC, CS%tracer_Reg) @@ -3104,10 +3092,13 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, & call open_boundary_register_restarts(HI, GV, US, CS%OBC, CS%tracer_Reg, & param_file, restart_CSp, use_temperature) + ! This call allocates the arrays on the segments for open boundary data, but it must occur + ! after any calls to call_tracer_register_obc_segments. + call initialize_segment_data(GV, US, CS%OBC, param_file, turns) + if (CS%debug_OBCs) call write_OBC_info(CS%OBC, G, GV, US) endif - if (present(waves_CSp)) then call waves_register_restarts(waves_CSp, HI, GV, US, param_file, restart_CSp) endif From 0fac41fa8b4d12c4f48a6b56fe937e2cdb001374 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Sun, 31 Aug 2025 07:15:38 -0400 Subject: [PATCH 6/6] Eliminate rotate_OBC_segment_data The temporary blocks of code allowing for initialize_segment_data either before or after rotate_OBC_config() were eliminated, as was the now unused subroutine rotate_OBC_segment_data(). All of the functionality of rotate_OBC_segment_data had already been merged into initialize_segment_data. All answers are bitwise identical, but there is one less internal subroutine, and the order of calls to initialize the open boundary conditions is now less flexible. --- src/core/MOM_open_boundary.F90 | 172 --------------------------------- 1 file changed, 172 deletions(-) diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index 3d42b74ac0..9470ad1ce8 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -314,7 +314,6 @@ module MOM_open_boundary logical :: zero_biharmonic = .false. !< If True, zeros the Laplacian of flow on open boundaries for !! use in the biharmonic viscosity term. logical :: brushcutter_mode = .false. !< If True, read data on supergrid. - logical :: initialize_segment_data_called = .false. !< This is set to true when initialize_segment_data is called. logical, allocatable :: tracer_x_reservoirs_used(:) !< Dimensioned by the number of tracers, set globally, !! true for those with x reservoirs (needed for restarts). logical, allocatable :: tracer_y_reservoirs_used(:) !< Dimensioned by the number of tracers, set globally, @@ -872,14 +871,6 @@ subroutine initialize_segment_data(GV, US, OBC, PF, turns) qturns = modulo(turns, 4) - OBC%initialize_segment_data_called = .true. - - ! ###IS THIS STILL TRUE? - ! There is a problem with the order of the OBC initialization - ! with respect to ALE_init. Currently handling this by copying the - ! param file so that I can use it later in step_MOM in order to finish - ! initializing segments on the first step. - call get_param(PF, mdl, "INPUTDIR", inputdir, default=".") inputdir = slasher(inputdir) @@ -5921,9 +5912,6 @@ subroutine rotate_OBC_config(OBC_in, G_in, OBC, G, turns) call rotate_OBC_segment_config(OBC_in%segment(l_seg), G_in, OBC%segment(l_seg), G, turns) ! Data stored in setup_[uv]_point_obc is needed for allocate_obc_segment_data call allocate_OBC_segment_data(OBC, OBC%segment(l_seg)) - ! Initialize the field-related data of a rotated segment. - if (OBC%initialize_segment_data_called) & - call rotate_OBC_segment_data(OBC_in%segment(l_seg), OBC%segment(l_seg), turns) enddo ! The horizontal segment map @@ -6025,20 +6013,6 @@ subroutine rotate_OBC_config(OBC_in, G_in, OBC, G, turns) OBC%tidal_longitudes = OBC_in%tidal_longitudes endif - ! remap_z_CS and remap_h_CS are set up by initialize_segment_data, so we copy the types here - ! if they are already associated. - if (OBC%initialize_segment_data_called) then - ! This block of code is unnecessary if initialize_segment_data is called after rotate_OBC_config. - 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 - endif - end subroutine rotate_OBC_config !> Rotate the OBC segment configuration data from the input to model index map. @@ -6251,152 +6225,6 @@ subroutine rotate_OBC_segment_values_needed(segment_in, segment, turns) end subroutine rotate_OBC_segment_values_needed -!> Rotate an OBC segment's fields from the input to the model index map. -subroutine rotate_OBC_segment_data(segment_in, segment, turns) - type(OBC_segment_type), intent(in) :: segment_in !< The unrotated segment to use as a source - type(OBC_segment_type), intent(inout) :: segment !< The rotated segment to initialize - integer, intent(in) :: turns !< The number of quarter turns of the grid to apply - - ! Local variables - logical :: flip_normal_vel_sign, flip_tang_vel_sign - integer :: n - integer :: num_fields - - num_fields = segment_in%num_fields - allocate(segment%field(num_fields)) - - call rotate_OBC_segment_values_needed(segment_in, segment, turns) - - segment%num_fields = segment_in%num_fields - do n = 1, num_fields - segment%field(n)%handle = segment_in%field(n)%handle - segment%field(n)%dz_handle = segment_in%field(n)%dz_handle - segment%field(n)%use_IO = segment_in%field(n)%use_IO - segment%field(n)%genre = segment_in%field(n)%genre - segment%field(n)%scale = segment_in%field(n)%scale - segment%field(n)%resrv_lfac_in = segment_in%field(n)%resrv_lfac_in - segment%field(n)%resrv_lfac_out = segment_in%field(n)%resrv_lfac_out - segment%field(n)%on_face = segment_in%field(n)%on_face - - if (allocated(segment_in%field(n)%buffer_dst)) then - call allocate_rotated_seg_data(segment_in%field(n)%buffer_dst, segment_in%HI, & - segment%field(n)%buffer_dst, segment) - call rotate_array(segment_in%field(n)%buffer_dst, turns, segment%field(n)%buffer_dst) - endif - - segment%field(n)%name = rotated_field_name(segment_in%field(n)%name, turns) - - if (allocated(segment_in%field(n)%buffer_src)) then - call allocate_rotated_seg_data(segment_in%field(n)%buffer_src, segment_in%HI, & - segment%field(n)%buffer_src, segment) - call rotate_array(segment_in%field(n)%buffer_src, turns, segment%field(n)%buffer_src) - endif - - segment%field(n)%nk_src = segment_in%field(n)%nk_src - - if (allocated(segment_in%field(n)%dz_src)) then - call allocate_rotated_seg_data(segment_in%field(n)%dz_src, segment_in%HI, segment%field(n)%dz_src, segment) - call rotate_array(segment_in%field(n)%dz_src, turns, segment%field(n)%dz_src) - endif - - segment%field(n)%value = segment_in%field(n)%value - enddo - - if (allocated(segment_in%SSH)) & - call rotate_array(segment_in%SSH, turns, segment%SSH) - if (allocated(segment_in%cg)) & - call rotate_array(segment_in%cg, turns, segment%cg) - if (allocated(segment_in%htot)) & - call rotate_array(segment_in%htot, turns, segment%htot) - if (allocated(segment_in%dztot)) & - call rotate_array(segment_in%dztot, turns, segment%dztot) - if (allocated(segment_in%h)) & - call rotate_array(segment_in%h, turns, segment%h) - if (allocated(segment_in%normal_vel)) & - call rotate_array(segment_in%normal_vel, turns, segment%normal_vel) - if (allocated(segment_in%normal_trans)) & - call rotate_array(segment_in%normal_trans, turns, segment%normal_trans) - if (allocated(segment_in%normal_vel_bt)) & - call rotate_array(segment_in%normal_vel_bt, turns, segment%normal_vel_bt) - if (allocated(segment_in%tangential_vel)) & - call rotate_array(segment_in%tangential_vel, turns, segment%tangential_vel) - if (allocated(segment_in%tangential_grad)) & - call rotate_array(segment_in%tangential_grad, turns, segment%tangential_grad) - if (allocated(segment_in%grad_normal)) & - call rotate_array(segment_in%grad_normal, turns, segment%grad_normal) - if (allocated(segment_in%grad_tan)) & - call rotate_array(segment_in%grad_tan, turns, segment%grad_tan) - if (allocated(segment_in%grad_gradient)) & - call rotate_array(segment_in%grad_gradient, turns, segment%grad_gradient) - if (modulo(turns, 2) /= 0) then - if (allocated(segment_in%rx_norm_rad)) & - call rotate_array(segment_in%rx_norm_rad, turns, segment%ry_norm_rad) - if (allocated(segment_in%ry_norm_rad)) & - call rotate_array(segment_in%ry_norm_rad, turns, segment%rx_norm_rad) - if (allocated(segment_in%rx_norm_obl)) & - call rotate_array(segment_in%rx_norm_obl, turns, segment%ry_norm_obl) - if (allocated(segment_in%ry_norm_obl)) & - call rotate_array(segment_in%ry_norm_obl, turns, segment%rx_norm_obl) - else - if (allocated(segment_in%rx_norm_rad)) & - call rotate_array(segment_in%rx_norm_rad, turns, segment%rx_norm_rad) - if (allocated(segment_in%ry_norm_rad)) & - call rotate_array(segment_in%ry_norm_rad, turns, segment%ry_norm_rad) - if (allocated(segment_in%rx_norm_obl)) & - call rotate_array(segment_in%rx_norm_obl, turns, segment%rx_norm_obl) - if (allocated(segment_in%ry_norm_obl)) & - call rotate_array(segment_in%ry_norm_obl, turns, segment%ry_norm_obl) - endif - if (allocated(segment_in%cff_normal)) & - call rotate_array(segment_in%cff_normal, turns, segment%cff_normal) - if (allocated(segment_in%nudged_normal_vel)) & - call rotate_array(segment_in%nudged_normal_vel, turns, segment%nudged_normal_vel) - if (allocated(segment_in%nudged_tangential_vel)) & - call rotate_array(segment_in%nudged_tangential_vel, turns, segment%nudged_tangential_vel) - if (allocated(segment_in%nudged_tangential_grad)) & - call rotate_array(segment_in%nudged_tangential_grad, turns, segment%nudged_tangential_grad) - - ! Change the sign of the normal or tangential velocities or transports that have been read in from - ! a file, depending on the orientation of the face and the number of quarter turns of the grid. - flip_normal_vel_sign = .false. ; flip_tang_vel_sign = .false. - do n = 1, num_fields - if (((segment%field(n)%name == 'U') .or. (segment%field(n)%name == 'Uamp')) .and. & - ((modulo(turns, 4) == 1) .or. (modulo(turns, 4) == 2)) ) then - if (allocated(segment%field(n)%buffer_dst)) & - segment%field(n)%buffer_dst(:,:,:) = -segment%field(n)%buffer_dst(:,:,:) - segment%field(n)%value = -segment%field(n)%value - if (segment%is_E_or_W) flip_normal_vel_sign = .true. - if (segment%is_N_or_S) flip_tang_vel_sign = .true. - elseif (((segment%field(n)%name == 'V') .or. (segment%field(n)%name == 'Vamp')) .and. & - ((modulo(turns, 4) == 3) .or. (modulo(turns, 4) == 2)) ) then - if (allocated(segment%field(n)%buffer_dst)) & - segment%field(n)%buffer_dst(:,:,:) = -segment%field(n)%buffer_dst(:,:,:) - segment%field(n)%value = -segment%field(n)%value - if (segment%is_N_or_S) flip_normal_vel_sign = .true. - if (segment%is_E_or_W) flip_tang_vel_sign = .true. - endif - enddo - - if (flip_normal_vel_sign) then - segment%normal_trans(:,:,:) = -segment%normal_trans(:,:,:) - segment%normal_vel(:,:,:) = -segment%normal_vel(:,:,:) - segment%normal_vel_bt(:,:) = -segment%normal_vel_bt(:,:) - if (allocated(segment%nudged_normal_vel)) & - segment%nudged_normal_vel(:,:,:) = -segment%nudged_normal_vel(:,:,:) - endif - - if (flip_tang_vel_sign) then - if (allocated(segment%tangential_vel)) & - segment%tangential_vel(:,:,:) = -segment%tangential_vel(:,:,:) - if (allocated(segment%nudged_tangential_vel)) & - segment%nudged_tangential_vel(:,:,:) = -segment%nudged_tangential_vel(:,:,:) - endif - - segment%temp_segment_data_exists = segment_in%temp_segment_data_exists - segment%salt_segment_data_exists = segment_in%salt_segment_data_exists - -end subroutine rotate_OBC_segment_data - !> Return the that the field would have after being rotated by the given number of quarter turns function rotated_field_name(input_name, turns)