From 22215cb0b778158cf6c56b309980900e0d9a2dca Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Tue, 14 Apr 2020 15:28:00 -0400 Subject: [PATCH] Internal field index rotation This patch introduces an index rotation feature applied to all internal fields. When this feature is enabled, model fields are initialized on the index map as specified by the configuration or input files. The values of these arrays are then moved to a rotated index map and saved to a new rotated array, which are then used for model simulation. The primary purpose of this feature is to detect regressions in rotationally symmetric calculations, such as fluid dynamic dynamic processes and other transport phenomena. Fields created or read by the drivers are unrotated and are handled as external user-defined field. External fields are rotated or unrotated as they pass between the driver and the ocean model. In order to maintain bit reproducibility, all field initialization is computed or assigned on the input grid. This patch also includes an extension of the test suite, which applies a single quarter-turn to each of the existing tests. The features of this patch are summarized below. Two new configuration features are added: - ROTATE_INDEX: Boolean flag which enables this feature. - INDEX_TURNS: Integer which descibes the number of counterclockwise quarter-turns. Default value is 1. Two new modules and several new functions are added: - MOM_array_transform: Fundamental operations for rotating arrays, vectors, rotationally equivalent pairs - MOM_transform_FMS: Wrappers to FMS operations which can only operate on the unrotated index pax. - Numerous rotation functions have been introduced in various modules associated with initialization (forcing, OBC, etc), which are documented in the source code. The following new features have been added to existing components: - paired checksums (`uvchksum`, `hchksum_pair`) now include a `scalar_pair` argument to distinguish between rotationally equivalent scalars (e.g. grid areas on u- and v-points) and vector components. - `register_restart_pair` was introduced to register rotationally associated pairs. - `turns` has been added to the horizontal index (HI) type, as a means to track the number of quarter turns in lower level operations. - `allocate_forcing_type` and `allocate_mech_forcing` have been extended to permit allocation based on an exisiting (unrotated) reference, alongside the older interface using control flags for each field type. - copy_MOM_domain includes a `turns` argument, to support rotational changes to the FMS domain. Reproducibility: - All results should be bitwise reproducible. A nonphysical exception exception is the calculation of CFL conditions in systems with zero velocity, such as the unit test case. The previous implementation could produce a "negative zero" CFL, whose sign would depend on the number of rotations. The new implementation eliminates these negative zeros, but may also change existing CFL calculations, and may report a regression in `ocean.stats`. - Rotated diagnostic checksum output is not tested, since there is not yet a way to associated rotationally equivalent diagnostics. See below for more information. Current limitations: - Only single PE runs can be rotated. Currently FMS assumes a specific ordering for assignment of PE IDs (east->west, south->north). An index rotation will typically re-orient these domains in a way which breaks this default assignment. This ID is used to manage MPI message passing targets, so any rearrangement of these IDs will cause errors. Since we currently cannot control the PE ID assigned by FMS, we are currently unable to support index rotation of parallel runs. - Certain fields and tracers cannot be rotated. There are a small subset of tracers and fields in various structures, such as boundary fluxes, which are currently unsupported. This is primarily beacuse these tracers are typically managed at the driver level. There is no reason why these fields cannot be supported, but it would require a larger test suite to detect and handle such fields. - Diagnostic output is not rotated. Diagnostics are computed on the rotated grid, but it would require additional FMS wrappers to support registering and writing diagnostics on the input index map. This is not yet implemented, so all diagnostics are on the rotated index map. One additional consequence is that there is no way to register pairs of diagnostics, and thus no way to automate the chksum.diag output. Currently, diagnostic checksums are disabled in the test suite. We hope to resolve this issue in a future patch. - Open boundary condition (OBC) rotation is restriced to INDEX_TURNS=1 Generalized OBC rotation is currently not supported. This is primarily due to the sensitivity of segment configuration, which is specified with a large number of flags denoting its orientation. It is possible to support generalized rotation, but it will require a refactor of the OBC segment code. We hope to support generalized rotation in a future patch. - Incorrect definitions of gridLonT and gridLatT Due to certain array length assumptions for 1d axes, we rotate and exchange the values of gridLatT and gridLonT during quarter-turn operations. This is not correct, since these arrays describe latitude and longitude, rather than the axes in the first and second dimension. However, several other components currently expect the first and second reference axis to have matching array lengths, so for now we swap these arrays during rotations. The 1d arrays are only used for initialization computation (which is on the input map) and I/O management, so this error does not affect any existing runs. But rotated runs which save these arrays will be incorrect. - Non-ALE sponge rotation is disabled. Only ALE sponge rotation is currently supported. It should be possible to extend rotation to include non-ALE sponges if needed. Further comments: - Fields are not rotated in-place, since many calculations require that both fields exist at any time. While it is unlikely that rotated and unrotated fields would persist for an entire run, there is likely to be an increase in memory overhead. - External fields from the drivers are rotated and unrotated at the beginning and end of every timestep, which can create additional computational overhead. Note that there are control flags, e.g. `cycle_start` and `cycle_end` which could be used to manage these operations, assuming that the driver does not change the fields inbetween coupler steps. - This patch fixes what appears to have been an error in the global h-point indexes of OBC segments. However, it also seems that these indices were unused, so this is unlikely to impact existing runs. This implementation was proposed by Robert Hallberg and follows an earlier implementation by Nic Hannah, whose efforts introduced several changes that implemented the rotational invariance of existing numerical calculations. --- .testing/Makefile | 12 +- .testing/tc2/MOM_input | 1 + src/core/MOM.F90 | 338 +++++++++-- src/core/MOM_barotropic.F90 | 34 +- src/core/MOM_dynamics_split_RK2.F90 | 40 +- src/core/MOM_forcing_type.F90 | 379 +++++++++++- src/core/MOM_open_boundary.F90 | 567 ++++++++++++++++-- src/core/MOM_transcribe_grid.F90 | 92 ++- src/core/MOM_variables.F90 | 76 +++ src/diagnostics/MOM_debugging.F90 | 4 +- src/diagnostics/MOM_sum_output.F90 | 21 +- src/framework/MOM_array_transform.F90 | 358 +++++++++++ src/framework/MOM_checksums.F90 | 492 ++++++++++++--- src/framework/MOM_domains.F90 | 73 ++- src/framework/MOM_dyn_horgrid.F90 | 2 +- src/framework/MOM_hor_index.F90 | 53 +- src/framework/MOM_horizontal_regridding.F90 | 10 +- src/framework/MOM_restart.F90 | 115 +++- src/framework/MOM_transform_FMS.F90 | 399 ++++++++++++ src/initialization/MOM_grid_initialize.F90 | 6 +- .../MOM_shared_initialization.F90 | 3 + src/parameterizations/lateral/MOM_MEKE.F90 | 9 +- .../lateral/MOM_lateral_mixing_coeffs.F90 | 9 +- .../lateral/MOM_thickness_diffuse.F90 | 3 +- .../vertical/MOM_ALE_sponge.F90 | 159 +++++ .../vertical/MOM_set_diffusivity.F90 | 13 +- .../vertical/MOM_set_viscosity.F90 | 11 +- .../vertical/MOM_vert_friction.F90 | 12 +- src/tracer/MOM_tracer_hor_diff.F90 | 6 +- 29 files changed, 2981 insertions(+), 316 deletions(-) create mode 100644 src/framework/MOM_array_transform.F90 create mode 100644 src/framework/MOM_transform_FMS.F90 diff --git a/.testing/Makefile b/.testing/Makefile index 8067e4218d..0d73979204 100644 --- a/.testing/Makefile +++ b/.testing/Makefile @@ -40,7 +40,7 @@ MKMF_TEMPLATE ?= build/mkmf/templates/ncrc-gnu.mk # Executables BUILDS = symmetric asymmetric repro openmp CONFIGS := $(wildcard tc*) -TESTS = grids layouts restarts nans dims openmps +TESTS = grids layouts restarts nans dims openmps rotations # REPRO tests enable reproducibility with optimization, and often do not match # the DEBUG results in older GCCs and vendor compilers, so we can optionally @@ -186,9 +186,15 @@ test: $(foreach t,$(TESTS),test.$(t)) # NOTE: We remove tc3 (OBC) from grid test since it cannot run asymmetric grids +# NOTE: rotation diag chksum disabled since we cannot yet compare rotationally +# equivalent diagnostics + +# TODO: restart checksum comparison is not yet implemented + .PHONY: $(foreach t,$(TESTS),test.$(t)) test.grids: $(foreach c,$(filter-out tc3,$(CONFIGS)),$(c).grid $(c).grid.diag) test.layouts: $(foreach c,$(CONFIGS),$(c).layout $(c).layout.diag) +test.rotations: $(foreach c,$(CONFIGS),$(c).rotate) test.restarts: $(foreach c,$(CONFIGS),$(c).restart) test.repros: $(foreach c,$(CONFIGS),$(c).repro $(c).repro.diag) test.openmps: $(foreach c,$(CONFIGS),$(c).openmp $(c).openmp.diag) @@ -210,6 +216,7 @@ endef $(eval $(call CMP_RULE,grid,symmetric asymmetric)) $(eval $(call CMP_RULE,layout,symmetric layout)) +$(eval $(call CMP_RULE,rotate,symmetric rotate)) $(eval $(call CMP_RULE,repro,symmetric repro)) $(eval $(call CMP_RULE,openmp,symmetric openmp)) $(eval $(call CMP_RULE,nan,symmetric nan)) @@ -260,7 +267,7 @@ results/%/ocean.stats.$(1): build/$(2)/MOM6 cp -rL $$*/* work/$$*/$(1) cd work/$$*/$(1) && if [ -f Makefile ]; then make; fi mkdir -p work/$$*/$(1)/RESTART - echo $(4) > work/$$*/$(1)/MOM_override + echo -e "$(4)" > work/$$*/$(1)/MOM_override cd work/$$*/$(1) && $$(call MPIRUN_CMD,$(5)) -n $(6) ../../../$$< 2> debug.out > std.out \ || ! sed 's/^/$$*.$(1): /' std.out debug.out \ && sed 's/^/$$*.$(1): /' std.out @@ -282,6 +289,7 @@ $(eval $(call STAT_RULE,target,target,,,,1)) $(eval $(call STAT_RULE,repro,repro,,,,1)) $(eval $(call STAT_RULE,openmp,openmp,,,,1)) $(eval $(call STAT_RULE,layout,symmetric,,LAYOUT=2$(,)1,,2)) +$(eval $(call STAT_RULE,rotate,symmetric,,ROTATE_INDEX=True\nINDEX_TURNS=1,,1)) $(eval $(call STAT_RULE,nan,symmetric,,,MALLOC_PERTURB_=256,1)) $(eval $(call STAT_RULE,dim.t,symmetric,,T_RESCALE_POWER=11,,1)) $(eval $(call STAT_RULE,dim.l,symmetric,,L_RESCALE_POWER=11,,1)) diff --git a/.testing/tc2/MOM_input b/.testing/tc2/MOM_input index c037648d95..285ee79e4b 100644 --- a/.testing/tc2/MOM_input +++ b/.testing/tc2/MOM_input @@ -600,3 +600,4 @@ ENERGYSAVEDAYS = 0.5 ! [days] default = 3600.0 ! energies of the run and other globally summed diagnostics. DIAG_AS_CHKSUM = True DEBUG = True +USE_GM_WORK_BUG = False diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 78d53e0b76..21e23d69cf 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -4,10 +4,12 @@ module MOM ! This file is part of MOM6. See LICENSE.md for the license. ! Infrastructure modules +use MOM_array_transform, only : rotate_array, rotate_vector use MOM_debugging, only : MOM_debugging_init, hchksum, uvchksum use MOM_debugging, only : check_redundant use MOM_checksum_packages, only : MOM_thermo_chksum, MOM_state_chksum use MOM_checksum_packages, only : MOM_accel_chksum, MOM_surface_chksum +use MOM_coms, only : num_PEs use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end use MOM_cpu_clock, only : CLOCK_COMPONENT, CLOCK_SUBCOMPONENT use MOM_cpu_clock, only : CLOCK_MODULE_DRIVER, CLOCK_MODULE, CLOCK_ROUTINE @@ -37,7 +39,8 @@ module MOM use MOM_io, only : MOM_io_init, vardesc, var_desc use MOM_io, only : slasher, file_exists, MOM_read_data use MOM_obsolete_params, only : find_obsolete_params -use MOM_restart, only : register_restart_field, query_initialized, save_restart +use MOM_restart, only : register_restart_field, register_restart_pair +use MOM_restart, only : query_initialized, save_restart use MOM_restart, only : restart_init, is_new_run, MOM_restart_CS use MOM_spatial_means, only : global_mass_integral use MOM_time_manager, only : time_type, real_to_time, time_type_to_real, operator(+) @@ -50,6 +53,7 @@ module MOM use MOM_ALE, only : ALE_init, ALE_end, ALE_main, ALE_CS, adjustGridForIntegrity use MOM_ALE, only : ALE_getCoordinate, ALE_getCoordinateUnits, ALE_writeCoordinateFile use MOM_ALE, only : ALE_updateVerticalGridType, ALE_remap_init_conds, ALE_register_diags +use MOM_ALE_sponge, only : rotate_ALE_sponge, update_ALE_sponge_field use MOM_barotropic, only : Barotropic_CS use MOM_boundary_update, only : call_OBC_register, OBC_register_end, update_OBC_CS use MOM_coord_initialization, only : MOM_initialize_coord @@ -72,9 +76,13 @@ module MOM use MOM_dyn_horgrid, only : dyn_horgrid_type, create_dyn_horgrid, destroy_dyn_horgrid use MOM_EOS, only : EOS_init, calculate_density, calculate_TFreeze use MOM_fixed_initialization, only : MOM_initialize_fixed +use MOM_forcing_type, only : allocate_forcing_type, allocate_mech_forcing +use MOM_forcing_type, only : deallocate_mech_forcing, deallocate_forcing_type +use MOM_forcing_type, only : rotate_forcing, rotate_mech_forcing use MOM_grid, only : ocean_grid_type, MOM_grid_init, MOM_grid_end use MOM_grid, only : set_first_direction, rescale_grid_bathymetry use MOM_hor_index, only : hor_index_type, hor_index_init +use MOM_hor_index, only : rotate_hor_index use MOM_interface_heights, only : find_eta use MOM_lateral_mixing_coeffs, only : calc_slope_functions, VarMix_init use MOM_lateral_mixing_coeffs, only : calc_resoln_function, calc_depth_function, VarMix_CS @@ -87,6 +95,7 @@ module MOM use MOM_open_boundary, only : register_temp_salt_segments use MOM_open_boundary, only : open_boundary_register_restarts use MOM_open_boundary, only : update_segment_tracer_reservoirs +use MOM_open_boundary, only : rotate_OBC_config, rotate_OBC_init use MOM_set_visc, only : set_viscous_BBL, set_viscous_ML, set_visc_init use MOM_set_visc, only : set_visc_register_restarts, set_visc_CS use MOM_sponge, only : init_sponge_diags, sponge_CS @@ -108,11 +117,13 @@ module MOM use MOM_tracer_flow_control, only : tracer_flow_control_init, call_tracer_surface_state use MOM_tracer_flow_control, only : tracer_flow_control_end use MOM_transcribe_grid, only : copy_dyngrid_to_MOM_grid, copy_MOM_grid_to_dyngrid +use MOM_transcribe_grid, only : rotate_dyngrid use MOM_unit_scaling, only : unit_scale_type, unit_scaling_init use MOM_unit_scaling, only : unit_scaling_end, fix_restart_unit_scaling use MOM_variables, only : surface, allocate_surface_state, deallocate_surface_state use MOM_variables, only : thermo_var_ptrs, vertvisc_type use MOM_variables, only : accel_diag_ptrs, cont_diag_ptrs, ocean_internal_state +use MOM_variables, only : rotate_surface_state use MOM_verticalGrid, only : verticalGrid_type, verticalGridInit, verticalGridEnd use MOM_verticalGrid, only : fix_restart_scaling use MOM_verticalGrid, only : get_thickness_units, get_flux_units, get_tr_flux_units @@ -180,7 +191,10 @@ module MOM real :: time_in_thermo_cycle !< The running time of the current time-stepping !! cycle in calls that step the thermodynamics [T ~> s]. - type(ocean_grid_type) :: G !< structure containing metrics and grid info + type(ocean_grid_type) :: G_in !< Input grid metric + type(ocean_grid_type), pointer :: G => NULL() !< Model grid metric + logical :: rotate_index = .false. !< True if index map is rotated + type(verticalGrid_type), pointer :: & GV => NULL() !< structure containing vertical grid info type(unit_scale_type), pointer :: & @@ -399,13 +413,13 @@ module MOM !! The action of lateral processes on tracers occur in calls to !! advect_tracer and tracer_hordiff. Vertical mixing and possibly remapping !! occur inside of diabatic. -subroutine step_MOM(forces, fluxes, sfc_state, Time_start, time_int_in, CS, & +subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS, & Waves, do_dynamics, do_thermodynamics, start_cycle, & end_cycle, cycle_length, reset_therm) - type(mech_forcing), intent(inout) :: forces !< A structure with the driving mechanical forces - type(forcing), intent(inout) :: fluxes !< A structure with pointers to themodynamic, + type(mech_forcing), target, intent(inout) :: forces_in !< A structure with the driving mechanical forces + type(forcing), target, intent(inout) :: fluxes_in !< A structure with pointers to themodynamic, !! tracer and mass exchange forcing fields - type(surface), intent(inout) :: sfc_state !< surface ocean state + type(surface), target, intent(inout) :: sfc_state !< surface ocean state type(time_type), intent(in) :: Time_start !< starting time of a segment, as a time type real, intent(in) :: time_int_in !< time interval covered by this run segment [s]. type(MOM_control_struct), pointer :: CS !< control structure from initialize_MOM @@ -430,6 +444,7 @@ subroutine step_MOM(forces, fluxes, sfc_state, Time_start, time_int_in, CS, & ! local variables type(ocean_grid_type), pointer :: G => NULL() ! pointer to a structure containing ! metrics and related information + type(ocean_grid_type), pointer :: G_in => NULL() ! Input grid metric type(verticalGrid_type), pointer :: GV => NULL() ! Pointer to the vertical grid structure type(unit_scale_type), pointer :: US => NULL() ! Pointer to a structure containing ! various unit conversion factors @@ -480,7 +495,13 @@ subroutine step_MOM(forces, fluxes, sfc_state, Time_start, time_int_in, CS, & type(group_pass_type) :: pass_tau_ustar_psurf logical :: showCallTree - G => CS%G ; GV => CS%GV ; US => CS%US + ! External forcing fields on the model index map + type(mech_forcing), pointer :: forces ! Mechanical forcing + type(forcing), pointer :: fluxes ! Boundary fluxes + type(surface), pointer :: sfc_state_diag ! Surface boundary fields + integer :: turns ! Number of quarter turns from input to model indexing + + G => CS%G ; G_in => CS%G_in ; GV => CS%GV ; US => CS%US is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed @@ -507,6 +528,21 @@ subroutine step_MOM(forces, fluxes, sfc_state, Time_start, time_int_in, CS, & showCallTree = callTree_showQuery() if (showCallTree) call callTree_enter("step_MOM(), MOM.F90") + ! Rotate the forces from G_in to G + if (CS%rotate_index) then + turns = G%HI%turns + allocate(forces) + call allocate_mech_forcing(forces_in, G, forces) + call rotate_mech_forcing(forces_in, turns, forces) + + allocate(fluxes) + call allocate_forcing_type(fluxes_in, G, fluxes) + call rotate_forcing(fluxes_in, fluxes, turns) + else + forces => forces_in + fluxes => fluxes_in + endif + ! First determine the time step that is consistent with this call and an ! integer fraction of time_interval. if (do_dyn) then @@ -838,19 +874,27 @@ subroutine step_MOM(forces, fluxes, sfc_state, Time_start, time_int_in, CS, & endif if (showCallTree) call callTree_waypoint("calling extract_surface_state (step_MOM)") + ! NOTE: sfc_state uses input indexing, since it is also used by drivers. call extract_surface_state(CS, sfc_state) ! Do diagnostics that only occur at the end of a complete forcing step. if (cycle_end) then + if (CS%rotate_index) then + allocate(sfc_state_diag) + call rotate_surface_state(sfc_state, G_in, sfc_state_diag, G, turns) + else + sfc_state_diag => sfc_state + endif + call cpu_clock_begin(id_clock_diagnostics) if (CS%time_in_cycle > 0.0) then call enable_averages(CS%time_in_cycle, Time_local, CS%diag) - call post_surface_dyn_diags(CS%sfc_IDs, G, CS%diag, sfc_state, ssh) + call post_surface_dyn_diags(CS%sfc_IDs, G, CS%diag, sfc_state_diag, ssh) endif if (CS%time_in_thermo_cycle > 0.0) then call enable_averages(CS%time_in_thermo_cycle, Time_local, CS%diag) call post_surface_thermo_diags(CS%sfc_IDs, G, GV, US, CS%diag, CS%time_in_thermo_cycle, & - sfc_state, CS%tv, ssh, CS%ave_ssh_ibc) + sfc_state_diag, CS%tv, ssh, CS%ave_ssh_ibc) endif call disable_averaging(CS%diag) call cpu_clock_end(id_clock_diagnostics) @@ -868,6 +912,17 @@ subroutine step_MOM(forces, fluxes, sfc_state, Time_start, time_int_in, CS, & call cpu_clock_end(id_clock_other) + ! De-rotate fluxes and copy back to the input, since they can be changed. + if (CS%rotate_index) then + call rotate_forcing(fluxes, fluxes_in, -turns) + + call deallocate_mech_forcing(forces) + deallocate(forces) + + call deallocate_forcing_type(fluxes) + deallocate(fluxes) + endif + if (showCallTree) call callTree_leave("step_MOM()") call cpu_clock_end(id_clock_ocean) @@ -1531,13 +1586,24 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & !! calls to step_MOM instead of the number of !! dynamics timesteps. ! local variables - type(ocean_grid_type), pointer :: G => NULL() ! A pointer to a structure with metrics and related - type(hor_index_type) :: HI ! A hor_index_type for array extents + type(ocean_grid_type), pointer :: G => NULL() ! A pointer to the metric grid use for the run + type(ocean_grid_type), pointer :: G_in => NULL() ! Pointer to the input grid + type(hor_index_type), pointer :: HI => NULL() ! A hor_index_type for array extents + type(hor_index_type), target :: HI_in ! HI on the input grid type(verticalGrid_type), pointer :: GV => NULL() type(dyn_horgrid_type), pointer :: dG => NULL() + type(dyn_horgrid_type), pointer :: dG_in => NULL() type(diag_ctrl), pointer :: diag => NULL() type(unit_scale_type), pointer :: US => NULL() character(len=4), parameter :: vers_num = 'v2.0' + integer :: turns ! Number of grid quarter-turns + + ! Initial state on the input index map + real, allocatable, dimension(:,:,:) :: u_in, v_in, h_in + real, allocatable, dimension(:,:,:), target :: T_in, S_in + type(ocean_OBC_type), pointer :: OBC_in => NULL() + type(sponge_CS), pointer :: sponge_in_CSp => NULL() + type(ALE_sponge_CS), pointer :: ALE_sponge_in_CSp => NULL() ! This include declares and sets the variable "version". # include "version_variable.h" @@ -1608,9 +1674,6 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & endif allocate(CS) - if (test_grid_copy) then ; allocate(G) - else ; G => CS%G ; endif - CS%Time => Time id_clock_init = cpu_clock_id('Ocean Initialization', grain=CLOCK_SUBCOMPONENT) @@ -1949,31 +2012,95 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & call callTree_waypoint("MOM parameters read (initialize_MOM)") + ! Grid rotation test + call get_param(param_file, "MOM", "ROTATE_INDEX", CS%rotate_index, & + "Enable rotation of the horizontal indices.", default=.false.) + if (CS%rotate_index) then + ! TODO: Index rotation currently only works when index rotation does not + ! change the MPI rank of each domain. Resolving this will require a + ! modification to FMS PE assignment. + ! For now, we only permit single-core runs. + + if (num_PEs() /= 1) & + call MOM_error(FATAL, "Index rotation is only supported on one PE.") + + call get_param(param_file, "MOM", "INDEX_TURNS", turns, & + "Number of counterclockwise quarter-turn index rotations.", default=1) + endif + ! Set up the model domain and grids. #ifdef SYMMETRIC_MEMORY_ symmetric = .true. #else symmetric = .false. #endif + G_in => CS%G_in #ifdef STATIC_MEMORY_ - call MOM_domains_init(G%domain, param_file, symmetric=symmetric, & + call MOM_domains_init(G_in%domain, param_file, symmetric=symmetric, & static_memory=.true., NIHALO=NIHALO_, NJHALO=NJHALO_, & NIGLOBAL=NIGLOBAL_, NJGLOBAL=NJGLOBAL_, NIPROC=NIPROC_, & NJPROC=NJPROC_) #else - call MOM_domains_init(G%domain, param_file, symmetric=symmetric) + call MOM_domains_init(G_in%domain, param_file, symmetric=symmetric, & + domain_name="MOM_in") #endif + + ! Copy input grid (G_in) domain to active grid G + ! Swap axes for quarter and 3-quarter turns + if (CS%rotate_index) then + allocate(CS%G) + call clone_MOM_domain(G_in%Domain, CS%G%Domain, turns=turns) + first_direction = modulo(first_direction + turns, 2) + else + CS%G => G_in + endif + + ! TODO: It is unlikey that test_grid_copy and rotate_index would work at the + ! same time. It may be possible to enable both but for now we prevent it. + if (test_grid_copy .and. CS%rotate_index) & + call MOM_error(FATAL, "Grid cannot be copied during index rotation.") + + if (test_grid_copy) then ; allocate(G) + else ; G => CS%G ; endif + call callTree_waypoint("domains initialized (initialize_MOM)") call MOM_debugging_init(param_file) call diag_mediator_infrastructure_init() call MOM_io_init(param_file) - call hor_index_init(G%Domain, HI, param_file, & + ! Create HI and dG on the input index map. + call hor_index_init(G_in%Domain, HI_in, param_file, & local_indexing=.not.global_indexing) + call create_dyn_horgrid(dG_in, HI_in, bathymetry_at_vel=bathy_at_vel) + call clone_MOM_domain(G_in%Domain, dG_in%Domain) + + ! Allocate initialize time-invariant MOM variables. + call MOM_initialize_fixed(dG_in, US, OBC_in, param_file, write_geom_files, & + dirs%output_directory) + + call callTree_waypoint("returned from MOM_initialize_fixed() (initialize_MOM)") - call create_dyn_horgrid(dG, HI, bathymetry_at_vel=bathy_at_vel) - call clone_MOM_domain(G%Domain, dG%Domain) + ! Determine HI and dG for the model index map. + if (CS%rotate_index) then + allocate(HI) + call rotate_hor_index(HI_in, turns, HI) + call create_dyn_horgrid(dG, HI, bathymetry_at_vel=bathy_at_vel) + call clone_MOM_domain(G%Domain, dG%Domain) + call rotate_dyngrid(dG_in, dG, US, turns) + if (associated(OBC_in)) then + ! TODO: General OBC index rotations is not yet supported. + if (modulo(turns, 4) /= 1) & + call MOM_error(FATAL, "OBC index rotation of 180 and 270 degrees is " & + // "not yet unsupported.") + allocate(CS%OBC) + call rotate_OBC_config(OBC_in, dG_in, CS%OBC, dG, turns) + endif + else + HI => HI_in + dG => dG_in + CS%OBC => OBC_in + endif call verticalGridInit( param_file, CS%GV, US ) GV => CS%GV @@ -1986,10 +2113,6 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & call MOM_timing_init(CS) - ! Allocate initialize time-invariant MOM variables. - call MOM_initialize_fixed(dG, US, CS%OBC, param_file, write_geom_files, dirs%output_directory) - call callTree_waypoint("returned from MOM_initialize_fixed() (initialize_MOM)") - if (associated(CS%OBC)) call call_OBC_register(param_file, CS%update_OBC_CSp, CS%OBC) call tracer_registry_init(param_file, CS%tracer_Reg) @@ -2045,6 +2168,18 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & flux_scale=conv2salt, convergence_units='kg m-2 s-1', & convergence_scale=0.001*GV%H_to_kg_m2, CMOR_tendprefix="osalt", diag_form=2) endif + ! NOTE: register_temp_salt_segments includes allocation of tracer fields + ! along segments. Bit reproducibility requires that MOM_initialize_state + ! be called on the input index map, so we must setup both OBC and OBC_in. + ! + ! XXX: This call on OBC_in allocates the tracer fields on the unrotated + ! grid, but also incorrectly stores a pointer to a tracer_type for the + ! rotated registry (e.g. segment%tr_reg%Tr(n)%Tr) from CS%tracer_reg. + ! + ! While incorrect and potentially dangerous, it does not seem that this + ! pointer is used during initialization, so we leave it for now. + if (CS%rotate_index .and. associated(OBC_in)) & + call register_temp_salt_segments(GV, OBC_in, CS%tracer_Reg, param_file) if (associated(CS%OBC)) & call register_temp_salt_segments(GV, CS%OBC, CS%tracer_Reg, param_file) endif @@ -2161,9 +2296,17 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & ! (potentially static) ocean-specific grid type. ! The next line would be needed if G%Domain had not already been init'd above: ! call clone_MOM_domain(dG%Domain, G%Domain) - call MOM_grid_init(G, param_file, US, HI, bathymetry_at_vel=bathy_at_vel) - call copy_dyngrid_to_MOM_grid(dG, G, US) - call destroy_dyn_horgrid(dG) + + ! NOTE: If indices are rotated, then G and G_in must both be initialized. + ! If not rotated, then G_in and G are the same grid. + if (CS%rotate_index) then + call MOM_grid_init(G, param_file, US, HI, bathymetry_at_vel=bathy_at_vel) + call copy_dyngrid_to_MOM_grid(dG, G, US) + call destroy_dyn_horgrid(dG) + endif + call MOM_grid_init(G_in, param_file, US, HI_in, bathymetry_at_vel=bathy_at_vel) + call copy_dyngrid_to_MOM_grid(dG_in, G_in, US) + call destroy_dyn_horgrid(dG_in) ! Set a few remaining fields that are specific to the ocean grid type. call set_first_direction(G, first_direction) @@ -2175,9 +2318,68 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & ! Consider removing this later? G%ke = GV%ke - 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%OBC, Time_in) + if (CS%rotate_index) then + G_in%ke = GV%ke + + allocate(u_in(G_in%IsdB:G_in%IedB, G_in%jsd:G_in%jed, nz)) + allocate(v_in(G_in%isd:G_in%ied, G_in%JsdB:G_in%JedB, nz)) + allocate(h_in(G_in%isd:G_in%ied, G_in%jsd:G_in%jed, nz)) + u_in(:,:,:) = 0.0 + v_in(:,:,:) = 0.0 + h_in(:,:,:) = GV%Angstrom_H + + if (use_temperature) then + allocate(T_in(G_in%isd:G_in%ied, G_in%jsd:G_in%jed, nz)) + allocate(S_in(G_in%isd:G_in%ied, G_in%jsd:G_in%jed, nz)) + T_in(:,:,:) = 0.0 + S_in(:,:,:) = 0.0 + + CS%tv%T => T_in + CS%tv%S => S_in + endif + + call MOM_initialize_state(u_in, v_in, h_in, CS%tv, Time, G_in, GV, US, & + param_file, dirs, restart_CSp, CS%ALE_CSp, CS%tracer_Reg, & + sponge_in_CSp, ALE_sponge_in_CSp, OBC_in, Time_in) + + if (use_temperature) then + CS%tv%T => CS%T + CS%tv%S => CS%S + endif + + call rotate_initial_state(u_in, v_in, h_in, T_in, S_in, use_temperature, & + turns, CS%u, CS%v, CS%h, CS%T, CS%S) + + if (associated(sponge_in_CSp)) then + ! TODO: Implementation and testing of non-ALE spong rotation + call MOM_error(FATAL, "Index rotation of non-ALE sponge is not yet " & + // "implemented.") + endif + + if (associated(ALE_sponge_in_CSp)) then + call rotate_ALE_sponge(ALE_sponge_in_CSp, G_in, CS%ALE_sponge_CSp, G, & + turns, param_file) + call update_ALE_sponge_field(CS%ALE_sponge_CSp, T_in, CS%T) + call update_ALE_sponge_field(CS%ALE_sponge_CSp, S_in, CS%S) + endif + + if (associated(OBC_in)) & + call rotate_OBC_init(OBC_in, G, GV, US, param_file, CS%tv, restart_CSp, & + CS%OBC) + + deallocate(u_in) + deallocate(v_in) + deallocate(h_in) + if (use_temperature) then + deallocate(T_in) + deallocate(S_in) + endif + 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%OBC, Time_in) + endif + call cpu_clock_end(id_clock_MOM_init) call callTree_waypoint("returned from MOM_initialize_state() (initialize_MOM)") @@ -2469,7 +2671,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & CS%nstep_tot = 0 if (present(count_calls)) CS%count_calls = count_calls - call MOM_sum_output_init(G, US, param_file, dirs%output_directory, & + call MOM_sum_output_init(G_in, US, param_file, dirs%output_directory, & CS%ntrunc, Time_init, CS%sum_output_CSp) ! Flag whether to save initial conditions in finish_MOM_initialization() or not. @@ -2526,7 +2728,7 @@ subroutine finish_MOM_initialization(Time, dirs, CS, restart_CSp) call register_restart_field(z_interface, "eta", .true., restart_CSp_tmp, & "Interface heights", "meter", z_grid='i') - call save_restart(dirs%output_directory, Time, G, & + call save_restart(dirs%output_directory, Time, CS%G_in, & restart_CSp_tmp, filename=CS%IC_file, GV=GV) deallocate(z_interface) deallocate(restart_CSp_tmp) @@ -2617,17 +2819,20 @@ subroutine set_restart_fields(GV, US, param_file, CS, restart_CSp) type(verticalGrid_type), intent(inout) :: GV !< ocean vertical grid structure type(unit_scale_type), intent(inout) :: US !< A dimensional unit scaling type type(param_file_type), intent(in) :: param_file !< opened file for parsing to get parameters - type(MOM_control_struct), intent(in) :: CS !< control structure set up by inialize_MOM + type(MOM_control_struct), intent(in) :: CS !< control structure set up by initialize_MOM type(MOM_restart_CS), pointer :: restart_CSp !< pointer to the restart control !! structure that will be used for MOM. ! Local variables logical :: use_ice_shelf ! Needed to determine whether to add CS%Hml to restarts character(len=48) :: thickness_units, flux_units - + type(vardesc) :: u_desc, v_desc thickness_units = get_thickness_units(GV) flux_units = get_flux_units(GV) + u_desc = var_desc("u", "m s-1", "Zonal velocity", hor_grid='Cu') + v_desc = var_desc("v", "m s-1", "Meridional velocity", hor_grid='Cv') + if (associated(CS%tv%T)) & call register_restart_field(CS%tv%T, "Temp", .true., restart_CSp, & "Potential Temperature", "degC") @@ -2638,11 +2843,7 @@ subroutine set_restart_fields(GV, US, param_file, CS, restart_CSp) call register_restart_field(CS%h, "h", .true., restart_CSp, & "Layer Thickness", thickness_units) - call register_restart_field(CS%u, "u", .true., restart_CSp, & - "Zonal velocity", "m s-1", hor_grid='Cu') - - call register_restart_field(CS%v, "v", .true., restart_CSp, & - "Meridional velocity", "m s-1", hor_grid='Cv') + call register_restart_pair(CS%u, CS%v, u_desc, v_desc, .true., restart_CSp) if (associated(CS%tv%frazil)) & call register_restart_field(CS%tv%frazil, "frazil", .false., restart_CSp, & @@ -2719,18 +2920,20 @@ end subroutine adjust_ssh_for_p_atm !> Set the surface (return) properties of the ocean model by !! setting the appropriate fields in sfc_state. Unused fields !! are set to NULL or are unallocated. -subroutine extract_surface_state(CS, sfc_state) - type(MOM_control_struct), pointer :: CS !< Master MOM control structure - type(surface), intent(inout) :: sfc_state !< transparent ocean surface state - !! structure shared with the calling routine - !! data in this structure is intent out. +subroutine extract_surface_state(CS, sfc_state_in) + type(MOM_control_struct), pointer :: CS !< Master MOM control structure + type(surface), target, intent(inout) :: sfc_state_in !< transparent ocean surface state + !! structure shared with the calling routine + !! data in this structure is intent out. ! Local variables real :: hu, hv ! Thicknesses interpolated to velocity points [H ~> m or kg m-2] - type(ocean_grid_type), pointer :: G => NULL() !< pointer to a structure containing - !! metrics and related information + type(ocean_grid_type), pointer :: G => NULL() !< pointer to a structure containing + !! metrics and related information + type(ocean_grid_type), pointer :: G_in => NULL() !< Input grid metric type(verticalGrid_type), pointer :: GV => NULL() !< structure containing vertical grid info type(unit_scale_type), pointer :: US => NULL() !< structure containing various unit conversion factors + type(surface), pointer :: sfc_state => NULL() ! surface state on the model grid real, dimension(:,:,:), pointer :: & h => NULL() !< h : layer thickness [H ~> m or kg m-2] real :: depth(SZI_(CS%G)) !< Distance from the surface in depth units [Z ~> m] or [H ~> m or kg m-2] @@ -2752,9 +2955,10 @@ subroutine extract_surface_state(CS, sfc_state) integer :: iscB, iecB, jscB, jecB, isdB, iedB, jsdB, jedB logical :: localError character(240) :: msg + integer :: turns ! Number of quarter turns call callTree_enter("extract_surface_state(), MOM.F90") - G => CS%G ; GV => CS%GV ; US => CS%US + G => CS%G ; G_in => CS%G_in ; GV => CS%GV ; US => CS%US 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 iscB = G%iscB ; iecB = G%iecB; jscB = G%jscB ; jecB = G%jecB @@ -2763,12 +2967,24 @@ subroutine extract_surface_state(CS, sfc_state) use_temperature = associated(CS%tv%T) - if (.not.sfc_state%arrays_allocated) then + turns = 0 + if (CS%rotate_index) & + turns = G%HI%turns + + if (.not.sfc_state_in%arrays_allocated) & ! Consider using a run-time flag to determine whether to do the vertical ! integrals, since the 3-d sums are not negligible in cost. - call allocate_surface_state(sfc_state, G, use_temperature, do_integrals=.true., & - omit_frazil=.not.associated(CS%tv%frazil)) + call allocate_surface_state(sfc_state_in, G_in, use_temperature, & + do_integrals=.true., omit_frazil=.not.associated(CS%tv%frazil)) + + if (CS%rotate_index) then + allocate(sfc_state) + call allocate_surface_state(sfc_state, G, use_temperature, & + do_integrals=.true., omit_frazil=.not.associated(CS%tv%frazil)) + else + sfc_state => sfc_state_in endif + sfc_state%T_is_conT = CS%tv%T_is_conT sfc_state%S_is_absS = CS%tv%S_is_absS @@ -3103,9 +3319,31 @@ subroutine extract_surface_state(CS, sfc_state) if (CS%debug) call MOM_surface_chksum("Post extract_sfc", sfc_state, G) + ! Rotate sfc_state back onto the input grid, sfc_state_in + if (CS%rotate_index) then + call rotate_surface_state(sfc_state, G, sfc_state_in, G_in, -turns) + call deallocate_surface_state(sfc_state) + endif + call callTree_leave("extract_surface_sfc_state()") end subroutine extract_surface_state +!> Rotate initialization fields from input to rotated arrays. +subroutine rotate_initial_state(u_in, v_in, h_in, T_in, S_in, & + use_temperature, turns, u, v, h, T, S) + real, dimension(:,:,:), intent(in) :: u_in, v_in, h_in, T_in, S_in + logical, intent(in) :: use_temperature + integer, intent(in) :: turns + real, dimension(:,:,:), intent(out) :: u, v, h, T, S + + call rotate_vector(u_in, v_in, turns, u, v) + call rotate_array(h_in, turns, h) + if (use_temperature) then + call rotate_array(T_in, turns, T) + call rotate_array(S_in, turns, S) + endif +end subroutine rotate_initial_state + !> Return true if all phases of step_MOM are at the same point in time. function MOM_state_is_synchronized(CS, adv_dyn) result(in_synch) type(MOM_control_struct), pointer :: CS !< MOM control structure @@ -3138,7 +3376,7 @@ subroutine get_MOM_state_elements(CS, G, GV, US, C_p, C_p_scaled, use_temp) !! units [Q degC-1 ~> J kg degC-1] logical, optional, intent(out) :: use_temp !< True if temperature is a state variable - if (present(G)) G => CS%G + if (present(G)) G => CS%G_in if (present(GV)) GV => CS%GV if (present(US)) US => CS%US if (present(C_p)) C_p = CS%US%Q_to_J_kg * CS%tv%C_p diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90 index 5998f08c16..4f7679b8e2 100644 --- a/src/core/MOM_barotropic.F90 +++ b/src/core/MOM_barotropic.F90 @@ -20,7 +20,8 @@ module MOM_barotropic use MOM_open_boundary, only : ocean_OBC_type, OBC_SIMPLE, OBC_NONE, open_boundary_query use MOM_open_boundary, only : OBC_DIRECTION_E, OBC_DIRECTION_W use MOM_open_boundary, only : OBC_DIRECTION_N, OBC_DIRECTION_S, OBC_segment_type -use MOM_restart, only : register_restart_field, query_initialized, MOM_restart_CS +use MOM_restart, only : register_restart_field, register_restart_pair +use MOM_restart, only : query_initialized, MOM_restart_CS use MOM_tidal_forcing, only : tidal_forcing_sensitivity, tidal_forcing_CS use MOM_time_manager, only : time_type, real_to_time, operator(+), operator(-) use MOM_unit_scaling, only : unit_scale_type @@ -1536,11 +1537,15 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, if (.not. use_BT_cont) then call uvchksum("BT Dat[uv]", Datu, Datv, CS%debug_BT_HI, haloshift=1, scale=US%L_to_m*GV%H_to_m) endif - call uvchksum("BT wt_[uv]", wt_u, wt_v, G%HI, 0, .true., .true.) - call uvchksum("BT frhat[uv]", CS%frhatu, CS%frhatv, G%HI, 0, .true., .true.) + call uvchksum("BT wt_[uv]", wt_u, wt_v, G%HI, haloshift=0, & + symmetric=.true., omit_corners=.true., scalar_pair=.true.) + call uvchksum("BT frhat[uv]", CS%frhatu, CS%frhatv, G%HI, haloshift=0, & + symmetric=.true., omit_corners=.true., scalar_pair=.true.) call uvchksum("BT bc_accel_[uv]", bc_accel_u, bc_accel_v, G%HI, haloshift=0, scale=US%L_T2_to_m_s2) - call uvchksum("BT IDat[uv]", CS%IDatu, CS%IDatv, G%HI, haloshift=0, scale=US%m_to_Z) - call uvchksum("BT visc_rem_[uv]", visc_rem_u, visc_rem_v, G%HI, haloshift=1) + call uvchksum("BT IDat[uv]", CS%IDatu, CS%IDatv, G%HI, haloshift=0, & + scale=US%m_to_Z, scalar_pair=.true.) + call uvchksum("BT visc_rem_[uv]", visc_rem_u, visc_rem_v, G%HI, & + haloshift=1, scalar_pair=.true.) endif if (query_averaging_enabled(CS%diag)) then @@ -3113,9 +3118,13 @@ subroutine btcalc(h, G, GV, CS, h_u, h_v, may_use_default, OBC) enddo ; endif if (CS%debug) then - call uvchksum("btcalc frhat[uv]", CS%frhatu, CS%frhatv, G%HI, 0, .true., .true.) + call uvchksum("btcalc frhat[uv]", CS%frhatu, CS%frhatv, G%HI, & + haloshift=0, symmetric=.true., omit_corners=.true., & + scalar_pair=.true.) if (present(h_u) .and. present(h_v)) & - call uvchksum("btcalc h_[uv]", h_u, h_v, G%HI, 0, .true., .true., scale=GV%H_to_m) + call uvchksum("btcalc h_[uv]", h_u, h_v, G%HI, haloshift=0, & + symmetric=.true., omit_corners=.true., scale=GV%H_to_m, & + scalar_pair=.true.) call hchksum(h, "btcalc h",G%HI, haloshift=1, scale=GV%H_to_m) endif @@ -4235,6 +4244,7 @@ subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS, CS%debug_BT_HI%IedB=CS%iedw CS%debug_BT_HI%JsdB=CS%jsdw-1 CS%debug_BT_HI%JedB=CS%jedw + CS%debug_BT_HI%turns = G%HI%turns endif ! IareaT, IdxCu, and IdyCv need to be allocated with wide halos. @@ -4607,6 +4617,7 @@ subroutine register_barotropic_restarts(HI, GV, param_file, CS, restart_CS) type(vardesc) :: vd(3) real :: slow_rate integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB + isd = HI%isd ; ied = HI%ied ; jsd = HI%jsd ; jed = HI%jed IsdB = HI%IsdB ; IedB = HI%IedB ; JsdB = HI%JsdB ; JedB = HI%JedB @@ -4628,8 +4639,7 @@ subroutine register_barotropic_restarts(HI, GV, param_file, CS, restart_CS) hor_grid='u', z_grid='1') vd(3) = var_desc("vbtav","m s-1","Time mean barotropic meridional velocity",& hor_grid='v', z_grid='1') - call register_restart_field(CS%ubtav, vd(2), .false., restart_CS) - call register_restart_field(CS%vbtav, vd(3), .false., restart_CS) + call register_restart_pair(CS%ubtav, CS%vbtav, vd(2), vd(3), .false., restart_CS) vd(2) = var_desc("ubt_IC", "m s-1", & longname="Next initial condition for the barotropic zonal velocity", & @@ -4637,8 +4647,7 @@ subroutine register_barotropic_restarts(HI, GV, param_file, CS, restart_CS) vd(3) = var_desc("vbt_IC", "m s-1", & longname="Next initial condition for the barotropic meridional velocity",& hor_grid='v', z_grid='1') - call register_restart_field(CS%ubt_IC, vd(2), .false., restart_CS) - call register_restart_field(CS%vbt_IC, vd(3), .false., restart_CS) + call register_restart_pair(CS%ubt_IC, CS%vbt_IC, vd(2), vd(3), .false., restart_CS) if (GV%Boussinesq) then vd(2) = var_desc("uhbt_IC", "m3 s-1", & @@ -4655,8 +4664,7 @@ subroutine register_barotropic_restarts(HI, GV, param_file, CS, restart_CS) longname="Next initial condition for the barotropic meridional transport",& hor_grid='v', z_grid='1') endif - call register_restart_field(CS%uhbt_IC, vd(2), .false., restart_CS) - call register_restart_field(CS%vhbt_IC, vd(3), .false., restart_CS) + call register_restart_pair(CS%uhbt_IC, CS%vhbt_IC, vd(2), vd(3), .false., restart_CS) call register_restart_field(CS%dtbt, "DTBT", .false., restart_CS, & longname="Barotropic timestep", units="seconds") diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index 8c0decd8c1..62cff69a8c 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -28,7 +28,8 @@ module MOM_dynamics_split_RK2 use MOM_file_parser, only : get_param, log_version, param_file_type use MOM_get_input, only : directories use MOM_io, only : MOM_io_init, vardesc, var_desc -use MOM_restart, only : register_restart_field, query_initialized, save_restart +use MOM_restart, only : register_restart_field, register_restart_pair +use MOM_restart, only : query_initialized, save_restart use MOM_restart, only : restart_init, is_new_run, MOM_restart_CS use MOM_time_manager, only : time_type, time_type_to_real, operator(+) use MOM_time_manager, only : operator(-), operator(>), operator(*), operator(/) @@ -886,11 +887,12 @@ subroutine register_restarts_dyn_split_RK2(HI, GV, param_file, CS, restart_CS, u real, dimension(SZI_(HI),SZJB_(HI),SZK_(GV)), & target, intent(inout) :: vh !< merid volume/mass transport [H L2 T-1 ~> m3 s-1 or kg s-1] - type(vardesc) :: vd + type(vardesc) :: vd(2) character(len=40) :: mdl = "MOM_dynamics_split_RK2" ! This module's name. character(len=48) :: thickness_units, flux_units integer :: isd, ied, jsd, jed, nz, IsdB, IedB, JsdB, JedB + isd = HI%isd ; ied = HI%ied ; jsd = HI%jsd ; jed = HI%jed ; nz = GV%ke IsdB = HI%IsdB ; IedB = HI%IedB ; JsdB = HI%JsdB ; JedB = HI%JedB @@ -918,32 +920,26 @@ subroutine register_restarts_dyn_split_RK2(HI, GV, param_file, CS, restart_CS, u flux_units = get_flux_units(GV) if (GV%Boussinesq) then - vd = var_desc("sfc",thickness_units,"Free surface Height",'h','1') + vd(1) = var_desc("sfc",thickness_units,"Free surface Height",'h','1') else - vd = var_desc("p_bot",thickness_units,"Bottom Pressure",'h','1') + vd(1) = var_desc("p_bot",thickness_units,"Bottom Pressure",'h','1') endif - call register_restart_field(CS%eta, vd, .false., restart_CS) - - vd = var_desc("u2","m s-1","Auxiliary Zonal velocity",'u','L') - call register_restart_field(CS%u_av, vd, .false., restart_CS) - - vd = var_desc("v2","m s-1","Auxiliary Meridional velocity",'v','L') - call register_restart_field(CS%v_av, vd, .false., restart_CS) - - vd = var_desc("h2",thickness_units,"Auxiliary Layer Thickness",'h','L') - call register_restart_field(CS%h_av, vd, .false., restart_CS) + call register_restart_field(CS%eta, vd(1), .false., restart_CS) - vd = var_desc("uh",flux_units,"Zonal thickness flux",'u','L') - call register_restart_field(uh, vd, .false., restart_CS) + vd(1) = var_desc("u2","m s-1","Auxiliary Zonal velocity",'u','L') + vd(2) = var_desc("v2","m s-1","Auxiliary Meridional velocity",'v','L') + call register_restart_pair(CS%u_av, CS%v_av, vd(1), vd(2), .false., restart_CS) - vd = var_desc("vh",flux_units,"Meridional thickness flux",'v','L') - call register_restart_field(vh, vd, .false., restart_CS) + vd(1) = var_desc("h2",thickness_units,"Auxiliary Layer Thickness",'h','L') + call register_restart_field(CS%h_av, vd(1), .false., restart_CS) - vd = var_desc("diffu","m s-2","Zonal horizontal viscous acceleration",'u','L') - call register_restart_field(CS%diffu, vd, .false., restart_CS) + vd(1) = var_desc("uh",flux_units,"Zonal thickness flux",'u','L') + vd(2) = var_desc("vh",flux_units,"Meridional thickness flux",'v','L') + call register_restart_pair(uh, vh, vd(1), vd(2), .false., restart_CS) - vd = var_desc("diffv","m s-2","Meridional horizontal viscous acceleration",'v','L') - call register_restart_field(CS%diffv, vd, .false., restart_CS) + vd(1) = var_desc("diffu","m s-2","Zonal horizontal viscous acceleration",'u','L') + vd(2) = var_desc("diffv","m s-2","Meridional horizontal viscous acceleration",'v','L') + call register_restart_pair(CS%diffu, CS%diffv, vd(1), vd(2), .false., restart_CS) call register_barotropic_restarts(HI, GV, param_file, CS%barotropic_CSp, & restart_CS) diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index b7260c2da6..0a624b93e6 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -3,6 +3,7 @@ module MOM_forcing_type ! This file is part of MOM6. See LICENSE.md for the license. +use MOM_array_transform, only : rotate_array, rotate_vector, rotate_array_pair use MOM_debugging, only : hchksum, uvchksum use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, CLOCK_ROUTINE use MOM_diag_mediator, only : post_data, register_diag_field, register_scalar_field @@ -35,6 +36,21 @@ module MOM_forcing_type public copy_common_forcing_fields, allocate_mech_forcing, deallocate_mech_forcing public set_derived_forcing_fields, copy_back_forcing_fields public set_net_mass_forcing, get_net_mass_forcing +public rotate_forcing, rotate_mech_forcing + +!> Allocate the fields of a (flux) forcing type, based on either a set of input +!! flags for each group of fields, or a pre-allocated reference forcing. +interface allocate_forcing_type + module procedure allocate_forcing_by_group + module procedure allocate_forcing_by_ref +end interface allocate_forcing_type + +!> Allocate the fields of a mechanical forcing type, based on either a set of +!! input flags for each group of fields, or a pre-allocated reference forcing. +interface allocate_mech_forcing + module procedure allocate_mech_forcing_by_group + module procedure allocate_mech_forcing_from_ref +end interface allocate_mech_forcing ! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional ! consistency testing. These are noted in comments with units like Z, H, L, and T, along with @@ -2201,8 +2217,8 @@ end subroutine copy_back_forcing_fields !> Offer mechanical forcing fields for diagnostics for those !! fields registered as part of register_forcing_type_diags. -subroutine mech_forcing_diags(forces, dt, G, time_end, diag, handles) - type(mech_forcing), intent(in) :: forces !< A structure with the driving mechanical forces +subroutine mech_forcing_diags(forces_in, dt, G, time_end, diag, handles) + type(mech_forcing), target, intent(in) :: forces_in !< mechanical forcing input fields real, intent(in) :: dt !< time step for the forcing [s] type(ocean_grid_type), intent(in) :: G !< grid type type(time_type), intent(in) :: time_end !< The end time of the diagnostic interval. @@ -2211,8 +2227,22 @@ subroutine mech_forcing_diags(forces, dt, G, time_end, diag, handles) integer :: i,j,is,ie,js,je + type(mech_forcing), pointer :: forces + integer :: turns + call cpu_clock_begin(handles%id_clock_forcing) + ! NOTE: post_data expects data to be on the input index map, so any rotations + ! must be undone before saving the output. + turns = diag%G%HI%turns + if (turns /= 0) then + allocate(forces) + call allocate_mech_forcing(forces_in, diag%G, forces) + call rotate_mech_forcing(forces_in, turns, forces) + else + forces => forces_in + endif + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec call enable_averaging(dt, time_end, diag) ! if (query_averaging_enabled(diag)) then @@ -2232,34 +2262,56 @@ subroutine mech_forcing_diags(forces, dt, G, time_end, diag, handles) ! endif call disable_averaging(diag) + + if (turns /= 0) then + call deallocate_mech_forcing(forces) + deallocate(forces) + endif + call cpu_clock_end(handles%id_clock_forcing) end subroutine mech_forcing_diags !> Offer buoyancy forcing fields for diagnostics for those !! fields registered as part of register_forcing_type_diags. -subroutine forcing_diagnostics(fluxes, sfc_state, G, US, time_end, diag, handles) - type(forcing), intent(in) :: fluxes !< A structure containing thermodynamic forcing fields +subroutine forcing_diagnostics(fluxes_in, sfc_state, G_in, US, time_end, diag, handles) + type(forcing), target, intent(in) :: fluxes_in !< A structure containing thermodynamic forcing fields type(surface), intent(in) :: sfc_state !< A structure containing fields that !! describe the surface state of the ocean. - type(ocean_grid_type), intent(in) :: G !< grid type + type(ocean_grid_type), target, intent(in) :: G_in !< Input grid type type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(time_type), intent(in) :: time_end !< The end time of the diagnostic interval. type(diag_ctrl), intent(inout) :: diag !< diagnostic regulator type(forcing_diags), intent(inout) :: handles !< diagnostic ids ! local - real, dimension(SZI_(G),SZJ_(G)) :: res + type(ocean_grid_type), pointer :: G ! Grid metric on model index map + type(forcing), pointer :: fluxes ! Fluxes on the model index map + real, dimension(SZI_(diag%G),SZJ_(diag%G)) :: res real :: total_transport ! for diagnosing integrated boundary transport real :: ave_flux ! for diagnosing averaged boundary flux real :: C_p ! seawater heat capacity [J degC-1 kg-1] real :: RZ_T_conversion ! A combination of scaling factors for mass fluxes [kg T m-2 s-1 R-1 Z-1 ~> 1] real :: I_dt ! inverse time step [s-1] real :: ppt2mks ! conversion between ppt and mks + integer :: turns ! Number of index quarter turns integer :: i,j,is,ie,js,je call cpu_clock_begin(handles%id_clock_forcing) + ! NOTE: post_data expects data to be on the input index map, so any rotations + ! must be undone before saving the output. + turns = diag%G%HI%turns + if (turns /= 0) then + G => diag%G + allocate(fluxes) + call allocate_forcing_type(fluxes_in, G, fluxes) + call rotate_forcing(fluxes_in, fluxes, turns) + else + G => G_in + fluxes => fluxes_in + endif + C_p = US%Q_to_J_kg*fluxes%C_p RZ_T_conversion = US%RZ_T_to_kg_m2s I_dt = 1.0 / (US%T_to_s*fluxes%dt_buoy_accum) @@ -2806,12 +2858,18 @@ subroutine forcing_diagnostics(fluxes, sfc_state, G, US, time_end, diag, handles ! endif ! query_averaging_enabled call disable_averaging(diag) + if (turns /= 0) then + call deallocate_forcing_type(fluxes) + deallocate(fluxes) + endif + call cpu_clock_end(handles%id_clock_forcing) end subroutine forcing_diagnostics !> Conditionally allocate fields within the forcing type -subroutine allocate_forcing_type(G, fluxes, water, heat, ustar, press, shelf, iceberg, salt, fix_accum_bug) +subroutine allocate_forcing_by_group(G, fluxes, water, heat, ustar, press, & + shelf, iceberg, salt, fix_accum_bug) type(ocean_grid_type), intent(in) :: G !< Ocean grid structure type(forcing), intent(inout) :: fluxes !< A structure containing thermodynamic forcing fields logical, optional, intent(in) :: water !< If present and true, allocate water fluxes @@ -2879,11 +2937,61 @@ subroutine allocate_forcing_type(G, fluxes, water, heat, ustar, press, shelf, ic call myAlloc(fluxes%mass_berg,isd,ied,jsd,jed, iceberg) if (present(fix_accum_bug)) fluxes%gustless_accum_bug = .not.fix_accum_bug +end subroutine allocate_forcing_by_group + -end subroutine allocate_forcing_type +subroutine allocate_forcing_by_ref(fluxes_ref, G, fluxes) + type(forcing), intent(in) :: fluxes_ref !< Reference fluxes + type(ocean_grid_type), intent(in) :: G !< Grid metric of target fluxes + type(forcing), intent(out) :: fluxes !< Target fluxes -!> Conditionally allocate fields within the mechanical forcing type -subroutine allocate_mech_forcing(G, forces, stress, ustar, shelf, press, iceberg) + logical :: do_ustar, do_water, do_heat, do_salt, do_press, do_shelf, & + do_iceberg, do_heat_added, do_buoy + + call get_forcing_groups(fluxes_ref, do_water, do_heat, do_ustar, do_press, & + do_shelf, do_iceberg, do_salt, do_heat_added, do_buoy) + + call allocate_forcing_type(G, fluxes, do_water, do_heat, do_ustar, & + do_press, do_shelf, do_iceberg, do_salt) + + ! The following fluxes would typically be allocated by the driver + call myAlloc(fluxes%sw_vis_dir, G%isd, G%ied, G%jsd, G%jed, & + associated(fluxes_ref%sw_vis_dir)) + call myAlloc(fluxes%sw_vis_dif, G%isd, G%ied, G%jsd, G%jed, & + associated(fluxes_ref%sw_vis_dif)) + call myAlloc(fluxes%sw_nir_dir, G%isd, G%ied, G%jsd, G%jed, & + associated(fluxes_ref%sw_nir_dir)) + call myAlloc(fluxes%sw_nir_dif, G%isd, G%ied, G%jsd, G%jed, & + associated(fluxes_ref%sw_nir_dif)) + + call myAlloc(fluxes%salt_flux_in, G%isd, G%ied, G%jsd, G%jed, & + associated(fluxes_ref%salt_flux_in)) + call myAlloc(fluxes%salt_flux_added, G%isd, G%ied, G%jsd, G%jed, & + associated(fluxes_ref%salt_flux_added)) + + call myAlloc(fluxes%p_surf_full, G%isd, G%ied, G%jsd, G%jed, & + associated(fluxes_ref%p_surf_full)) + + call myAlloc(fluxes%heat_added, G%isd, G%ied, G%jsd, G%jed, & + associated(fluxes_ref%heat_added)) + call myAlloc(fluxes%buoy, G%isd, G%ied, G%jsd, G%jed, & + associated(fluxes_ref%buoy)) + + call myAlloc(fluxes%TKE_tidal, G%isd, G%ied, G%jsd, G%jed, & + associated(fluxes_ref%TKE_tidal)) + call myAlloc(fluxes%ustar_tidal, G%isd, G%ied, G%jsd, G%jed, & + associated(fluxes_ref%ustar_tidal)) + + ! This flag would normally be set by a control flag in allocate_forcing_type. + ! Here we copy the flag from the reference forcing. + fluxes%gustless_accum_bug = fluxes_ref%gustless_accum_bug +end subroutine allocate_forcing_by_ref + + +!> Conditionally allocate fields within the mechanical forcing type using +!! control flags. +subroutine allocate_mech_forcing_by_group(G, forces, stress, ustar, shelf, & + press, iceberg) type(ocean_grid_type), intent(in) :: G !< Ocean grid structure type(mech_forcing), intent(inout) :: forces !< Forcing fields structure @@ -2917,8 +3025,82 @@ subroutine allocate_mech_forcing(G, forces, stress, ustar, shelf, press, iceberg !These fields should only on allocated when iceberg area is being passed through the coupler. call myAlloc(forces%area_berg,isd,ied,jsd,jed, iceberg) call myAlloc(forces%mass_berg,isd,ied,jsd,jed, iceberg) +end subroutine allocate_mech_forcing_by_group + + +!> Conditionally allocate fields within the mechanical forcing type based on a +!! reference forcing. +subroutine allocate_mech_forcing_from_ref(forces_ref, G, forces) + type(mech_forcing), intent(in) :: forces_ref !< Reference forcing fields + type(ocean_grid_type), intent(in) :: G !< Grid metric of target forcing + type(mech_forcing), intent(out) :: forces !< Mechanical forcing fields + + logical :: do_stress, do_ustar, do_shelf, do_press, do_iceberg + + ! Identify the active fields in the reference forcing + call get_mech_forcing_groups(forces_ref, do_stress, do_ustar, do_shelf, & + do_press, do_iceberg) + + call allocate_mech_forcing(G, forces, do_stress, do_ustar, do_shelf, & + do_press, do_iceberg) +end subroutine allocate_mech_forcing_from_ref + + +!> Return flags indicating which groups of forcings are allocated +subroutine get_forcing_groups(fluxes, water, heat, ustar, press, shelf, & + iceberg, salt, heat_added, buoy) + type(forcing), intent(in) :: fluxes !< Reference flux fields + logical, intent(out) :: water !< True if fluxes contains water-based fluxes + logical, intent(out) :: heat !< True if fluxes contains heat-based fluxes + logical, intent(out) :: ustar !< True if fluxes contains ustar fluxes + logical, intent(out) :: press !< True if fluxes contains surface pressure + logical, intent(out) :: shelf !< True if fluxes contains ice shelf fields + logical, intent(out) :: iceberg !< True if fluxes contains iceberg fluxes + logical, intent(out) :: salt !< True if fluxes contains salt flux + logical, intent(out) :: heat_added !< True if fluxes contains explicit heat + logical, intent(out) :: buoy !< True if fluxes contains buoyancy fluxes + + ! NOTE: heat, salt, heat_added, and buoy would typically depend on each other + ! to some degree. But since this would be enforced at the driver level, + ! we handle them here as independent flags. + + ustar = associated(fluxes%ustar) & + .and. associated(fluxes%ustar_gustless) + ! TODO: Check for all associated fields, but for now just check one as a marker + water = associated(fluxes%evap) + heat = associated(fluxes%seaice_melt_heat) + salt = associated(fluxes%salt_flux) + press = associated(fluxes%p_surf) + shelf = associated(fluxes%frac_shelf_h) + iceberg = associated(fluxes%ustar_berg) + heat_added = associated(fluxes%heat_added) + buoy = associated(fluxes%buoy) +end subroutine get_forcing_groups + + +!> Return flags indicating which groups of mechanical forcings are allocated +subroutine get_mech_forcing_groups(forces, stress, ustar, shelf, press, iceberg) + type(mech_forcing), intent(in) :: forces !< Reference forcing fields + logical, intent(out) :: stress !< True if forces contains wind stress fields + logical, intent(out) :: ustar !< True if forces contains ustar field + logical, intent(out) :: shelf !< True if forces contains ice shelf fields + logical, intent(out) :: press !< True if forces contains pressure fields + logical, intent(out) :: iceberg !< True if forces contains iceberg fields + + stress = associated(forces%taux) & + .and. associated(forces%tauy) + ustar = associated(forces%ustar) + shelf = associated(forces%rigidity_ice_u) & + .and. associated(forces%rigidity_ice_v) & + .and. associated(forces%frac_shelf_u) & + .and. associated(forces%frac_shelf_v) + press = associated(forces%p_surf) & + .and. associated(forces%p_surf_full) & + .and. associated(forces%net_mass_src) + iceberg = associated(forces%area_berg) & + .and. associated(forces%mass_berg) +end subroutine get_mech_forcing_groups -end subroutine allocate_mech_forcing !> Allocates and zeroes-out array. subroutine myAlloc(array, is, ie, js, je, flag) @@ -3006,6 +3188,181 @@ subroutine deallocate_mech_forcing(forces) end subroutine deallocate_mech_forcing +!< Rotate the fluxes by a set number of quarter turns +subroutine rotate_forcing(fluxes_in, fluxes, turns) + type(forcing), intent(in) :: fluxes_in !< Input forcing struct + type(forcing), intent(inout) :: fluxes !< Rotated forcing struct + integer, intent(in) :: turns !< Number of quarter turns + + logical :: do_ustar, do_water, do_heat, do_salt, do_press, do_shelf, & + do_iceberg, do_heat_added, do_buoy + + call get_forcing_groups(fluxes_in, do_water, do_heat, do_ustar, do_press, & + do_shelf, do_iceberg, do_salt, do_heat_added, do_buoy) + + if (do_ustar) then + call rotate_array(fluxes_in%ustar, turns, fluxes%ustar) + call rotate_array(fluxes_in%ustar_gustless, turns, fluxes%ustar_gustless) + endif + + if (do_water) then + call rotate_array(fluxes_in%evap, turns, fluxes%evap) + call rotate_array(fluxes_in%lprec, turns, fluxes%lprec) + call rotate_array(fluxes_in%fprec, turns, fluxes%fprec) + call rotate_array(fluxes_in%vprec, turns, fluxes%vprec) + call rotate_array(fluxes_in%lrunoff, turns, fluxes%lrunoff) + call rotate_array(fluxes_in%frunoff, turns, fluxes%frunoff) + call rotate_array(fluxes_in%seaice_melt, turns, fluxes%seaice_melt) + call rotate_array(fluxes_in%netMassOut, turns, fluxes%netMassOut) + call rotate_array(fluxes_in%netMassIn, turns, fluxes%netMassIn) + call rotate_array(fluxes_in%netSalt, turns, fluxes%netSalt) + endif + + if (do_heat) then + call rotate_array(fluxes_in%seaice_melt_heat, turns, fluxes%seaice_melt_heat) + call rotate_array(fluxes_in%sw, turns, fluxes%sw) + call rotate_array(fluxes_in%lw, turns, fluxes%lw) + call rotate_array(fluxes_in%latent, turns, fluxes%latent) + call rotate_array(fluxes_in%sens, turns, fluxes%sens) + call rotate_array(fluxes_in%latent_evap_diag, turns, fluxes%latent_evap_diag) + call rotate_array(fluxes_in%latent_fprec_diag, turns, fluxes%latent_fprec_diag) + call rotate_array(fluxes_in%latent_frunoff_diag, turns, fluxes%latent_frunoff_diag) + endif + + if (do_salt) then + call rotate_array(fluxes_in%salt_flux, turns, fluxes%salt_flux) + endif + + if (do_heat .and. do_water) then + call rotate_array(fluxes_in%heat_content_cond, turns, fluxes%heat_content_cond) + call rotate_array(fluxes_in%heat_content_icemelt, turns, fluxes%heat_content_icemelt) + call rotate_array(fluxes_in%heat_content_lprec, turns, fluxes%heat_content_lprec) + call rotate_array(fluxes_in%heat_content_fprec, turns, fluxes%heat_content_fprec) + call rotate_array(fluxes_in%heat_content_vprec, turns, fluxes%heat_content_vprec) + call rotate_array(fluxes_in%heat_content_lrunoff, turns, fluxes%heat_content_lrunoff) + call rotate_array(fluxes_in%heat_content_frunoff, turns, fluxes%heat_content_frunoff) + call rotate_array(fluxes_in%heat_content_massout, turns, fluxes%heat_content_massout) + call rotate_array(fluxes_in%heat_content_massin, turns, fluxes%heat_content_massin) + endif + + if (do_press) then + call rotate_array(fluxes_in%p_surf, turns, fluxes%p_surf) + endif + + if (do_shelf) then + call rotate_array(fluxes_in%frac_shelf_h, turns, fluxes%frac_shelf_h) + call rotate_array(fluxes_in%ustar_shelf, turns, fluxes%ustar_shelf) + call rotate_array(fluxes_in%iceshelf_melt, turns, fluxes%iceshelf_melt) + endif + + if (do_iceberg) then + call rotate_array(fluxes_in%ustar_berg, turns, fluxes%ustar_berg) + call rotate_array(fluxes_in%area_berg, turns, fluxes%area_berg) + call rotate_array(fluxes_in%iceshelf_melt, turns, fluxes%iceshelf_melt) + endif + + if (do_heat_added) then + call rotate_array(fluxes_in%heat_added, turns, fluxes%heat_added) + endif + + ! The following fields are handled by drivers rather than control flags. + if (associated(fluxes_in%sw_vis_dir)) & + call rotate_array(fluxes_in%sw_vis_dir, turns, fluxes%sw_vis_dir) + if (associated(fluxes_in%sw_vis_dif)) & + call rotate_array(fluxes_in%sw_vis_dif, turns, fluxes%sw_vis_dif) + if (associated(fluxes_in%sw_nir_dir)) & + call rotate_array(fluxes_in%sw_nir_dir, turns, fluxes%sw_nir_dir) + if (associated(fluxes_in%sw_nir_dif)) & + call rotate_array(fluxes_in%sw_nir_dif, turns, fluxes%sw_nir_dif) + + if (associated(fluxes_in%salt_flux_in)) & + call rotate_array(fluxes_in%salt_flux_in, turns, fluxes%salt_flux_in) + if (associated(fluxes_in%salt_flux_added)) & + call rotate_array(fluxes_in%salt_flux_added, turns, fluxes%salt_flux_added) + + if (associated(fluxes_in%p_surf_full)) & + call rotate_array(fluxes_in%p_surf_full, turns, fluxes%p_surf_full) + + if (associated(fluxes_in%buoy)) & + call rotate_array(fluxes_in%buoy, turns, fluxes%buoy) + + if (associated(fluxes_in%TKE_tidal)) & + call rotate_array(fluxes_in%TKE_tidal, turns, fluxes%TKE_tidal) + if (associated(fluxes_in%ustar_tidal)) & + call rotate_array(fluxes_in%ustar_tidal, turns, fluxes%ustar_tidal) + + ! TODO: tracer flux rotation + if (coupler_type_initialized(fluxes%tr_fluxes)) & + call MOM_error(FATAL, "Rotation of tracer BC fluxes not yet implemented.") + + ! Scalars and flags + fluxes%accumulate_p_surf = fluxes_in%accumulate_p_surf + + fluxes%vPrecGlobalAdj = fluxes_in%vPrecGlobalAdj + fluxes%saltFluxGlobalAdj = fluxes_in%saltFluxGlobalAdj + fluxes%netFWGlobalAdj = fluxes_in%netFWGlobalAdj + fluxes%vPrecGlobalScl = fluxes_in%vPrecGlobalScl + fluxes%saltFluxGlobalScl = fluxes_in%saltFluxGlobalScl + fluxes%netFWGlobalScl = fluxes_in%netFWGlobalScl + + fluxes%fluxes_used = fluxes_in%fluxes_used + fluxes%dt_buoy_accum = fluxes_in%dt_buoy_accum + fluxes%C_p = fluxes_in%C_p + ! NOTE: gustless_accum_bug is set during allocation + + fluxes%num_msg = fluxes_in%num_msg + fluxes%max_msg = fluxes_in%max_msg +end subroutine rotate_forcing + +!< Rotate the forcing fields from the input domain +subroutine rotate_mech_forcing(forces_in, turns, forces) + type(mech_forcing), intent(in) :: forces_in !< Forcing on the input domain + integer, intent(in) :: turns !< Number of quarter-turns + type(mech_forcing), intent(inout) :: forces !< Forcing on the rotated domain + + logical :: do_stress, do_ustar, do_shelf, do_press, do_iceberg + + call get_mech_forcing_groups(forces_in, do_stress, do_ustar, do_shelf, & + do_press, do_iceberg) + + if (do_stress) & + call rotate_vector(forces_in%taux, forces_in%tauy, turns, & + forces%taux, forces%tauy) + + if (do_ustar) & + call rotate_array(forces_in%ustar, turns, forces%ustar) + + if (do_shelf) then + call rotate_array_pair( & + forces_in%rigidity_ice_u, forces_in%rigidity_ice_v, turns, & + forces%rigidity_ice_u, forces%rigidity_ice_v & + ) + call rotate_array_pair( & + forces_in%frac_shelf_u, forces_in%frac_shelf_v, turns, & + forces%frac_shelf_u, forces%frac_shelf_v & + ) + endif + + if (do_press) then + ! NOTE: p_surf_SSH either points to p_surf or p_surf_full + call rotate_array(forces_in%p_surf, turns, forces%p_surf) + call rotate_array(forces_in%p_surf_full, turns, forces%p_surf_full) + call rotate_array(forces_in%net_mass_src, turns, forces%net_mass_src) + endif + + if (do_iceberg) then + call rotate_array(forces_in%area_berg, turns, forces%area_berg) + call rotate_array(forces_in%mass_berg, turns, forces%mass_berg) + endif + + ! Copy fields + forces%dt_force_accum = forces_in%dt_force_accum + forces%net_mass_src_set = forces_in%net_mass_src_set + forces%accumulate_p_surf = forces_in%accumulate_p_surf + forces%accumulate_rigidity = forces_in%accumulate_rigidity + forces%initialized = forces_in%initialized +end subroutine rotate_mech_forcing + !> \namespace mom_forcing_type !! !! \section section_fluxes Boundary fluxes diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index 3b1559ab81..28e50f4be5 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -3,6 +3,8 @@ module MOM_open_boundary ! This file is part of MOM6. See LICENSE.md for the license. +use MOM_array_transform, only : rotate_array, rotate_array_pair +use MOM_array_transform, only : allocate_rotated_array use MOM_coms, only : sum_across_PEs use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, CLOCK_ROUTINE use MOM_diag_mediator, only : diag_ctrl, time_type @@ -16,7 +18,8 @@ module MOM_open_boundary use MOM_io, only : EAST_FACE, NORTH_FACE use MOM_io, only : slasher, read_data, field_size, SINGLE_FILE use MOM_io, only : vardesc, query_vardesc, var_desc -use MOM_restart, only : register_restart_field, query_initialized, MOM_restart_CS +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 use MOM_time_manager, only : time_type, time_type_to_real, operator(-) @@ -57,6 +60,8 @@ module MOM_open_boundary public open_boundary_register_restarts public update_segment_tracer_reservoirs public update_OBC_ramp +public rotate_OBC_config +public rotate_OBC_init integer, parameter, public :: OBC_NONE = 0 !< Indicates the use of no open boundary integer, parameter, public :: OBC_SIMPLE = 1 !< Indicates the use of a simple inflow open boundary @@ -74,11 +79,11 @@ module MOM_open_boundary integer :: fid !< handle from FMS associated with segment data on disk integer :: fid_dz !< handle from FMS associated with segment thicknesses on disk character(len=8) :: name !< a name identifier for the segment data - real, pointer, dimension(:,:,:) :: buffer_src=>NULL() !< buffer for segment data located at cell faces + real, dimension(:,:,:), allocatable :: buffer_src !< buffer for segment data located at cell faces !! and on the original vertical grid integer :: nk_src !< Number of vertical levels in the source data - real, dimension(:,:,:), pointer :: dz_src=>NULL() !< vertical grid cell spacing of the incoming segment - !! data, set in [Z ~> m] then scaled to [H ~> m or kg m-2] + real, dimension(:,:,:), allocatable :: dz_src !< vertical grid cell spacing of the incoming segment + !! data, set in [Z ~> m] then scaled to [H ~> m or kg m-2] real, dimension(:,:,:), pointer :: buffer_dst=>NULL() !< buffer src data remapped to the target vertical grid real, dimension(:,:), pointer :: bt_vel=>NULL() !< barotropic velocity [L T-1 ~> m s-1] real :: value !< constant value if fid is equal to -1 @@ -836,53 +841,116 @@ subroutine setup_segment_indices(G, seg, Is_obc, Ie_obc, Js_obc, Je_obc) integer, intent(in) :: Js_obc !< Q-point global j-index of start of segment integer, intent(in) :: Je_obc !< Q-point global j-index of end of segment ! Local variables - integer :: Isg,Ieg,Jsg,Jeg + integer :: IsgB, IegB, JsgB, JegB + integer :: isg, ieg, jsg, jeg ! Isg, Ieg will be I*_obc in global space - if (Ie_obc Ie_obc) then + ! Northern boundary + isg = IsgB + 1 + jsg = JsgB + ieg = IegB + jeg = JegB + endif + + if (Is_obc < Ie_obc) then + ! Southern boundary + isg = IsgB + 1 + jsg = JsgB + 1 + ieg = IegB + jeg = JegB + 1 + endif + + if (Js_obc < Je_obc) then + ! Eastern boundary + isg = IsgB + jsg = JsgB + 1 + ieg = IegB + jeg = JegB + endif + + if (Js_obc > Je_obc) then + ! Western boundary + isg = IsgB + 1 + jsg = JsgB + 1 + ieg = IegB + 1 + jeg = JegB endif ! Global space I*_obc but sorted - seg%HI%IsgB = Isg ; seg%HI%IegB = Ieg - seg%HI%isg = Isg+1 ; seg%HI%ieg = Ieg - seg%HI%JsgB = Jsg ; seg%HI%JegB = Jeg - seg%HI%jsg = Jsg+1 ; seg%HI%Jeg = Jeg + seg%HI%IsgB = IsgB + seg%HI%JegB = JegB + seg%HI%IegB = IegB + seg%HI%JsgB = JsgB + + seg%HI%isg = isg + seg%HI%jsg = jsg + seg%HI%ieg = ieg + seg%HI%jeg = jeg ! Move into local index space - Isg = Isg - G%idg_offset - Jsg = Jsg - G%jdg_offset - Ieg = Ieg - G%idg_offset - Jeg = Jeg - G%jdg_offset + IsgB = IsgB - G%idg_offset + JsgB = JsgB - G%jdg_offset + IegB = IegB - G%idg_offset + JegB = JegB - G%jdg_offset + + isg = isg - G%idg_offset + jsg = jsg - G%jdg_offset + ieg = ieg - G%idg_offset + jeg = jeg - G%jdg_offset ! This is the i-extent of the segment on this PE. ! The values are nonsense if the segment is not on this PE. - seg%HI%IsdB = min( max(Isg, G%HI%IsdB), G%HI%IedB) - seg%HI%IedB = min( max(Ieg, G%HI%IsdB), G%HI%IedB) - seg%HI%isd = min( max(Isg+1, G%HI%isd), G%HI%ied) - seg%HI%ied = min( max(Ieg, G%HI%isd), G%HI%ied) - seg%HI%IscB = min( max(Isg, G%HI%IscB), G%HI%IecB) - seg%HI%IecB = min( max(Ieg, G%HI%IscB), G%HI%IecB) - seg%HI%isc = min( max(Isg+1, G%HI%isc), G%HI%iec) - seg%HI%iec = min( max(Ieg, G%HI%isc), G%HI%iec) + seg%HI%IsdB = min(max(IsgB, G%HI%IsdB), G%HI%IedB) + seg%HI%IedB = min(max(IegB, G%HI%IsdB), G%HI%IedB) + seg%HI%isd = min(max(isg, G%HI%isd), G%HI%ied) + seg%HI%ied = min(max(ieg, G%HI%isd), G%HI%ied) + seg%HI%IscB = min(max(IsgB, G%HI%IscB), G%HI%IecB) + seg%HI%IecB = min(max(IegB, G%HI%IscB), G%HI%IecB) + seg%HI%isc = min(max(isg, G%HI%isc), G%HI%iec) + seg%HI%iec = min(max(ieg, G%HI%isc), G%HI%iec) ! This is the j-extent of the segment on this PE. ! The values are nonsense if the segment is not on this PE. - seg%HI%JsdB = min( max(Jsg, G%HI%JsdB), G%HI%JedB) - seg%HI%JedB = min( max(Jeg, G%HI%JsdB), G%HI%JedB) - seg%HI%jsd = min( max(Jsg+1, G%HI%jsd), G%HI%jed) - seg%HI%jed = min( max(Jeg, G%HI%jsd), G%HI%jed) - seg%HI%JscB = min( max(Jsg, G%HI%JscB), G%HI%JecB) - seg%HI%JecB = min( max(Jeg, G%HI%JscB), G%HI%JecB) - seg%HI%jsc = min( max(Jsg+1, G%HI%jsc), G%HI%jec) - seg%HI%jec = min( max(Jeg, G%HI%jsc), G%HI%jec) + seg%HI%JsdB = min(max(JsgB, G%HI%JsdB), G%HI%JedB) + seg%HI%JedB = min(max(JegB, G%HI%JsdB), G%HI%JedB) + seg%HI%jsd = min(max(jsg, G%HI%jsd), G%HI%jed) + seg%HI%jed = min(max(jeg, G%HI%jsd), G%HI%jed) + seg%HI%JscB = min(max(JsgB, G%HI%JscB), G%HI%JecB) + seg%HI%JecB = min(max(JegB, G%HI%JscB), G%HI%JecB) + seg%HI%jsc = min(max(jsg, G%HI%jsc), G%HI%jec) + seg%HI%jec = min(max(jeg, G%HI%jsc), G%HI%jec) end subroutine setup_segment_indices @@ -1787,7 +1855,7 @@ end subroutine open_boundary_impose_land_mask !> Make sure the OBC tracer reservoirs are initialized. subroutine setup_OBC_tracer_reservoirs(G, OBC) - type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure + type(ocean_grid_type), intent(in) :: G !< Ocean grid structure type(ocean_OBC_type), pointer :: OBC !< Open boundary control structure ! Local variables type(OBC_segment_type), pointer :: segment => NULL() @@ -3453,23 +3521,28 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) character(len=200) :: filename, OBC_file, inputdir ! Strings for file/path type(OBC_segment_type), pointer :: segment => NULL() integer, dimension(4) :: siz + real, dimension(:,:,:), pointer :: tmp_buffer_in => NULL() ! Unrotated input integer :: ni_seg, nj_seg ! number of src gridpoints along the segments + integer :: ni_buf, nj_buf ! Number of filled values in tmp_buffer integer :: i2, j2 ! indices for referencing local domain array integer :: is_obc, ie_obc, js_obc, je_obc ! segment indices within local domain integer :: ishift, jshift ! offsets for staggered locations real, dimension(:,:), pointer :: seg_vel => NULL() ! pointer to segment velocity array real, dimension(:,:), pointer :: seg_trans => NULL() ! pointer to segment transport array - real, dimension(:,:,:), allocatable :: tmp_buffer + real, dimension(:,:,:), allocatable, target :: tmp_buffer real, dimension(:), allocatable :: h_stack integer :: is_obc2, js_obc2 real :: net_H_src, net_H_int, scl_fac real, pointer, dimension(:,:) :: normal_trans_bt=>NULL() ! barotropic transport + integer :: turns ! Number of index quarter turns is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB nz=G%ke + turns = G%HI%turns + if (.not. associated(OBC)) return do n = 1, OBC%number_of_segments @@ -3477,6 +3550,7 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) if (.not. segment%on_pe) cycle ! continue to next segment if not in computational domain + ! NOTE: These are in segment%HI, but defined slightly differently ni_seg = segment%ie_obc-segment%is_obc+1 nj_seg = segment%je_obc-segment%js_obc+1 is_obc = max(segment%is_obc,isd-1) @@ -3580,6 +3654,7 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) segment%field(m)%buffer_dst(:,:,:)=0.0 endif ! read source data interpolated to the current model time + ! NOTE: buffer is sized for vertex points, but may be used for faces if (siz(1)==1) then if (OBC%brushcutter_mode) then allocate(tmp_buffer(1,nj_seg*2-1,segment%field(m)%nk_src)) ! segment data is currrently on supergrid @@ -3594,7 +3669,44 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) endif endif - call time_interp_external(segment%field(m)%fid,Time, tmp_buffer) + ! TODO: Since we conditionally rotate a subset of tmp_buffer_in after + ! reading the value, it is currently not possible to use the rotated + ! implementation of time_interp_external. + ! For now, we must explicitly allocate and rotate this array. + if (turns /= 0) then + if (modulo(turns, 2) /= 0) then + allocate(tmp_buffer_in(size(tmp_buffer, 2), size(tmp_buffer, 1), size(tmp_buffer, 3))) + else + allocate(tmp_buffer_in(size(tmp_buffer, 1), size(tmp_buffer, 2), size(tmp_buffer, 3))) + endif + else + tmp_buffer_in => tmp_buffer + endif + + call time_interp_external(segment%field(m)%fid,Time, tmp_buffer_in) + ! NOTE: Rotation of face-points require that we skip the final value + if (turns /= 0) then + ! TODO: This is hardcoded for 90 degrees, and needs to be generalized. + if (segment%is_E_or_W & + .and. .not. (segment%field(m)%name == 'V' .or. segment%field(m)%name == 'DVDX')) then + nj_buf = size(tmp_buffer, 2) - 1 + call rotate_array(tmp_buffer_in(:nj_buf,:,:), turns, tmp_buffer(:,:nj_buf,:)) + elseif (segment%is_N_or_S & + .and. .not. (segment%field(m)%name == 'U' .or. segment%field(m)%name == 'DUDY')) then + ni_buf = size(tmp_buffer, 1) - 1 + call rotate_array(tmp_buffer_in(:,:ni_buf,:), turns, tmp_buffer(:ni_buf,:,:)) + else + call rotate_array(tmp_buffer_in, turns, tmp_buffer) + endif + + ! TODO: This is hardcoded for 90 degrees, and needs to be generalized. + if (segment%field(m)%name == 'U' & + .or. segment%field(m)%name == 'DVDX' & + .or. segment%field(m)%name == 'DUDY') then + tmp_buffer(:,:,:) = -tmp_buffer(:,:,:) + endif + endif + if (OBC%brushcutter_mode) then if (segment%is_E_or_W) then if (segment%field(m)%name == 'V' .or. segment%field(m)%name == 'DVDX') then @@ -3629,7 +3741,21 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) endif endif if (segment%field(m)%nk_src > 1) then - call time_interp_external(segment%field(m)%fid_dz,Time, tmp_buffer) + call time_interp_external(segment%field(m)%fid_dz,Time, tmp_buffer_in) + if (turns /= 0) then + ! TODO: This is hardcoded for 90 degrees, and needs to be generalized. + if (segment%is_E_or_W & + .and. .not. (segment%field(m)%name == 'V' .or. segment%field(m)%name == 'DVDX')) then + nj_buf = size(tmp_buffer, 2) - 1 + call rotate_array(tmp_buffer_in(:nj_buf,:,:), turns, tmp_buffer(:,:nj_buf,:)) + elseif (segment%is_N_or_S & + .and. .not. (segment%field(m)%name == 'U' .or. segment%field(m)%name == 'DUDY')) then + ni_buf = size(tmp_buffer, 1) - 1 + call rotate_array(tmp_buffer_in(:,:ni_buf,:), turns, tmp_buffer(:ni_buf,:,:)) + else + call rotate_array(tmp_buffer_in, turns, tmp_buffer) + endif + endif if (OBC%brushcutter_mode) then if (segment%is_E_or_W) then if (segment%field(m)%name == 'V' .or. segment%field(m)%name == 'DVDX') then @@ -3763,6 +3889,8 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) segment%field(m)%buffer_dst(:,:,1) = segment%field(m)%buffer_src(:,:,1) ! initialize remap destination buffer endif deallocate(tmp_buffer) + if (turns /= 0) & + deallocate(tmp_buffer_in) else ! fid <= 0 (Uniform value) if (.not. associated(segment%field(m)%buffer_dst)) then if (segment%is_E_or_W) then @@ -4214,7 +4342,7 @@ subroutine register_temp_salt_segments(GV, OBC, tr_Reg, param_file) end subroutine register_temp_salt_segments subroutine fill_temp_salt_segments(G, OBC, tv) - type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure + type(ocean_grid_type), intent(in) :: G !< Ocean grid structure type(ocean_OBC_type), pointer :: OBC !< Open boundary structure type(thermo_var_ptrs), intent(inout) :: tv !< Thermodynamics structure @@ -4268,6 +4396,7 @@ subroutine fill_temp_salt_segments(G, OBC, tv) segment%tr_Reg%Tr(1)%tres(:,:,:) = segment%tr_Reg%Tr(1)%t(:,:,:) segment%tr_Reg%Tr(2)%tres(:,:,:) = segment%tr_Reg%Tr(2)%t(:,:,:) enddo + call setup_OBC_tracer_reservoirs(G, OBC) end subroutine fill_temp_salt_segments @@ -4513,7 +4642,7 @@ subroutine open_boundary_register_restarts(HI, GV, OBC, Reg, param_file, restart type(MOM_restart_CS), pointer :: restart_CSp !< Restart structure, data intent(inout) logical, intent(in) :: use_temperature !< If true, T and S are used ! Local variables - type(vardesc) :: vd + type(vardesc) :: vd(2) integer :: m, n character(len=100) :: mesg type(OBC_segment_type), pointer :: segment=>NULL() @@ -4537,27 +4666,31 @@ subroutine open_boundary_register_restarts(HI, GV, OBC, Reg, param_file, restart ! so much memory and disk space. *** if (OBC%radiation_BCs_exist_globally) then allocate(OBC%rx_normal(HI%isdB:HI%iedB,HI%jsd:HI%jed,GV%ke)) - OBC%rx_normal(:,:,:) = 0.0 - vd = var_desc("rx_normal", "m s-1", "Normal Phase Speed for EW radiation OBCs", 'u', 'L') - call register_restart_field(OBC%rx_normal, vd, .false., restart_CSp) allocate(OBC%ry_normal(HI%isd:HI%ied,HI%jsdB:HI%jedB,GV%ke)) + OBC%rx_normal(:,:,:) = 0.0 OBC%ry_normal(:,:,:) = 0.0 - vd = var_desc("ry_normal", "m s-1", "Normal Phase Speed for NS radiation OBCs", 'v', 'L') - call register_restart_field(OBC%ry_normal, vd, .false., restart_CSp) + + vd(1) = var_desc("rx_normal", "m s-1", "Normal Phase Speed for EW radiation OBCs", 'u', 'L') + vd(2) = var_desc("ry_normal", "m s-1", "Normal Phase Speed for NS radiation OBCs", 'v', 'L') + call register_restart_pair(OBC%rx_normal, OBC%ry_normal, vd(1), vd(2), & + .false., restart_CSp) endif + if (OBC%oblique_BCs_exist_globally) then allocate(OBC%rx_oblique(HI%isdB:HI%iedB,HI%jsd:HI%jed,GV%ke)) - OBC%rx_oblique(:,:,:) = 0.0 - vd = var_desc("rx_oblique", "m2 s-2", "Radiation Speed Squared for EW oblique OBCs", 'u', 'L') - call register_restart_field(OBC%rx_oblique, vd, .false., restart_CSp) allocate(OBC%ry_oblique(HI%isd:HI%ied,HI%jsdB:HI%jedB,GV%ke)) + OBC%rx_oblique(:,:,:) = 0.0 OBC%ry_oblique(:,:,:) = 0.0 - vd = var_desc("ry_oblique", "m2 s-2", "Radiation Speed Squared for NS oblique OBCs", 'v', 'L') - call register_restart_field(OBC%ry_oblique, vd, .false., restart_CSp) + + vd(1) = var_desc("rx_oblique", "m2 s-2", "Radiation Speed Squared for EW oblique OBCs", 'u', 'L') + vd(2) = var_desc("ry_oblique", "m2 s-2", "Radiation Speed Squared for NS oblique OBCs", 'v', 'L') + call register_restart_pair(OBC%rx_oblique, OBC%ry_oblique, vd(1), vd(2), & + .false., restart_CSp) + allocate(OBC%cff_normal(HI%IsdB:HI%IedB,HI%jsdB:HI%jedB,GV%ke)) OBC%cff_normal(:,:,:) = 0.0 - vd = var_desc("cff_normal", "m2 s-2", "denominator for oblique OBCs", 'q', 'L') - call register_restart_field(OBC%cff_normal, vd, .false., restart_CSp) + vd(1) = var_desc("cff_normal", "m2 s-2", "denominator for oblique OBCs", 'q', 'L') + call register_restart_field(OBC%cff_normal, vd(1), .false., restart_CSp) endif if (Reg%ntr == 0) return @@ -4583,9 +4716,15 @@ subroutine open_boundary_register_restarts(HI, GV, OBC, Reg, param_file, restart OBC%tres_x(:,:,:,:) = 0.0 do m=1,OBC%ntr if (OBC%tracer_x_reservoirs_used(m)) then - write(mesg,'("tres_x_",I3.3)') m - vd = var_desc(mesg,"Conc", "Tracer concentration for EW OBCs",'u','L') - call register_restart_field(OBC%tres_x(:,:,:,m), vd, .false., restart_CSp) + if (modulo(HI%turns, 2) /= 0) then + write(mesg,'("tres_y_",I3.3)') m + vd(1) = var_desc(mesg,"Conc", "Tracer concentration for NS OBCs",'v','L') + call register_restart_field(OBC%tres_x(:,:,:,m), vd(1), .false., restart_CSp) + else + write(mesg,'("tres_x_",I3.3)') m + vd(1) = var_desc(mesg,"Conc", "Tracer concentration for EW OBCs",'u','L') + call register_restart_field(OBC%tres_x(:,:,:,m), vd(1), .false., restart_CSp) + endif endif enddo endif @@ -4594,13 +4733,18 @@ subroutine open_boundary_register_restarts(HI, GV, OBC, Reg, param_file, restart OBC%tres_y(:,:,:,:) = 0.0 do m=1,OBC%ntr if (OBC%tracer_y_reservoirs_used(m)) then - write(mesg,'("tres_y_",I3.3)') m - vd = var_desc(mesg,"Conc", "Tracer concentration for NS OBCs",'v','L') - call register_restart_field(OBC%tres_y(:,:,:,m), vd, .false., restart_CSp) + if (modulo(HI%turns, 2) /= 0) then + write(mesg,'("tres_x_",I3.3)') m + vd(1) = var_desc(mesg,"Conc", "Tracer concentration for EW OBCs",'u','L') + call register_restart_field(OBC%tres_y(:,:,:,m), vd(1), .false., restart_CSp) + else + write(mesg,'("tres_y_",I3.3)') m + vd(1) = var_desc(mesg,"Conc", "Tracer concentration for NS OBCs",'v','L') + call register_restart_field(OBC%tres_y(:,:,:,m), vd(1), .false., restart_CSp) + endif endif enddo endif - end subroutine open_boundary_register_restarts !> Update the OBC tracer reservoirs after the tracers have been updated. @@ -4783,6 +4927,309 @@ subroutine adjustSegmentEtaToFitBathymetry(G, GV, US, segment,fld) end subroutine adjustSegmentEtaToFitBathymetry +!> This is more of a rotate initialization than an actual rotate +subroutine rotate_OBC_config(OBC_in, G_in, OBC, G, turns) + type(ocean_OBC_type), pointer, intent(in) :: OBC_in !< Input OBC + type(dyn_horgrid_type), intent(in) :: G_in !< Input grid metric + type(ocean_OBC_type), pointer, intent(inout) :: OBC !< Rotated OBC + type(dyn_horgrid_type), intent(in) :: G !< Rotated grid metric + integer, intent(in) :: turns !< Number of quarter turns + + integer :: l + + ! Scalar and logical transfer + OBC%number_of_segments = OBC_in%number_of_segments + OBC%g_Earth = OBC_in%g_Earth + OBC%ke = OBC_in%ke + OBC%user_BCs_set_globally = OBC_in%user_BCs_set_globally + + ! These are conditionally read and set if number_of_segments > 0 + OBC%zero_vorticity = OBC_in%zero_vorticity + OBC%freeslip_vorticity = OBC_in%freeslip_vorticity + OBC%computed_vorticity = OBC_in%computed_vorticity + OBC%specified_vorticity = OBC_in%specified_vorticity + OBC%zero_strain = OBC_in%zero_strain + OBC%freeslip_strain = OBC_in%freeslip_strain + OBC%computed_strain = OBC_in%computed_strain + OBC%specified_strain = OBC_in%specified_strain + OBC%zero_biharmonic = OBC_in%zero_biharmonic + OBC%silly_h = OBC_in%silly_h + OBC%silly_u = OBC_in%silly_u + + ! Segment rotation + allocate(OBC%segment(0:OBC%number_of_segments)) + do l = 0, OBC%number_of_segments + call rotate_OBC_segment_config(OBC_in%segment(l), G_in, OBC%segment(l), G, turns) + ! Data up to setup_[uv]_point_obc is needed for allocate_obc_segment_data! + call allocate_OBC_segment_data(OBC, OBC%segment(l)) + call rotate_OBC_segment_data(OBC_in%segment(l), OBC%segment(l), turns) + enddo + + ! The horizontal segment map + allocate(OBC%segnum_u(G%IsdB:G%IedB,G%jsd:G%jed)) + allocate(OBC%segnum_v(G%isd:G%ied,G%JsdB:G%JedB)) + call rotate_array_pair(OBC_in%segnum_u, OBC_in%segnum_v, turns, & + OBC%segnum_u, OBC%segnum_v) + + ! These are conditionally enabled during segment configuration + OBC%open_u_BCs_exist_globally = OBC_in%open_v_BCs_exist_globally + OBC%open_v_BCs_exist_globally = OBC_in%open_u_BCs_exist_globally + OBC%Flather_u_BCs_exist_globally = OBC_in%Flather_v_BCs_exist_globally + OBC%Flather_v_BCs_exist_globally = OBC_in%Flather_u_BCs_exist_globally + OBC%oblique_BCs_exist_globally = OBC_in%oblique_BCs_exist_globally + OBC%nudged_u_BCs_exist_globally = OBC_in%nudged_v_BCs_exist_globally + OBC%nudged_v_BCs_exist_globally = OBC_in%nudged_u_BCs_exist_globally + OBC%specified_u_BCs_exist_globally= OBC_in%specified_v_BCs_exist_globally + OBC%specified_v_BCs_exist_globally= OBC_in%specified_u_BCs_exist_globally + OBC%radiation_BCs_exist_globally = OBC_in%radiation_BCs_exist_globally + + ! These are set by initialize_segment_data + OBC%brushcutter_mode = OBC_in%brushcutter_mode + OBC%update_OBC = OBC_in%update_OBC + OBC%needs_IO_for_data = OBC_in%needs_IO_for_data + + OBC%ntr = OBC_in%ntr + + OBC%gamma_uv = OBC_in%gamma_uv + 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. + allocate(OBC%remap_CS) + OBC%remap_CS = OBC_in%remap_CS + + ! TODO: The OBC registry seems to be a list of "registered" OBC types. + ! It does not appear to be used, so for now we skip this record. + !OBC%OBC_Reg => OBC_in%OBC_Reg +end subroutine rotate_OBC_config + +!> Rotate the OBC segment configuration data from the input to model index map. +subroutine rotate_OBC_segment_config(segment_in, G_in, segment, G, turns) + type(OBC_segment_type), intent(in) :: segment_in !< Input OBC segment + type(dyn_horgrid_type), intent(in) :: G_in !< Input grid metric + type(OBC_segment_type), intent(inout) :: segment !< Rotated OBC segment + type(dyn_horgrid_type), intent(in) :: G !< Rotated grid metric + integer, intent(in) :: turns !< Number of quarter turns + + ! Global segment indices + integer :: Is_obc_in, Ie_obc_in, Js_obc_in, Je_obc_in ! Input domain + integer :: Is_obc, Ie_obc, Js_obc, Je_obc ! Rotated domain + + ! NOTE: A "rotation" of the OBC segment string would allow us to use + ! setup_[uv]_point_obc to set up most of this. For now, we just copy/swap + ! flags and manually rotate the indices. + + ! This is set if the segment is in the local grid + segment%on_pe = segment_in%on_pe + + ! Transfer configuration flags + segment%Flather = segment_in%Flather + segment%radiation = segment_in%radiation + segment%radiation_tan = segment_in%radiation_tan + segment%radiation_grad = segment_in%radiation_grad + segment%oblique = segment_in%oblique + segment%oblique_tan = segment_in%oblique_tan + segment%oblique_grad = segment_in%oblique_grad + segment%nudged = segment_in%nudged + segment%nudged_tan = segment_in%nudged_tan + segment%nudged_grad = segment_in%nudged_grad + segment%specified = segment_in%specified + segment%specified_tan = segment_in%specified_tan + segment%specified_grad = segment_in%specified_grad + segment%open = segment_in%open + segment%gradient = segment_in%gradient + + ! 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%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%values_needed = segment_in%values_needed + + ! These are conditionally set if nudged + segment%Velocity_nudging_timescale_in = segment_in%Velocity_nudging_timescale_in + segment%Velocity_nudging_timescale_out= segment_in%Velocity_nudging_timescale_out + + ! Rotate segment indices + + ! Reverse engineer the input [IJ][se]_obc segment indices + ! NOTE: The values stored in the segment are always saved in ascending order, + ! e.g. (is < ie). In order to use setup_segment_indices, we reorder the + ! indices here to indicate face direction. + ! Segment indices are also indexed locally, so we remove the halo offset. + if (segment_in%direction == OBC_DIRECTION_N) then + Is_obc_in = segment_in%Ie_obc + G_in%idg_offset + Ie_obc_in = segment_in%Is_obc + G_in%idg_offset + else + Is_obc_in = segment_in%Is_obc + G_in%idg_offset + Ie_obc_in = segment_in%Ie_obc + G_in%idg_offset + endif + + if (segment_in%direction == OBC_DIRECTION_W) then + Js_obc_in = segment_in%Je_obc + G_in%jdg_offset + Je_obc_in = segment_in%Js_obc + G_in%jdg_offset + else + Js_obc_in = segment_in%Js_obc + G_in%jdg_offset + Je_obc_in = segment_in%Je_obc + G_in%jdg_offset + endif + + ! TODO: This is hardcoded for 90 degrees, and needs to be generalized. + Is_obc = G_in%jegB - Js_obc_in + Ie_obc = G_in%JegB - Je_obc_in + Js_obc = Is_obc_in + Je_obc = Ie_obc_in + + ! Orientation is based on the index ordering, [IJ][se]_obc are re-ordered + ! after the index is set. So we now need to restore the original order + + call setup_segment_indices(G, segment, Is_obc, Ie_obc, Js_obc, Je_obc) + + ! Re-order [IJ][se]_obc back to ascending, and remove the halo offset. + if (Is_obc > Ie_obc) then + segment%Is_obc = Ie_obc - G%idg_offset + segment%Ie_obc = Is_obc - G%idg_offset + else + segment%Is_obc = Is_obc - G%idg_offset + segment%Ie_obc = Ie_obc - G%idg_offset + endif + + if (Js_obc > Je_obc) then + segment%Js_obc = Je_obc - G%jdg_offset + segment%Je_obc = Js_obc - G%jdg_offset + else + segment%Js_obc = Js_obc - G%jdg_offset + segment%Je_obc = Je_obc - G%jdg_offset + endif + + ! Reconfigure the directional flags + ! TODO: This is hardcoded for 90 degrees, and needs to be generalized. + select case (segment_in%direction) + case (OBC_DIRECTION_N) + segment%direction = OBC_DIRECTION_W + segment%is_E_or_W_2 = segment_in%is_N_or_S + segment%is_E_or_W = segment_in%is_N_or_S .and. segment_in%on_pe + segment%is_N_or_S = .false. + case (OBC_DIRECTION_W) + segment%direction = OBC_DIRECTION_S + segment%is_N_or_S = segment_in%is_E_or_W + segment%is_E_or_W = .false. + segment%is_E_or_W_2 = .false. + case (OBC_DIRECTION_S) + segment%direction = OBC_DIRECTION_E + segment%is_E_or_W_2 = segment_in%is_N_or_S + segment%is_E_or_W = segment_in%is_N_or_S .and. segment_in%on_pe + segment%is_N_or_S = .false. + case (OBC_DIRECTION_E) + segment%direction = OBC_DIRECTION_N + segment%is_N_or_S = segment_in%is_E_or_W + segment%is_E_or_W = .false. + segment%is_E_or_W_2 = .false. + case (OBC_NONE) + segment%direction = OBC_NONE + end select + + ! These are conditionally set if Lscale_{in,out} are present + segment%Tr_InvLscale_in = segment_in%Tr_InvLscale_in + segment%Tr_InvLscale_out = segment_in%Tr_InvLscale_out +end subroutine rotate_OBC_segment_config + + +!> Initialize the segments and field-related data of a rotated OBC. +subroutine rotate_OBC_init(OBC_in, G, GV, US, param_file, tv, restart_CSp, OBC) + type(ocean_OBC_type), pointer, intent(in) :: OBC_in !< OBC on input map + 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 + type(param_file_type), intent(in) :: param_file !< Input parameters + type(thermo_var_ptrs), intent(inout) :: tv !< Tracer fields + type(MOM_restart_CS), pointer, intent(in) :: restart_CSp !< Restart CS + type(ocean_OBC_type), pointer, intent(inout) :: OBC !< Rotated OBC + + logical :: use_temperature + integer :: l + + call get_param(param_file, "MOM", "ENABLE_THERMODYNAMICS", use_temperature, & + "If true, Temperature and salinity are used as state "//& + "variables.", default=.true., do_not_log=.true.) + + do l = 0, OBC%number_of_segments + call rotate_OBC_segment_data(OBC_in%segment(l), OBC%segment(l), G%HI%turns) + enddo + + if (use_temperature) & + call fill_temp_salt_segments(G, OBC, tv) + + call open_boundary_init(G, GV, US, param_file, OBC, restart_CSp) +end subroutine rotate_OBC_init + + +!> 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 + type(OBC_segment_type), intent(inout) :: segment + integer, intent(in) :: turns + + integer :: n + integer :: is, ie, js, je, nk + integer :: num_fields + + + num_fields = segment_in%num_fields + allocate(segment%field(num_fields)) + + segment%num_fields = segment_in%num_fields + do n = 1, num_fields + segment%field(n)%fid = segment_in%field(n)%fid + segment%field(n)%fid_dz = segment_in%field(n)%fid_dz + + if (modulo(turns, 2) /= 0) then + select case (segment_in%field(n)%name) + case ('U') + segment%field(n)%name = 'V' + case ('V') + segment%field(n)%name = 'U' + 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 + + if (allocated(segment_in%field(n)%buffer_src)) then + call allocate_rotated_array(segment_in%field(n)%buffer_src, & + lbound(segment_in%field(n)%buffer_src), turns, & + segment%field(n)%buffer_src) + 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_array(segment_in%field(n)%dz_src, & + lbound(segment_in%field(n)%dz_src), turns, & + segment%field(n)%dz_src) + call rotate_array(segment_in%field(n)%dz_src, turns, & + segment%field(n)%dz_src) + endif + + segment%field(n)%buffer_dst => NULL() + segment%field(n)%bt_vel => NULL() + + segment%field(n)%value = segment_in%field(n)%value + enddo + + 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 + !> \namespace mom_open_boundary !! This module implements some aspects of internal open boundary !! conditions in MOM. diff --git a/src/core/MOM_transcribe_grid.F90 b/src/core/MOM_transcribe_grid.F90 index 045fc9261c..51d44c1041 100644 --- a/src/core/MOM_transcribe_grid.F90 +++ b/src/core/MOM_transcribe_grid.F90 @@ -4,6 +4,7 @@ module MOM_transcribe_grid ! This file is part of MOM6. See LICENSE.md for the license. +use MOM_array_transform, only: rotate_array, rotate_array_pair use MOM_domains, only : pass_var, pass_vector use MOM_domains, only : To_All, SCALAR_PAIR, CGRID_NE, AGRID, BGRID_NE, CORNER use MOM_dyn_horgrid, only : dyn_horgrid_type, set_derived_dyn_horgrid @@ -11,9 +12,10 @@ module MOM_transcribe_grid use MOM_grid, only : ocean_grid_type, set_derived_metrics use MOM_unit_scaling, only : unit_scale_type + implicit none ; private -public copy_dyngrid_to_MOM_grid, copy_MOM_grid_to_dyngrid +public copy_dyngrid_to_MOM_grid, copy_MOM_grid_to_dyngrid, rotate_dyngrid contains @@ -305,4 +307,92 @@ subroutine copy_MOM_grid_to_dyngrid(oG, dG, US) end subroutine copy_MOM_grid_to_dyngrid +subroutine rotate_dyngrid(G_in, G, US, turns) + type(dyn_horgrid_type), intent(in) :: G_in !< Common horizontal grid type + type(dyn_horgrid_type), intent(inout) :: G !< Ocean grid type + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + integer, intent(in) :: turns !< Number of quarter turns + + integer :: jsc, jec, jscB, jecB + integer :: qturn + + ! Center point + call rotate_array(G_in%geoLonT, turns, G%geoLonT) + call rotate_array(G_in%geoLatT, turns, G%geoLatT) + call rotate_array_pair(G_in%dxT, G_in%dyT, turns, G%dxT, G%dyT) + call rotate_array(G_in%areaT, turns, G%areaT) + call rotate_array(G_in%bathyT, turns, G%bathyT) + + call rotate_array_pair(G_in%df_dx, G_in%df_dy, turns, G%df_dx, G%df_dy) + call rotate_array(G_in%sin_rot, turns, G%sin_rot) + call rotate_array(G_in%cos_rot, turns, G%cos_rot) + call rotate_array(G_in%mask2dT, turns, G%mask2dT) + + ! Face point + call rotate_array_pair(G_in%geoLonCu, G_in%geoLonCv, turns, & + G%geoLonCu, G%geoLonCv) + call rotate_array_pair(G_in%geoLatCu, G_in%geoLatCv, turns, & + G%geoLatCu, G%geoLatCv) + call rotate_array_pair(G_in%dxCu, G_in%dyCv, turns, G%dxCu, G%dyCv) + call rotate_array_pair(G_in%dxCv, G_in%dyCu, turns, G%dxCv, G%dyCu) + call rotate_array_pair(G_in%dx_Cv, G_in%dy_Cu, turns, G%dx_Cv, G%dy_Cu) + + call rotate_array_pair(G_in%mask2dCu, G_in%mask2dCv, turns, & + G%mask2dCu, G%mask2dCv) + call rotate_array_pair(G_in%areaCu, G_in%areaCv, turns, & + G%areaCu, G%areaCv) + call rotate_array_pair(G_in%IareaCu, G_in%IareaCv, turns, & + G%IareaCu, G%IareaCv) + + ! Vertex point + call rotate_array(G_in%geoLonBu, turns, G%geoLonBu) + call rotate_array(G_in%geoLatBu, turns, G%geoLatBu) + call rotate_array_pair(G_in%dxBu, G_in%dyBu, turns, G%dxBu, G%dyBu) + call rotate_array(G_in%areaBu, turns, G%areaBu) + call rotate_array(G_in%CoriolisBu, turns, G%CoriolisBu) + call rotate_array(G_in%mask2dBu, turns, G%mask2dBu) + + ! Topographic + G%bathymetry_at_vel = G_in%bathymetry_at_vel + if (G%bathymetry_at_vel) then + call rotate_array_pair(G_in%Dblock_u, G_in%Dblock_v, turns, & + G%Dblock_u, G%Dblock_v) + call rotate_array_pair(G_in%Dopen_u, G_in%Dopen_v, turns, & + G%Dopen_u, G%Dopen_v) + endif + + ! Nominal grid axes + ! TODO: We should not assign lat values to the lon axis, and vice versa. + ! We temporarily copy lat <-> lon since several components still expect + ! lat and lon sizes to match the first and second dimension sizes. + ! But we ought to instead leave them unchanged and adjust the references to + ! these axes. + if (modulo(turns, 2) /= 0) then + G%gridLonT(:) = G_in%gridLatT(G_in%jeg:G_in%jsg:-1) + G%gridLatT(:) = G_in%gridLonT(:) + G%gridLonB(:) = G_in%gridLatB(G_in%jeg:(G_in%jsg-1):-1) + G%gridLatB(:) = G_in%gridLonB(:) + else + G%gridLonT(:) = G_in%gridLonT(:) + G%gridLatT(:) = G_in%gridLatT(:) + G%gridLonB(:) = G_in%gridLonB(:) + G%gridLatB(:) = G_in%gridLatB(:) + endif + + G%x_axis_units = G_in%y_axis_units + G%y_axis_units = G_in%x_axis_units + G%south_lat = G_in%south_lat + G%west_lon = G_in%west_lon + G%len_lat = G_in%len_lat + G%len_lon = G_in%len_lon + + ! Rotation-invariant fields + G%areaT_global = G_in%areaT_global + G%IareaT_global = G_in%IareaT_global + G%Rad_Earth = G_in%Rad_Earth + G%max_depth = G_in%max_depth + + call set_derived_dyn_horgrid(G, US) +end subroutine rotate_dyngrid + end module MOM_transcribe_grid diff --git a/src/core/MOM_variables.F90 b/src/core/MOM_variables.F90 index 09cbd14c60..8e7dad1e51 100644 --- a/src/core/MOM_variables.F90 +++ b/src/core/MOM_variables.F90 @@ -3,6 +3,7 @@ module MOM_variables ! This file is part of MOM6. See LICENSE.md for the license. +use MOM_array_transform, only : rotate_array, rotate_vector use MOM_domains, only : MOM_domain_type, get_domain_extent, group_pass_type use MOM_debugging, only : hchksum use MOM_error_handler, only : MOM_error, FATAL @@ -11,6 +12,7 @@ module MOM_variables use coupler_types_mod, only : coupler_1d_bc_type, coupler_2d_bc_type use coupler_types_mod, only : coupler_type_spawn, coupler_type_destructor +use coupler_types_mod, only : coupler_type_initialized implicit none ; private @@ -18,6 +20,7 @@ module MOM_variables public allocate_surface_state, deallocate_surface_state, MOM_thermovar_chksum public ocean_grid_type, alloc_BT_cont_type, dealloc_BT_cont_type +public rotate_surface_state ! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional ! consistency testing. These are noted in comments with units like Z, H, L, and T, along with @@ -395,6 +398,79 @@ subroutine deallocate_surface_state(sfc_state) end subroutine deallocate_surface_state +!> Rotate the surface state fields from the input to the model indices. +subroutine rotate_surface_state(sfc_state_in, G_in, sfc_state, G, turns) + type(surface), intent(in) :: sfc_state_in + type(ocean_grid_type), intent(in) :: G_in + type(surface), intent(inout) :: sfc_state + type(ocean_grid_type), intent(in) :: G + integer, intent(in) :: turns + + logical :: use_temperature, do_integrals, use_melt_potential, use_iceshelves + + ! NOTE: Many of these are weak tests, since only one is checked + use_temperature = allocated(sfc_state_in%SST) & + .and. allocated(sfc_state_in%SSS) + use_melt_potential = allocated(sfc_state_in%melt_potential) + do_integrals = allocated(sfc_state_in%ocean_mass) + use_iceshelves = allocated(sfc_state_in%taux_shelf) & + .and. allocated(sfc_state_in%tauy_shelf) + + if (.not. sfc_state%arrays_allocated) then + call allocate_surface_state(sfc_state, G, & + use_temperature=use_temperature, & + do_integrals=do_integrals, & + use_meltpot=use_melt_potential, & + use_iceshelves=use_iceshelves & + ) + sfc_state%arrays_allocated = .true. + endif + + if (use_temperature) then + call rotate_array(sfc_state_in%SST, turns, sfc_state%SST) + call rotate_array(sfc_state_in%SSS, turns, sfc_state%SSS) + else + call rotate_array(sfc_state_in%sfc_density, turns, sfc_state%sfc_density) + endif + + call rotate_array(sfc_state_in%Hml, turns, sfc_state%Hml) + call rotate_vector(sfc_state_in%u, sfc_state_in%v, turns, & + sfc_state%u, sfc_state%v) + call rotate_array(sfc_state_in%sea_lev, turns, sfc_state%sea_lev) + + if (use_melt_potential) then + call rotate_array(sfc_state_in%melt_potential, turns, sfc_state%melt_potential) + endif + + if (do_integrals) then + call rotate_array(sfc_state_in%ocean_mass, turns, sfc_state%ocean_mass) + if (use_temperature) then + call rotate_array(sfc_state_in%ocean_heat, turns, sfc_state%ocean_heat) + call rotate_array(sfc_state_in%ocean_salt, turns, sfc_state%ocean_salt) + call rotate_array(sfc_state_in%SSS, turns, sfc_state%TempxPmE) + call rotate_array(sfc_state_in%salt_deficit, turns, sfc_state%salt_deficit) + call rotate_array(sfc_state_in%internal_heat, turns, sfc_state%internal_heat) + endif + endif + + if (use_iceshelves) then + call rotate_vector(sfc_state_in%taux_shelf, sfc_state_in%tauy_shelf, turns, & + sfc_state%taux_shelf, sfc_state%tauy_shelf) + endif + + if (use_temperature .and. allocated(sfc_state_in%frazil)) & + call rotate_array(sfc_state_in%frazil, turns, sfc_state%frazil) + + ! Scalar transfers + sfc_state%T_is_conT = sfc_state_in%T_is_conT + sfc_state%S_is_absS = sfc_state_in%S_is_absS + + ! TODO: tracer field rotation + if (coupler_type_initialized(sfc_state_in%tr_fields)) & + call MOM_error(FATAL, "Rotation of surface state tracers is not yet " & + // "implemented.") +end subroutine rotate_surface_state + !> Allocates the arrays contained within a BT_cont_type and initializes them to 0. subroutine alloc_BT_cont_type(BT_cont, G, alloc_faces) type(BT_cont_type), pointer :: BT_cont !< The BT_cont_type whose elements will be allocated diff --git a/src/diagnostics/MOM_debugging.F90 b/src/diagnostics/MOM_debugging.F90 index d4d267d50d..611e6da2fc 100644 --- a/src/diagnostics/MOM_debugging.F90 +++ b/src/diagnostics/MOM_debugging.F90 @@ -8,7 +8,7 @@ module MOM_debugging ! This file is part of MOM6. See LICENSE.md for the license. -use MOM_checksums, only : hchksum, Bchksum, qchksum, uvchksum +use MOM_checksums, only : hchksum, Bchksum, qchksum, uvchksum, hchksum_pair use MOM_checksums, only : is_NaN, chksum, MOM_checksums_init use MOM_coms, only : PE_here, root_PE, num_PEs, sum_across_PEs use MOM_coms, only : min_across_PEs, max_across_PEs, reproducing_sum @@ -27,7 +27,7 @@ module MOM_debugging public :: check_column_integral, check_column_integrals ! These interfaces come from MOM_checksums. -public :: hchksum, Bchksum, qchksum, is_NaN, chksum, uvchksum +public :: hchksum, Bchksum, qchksum, is_NaN, chksum, uvchksum, hchksum_pair !> Check for consistency between the duplicated points of a C-grid vector interface check_redundant diff --git a/src/diagnostics/MOM_sum_output.F90 b/src/diagnostics/MOM_sum_output.F90 index 775cf39c22..de0f4eb8cd 100644 --- a/src/diagnostics/MOM_sum_output.F90 +++ b/src/diagnostics/MOM_sum_output.F90 @@ -376,6 +376,7 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, OBC, dt_ mass_EFP, & ! Extended fixed point sums of total mass, etc. salt_EFP, heat_EFP, salt_chg_EFP, heat_chg_EFP, mass_chg_EFP, & mass_anom_EFP, salt_anom_EFP, heat_anom_EFP + real :: CFL_Iarea ! Direction-based inverse area used in CFL test [L-2]. real :: CFL_trans ! A transport-based definition of the CFL number [nondim]. real :: CFL_lin ! A simpler definition of the CFL number [nondim]. real :: max_CFL(2) ! The maxima of the CFL numbers [nondim]. @@ -719,21 +720,21 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, OBC, dt_ ! Calculate the maximum CFL numbers. max_CFL(1:2) = 0.0 do k=1,nz ; do j=js,je ; do I=Isq,Ieq - if (u(I,j,k) < 0.0) then - CFL_trans = (-u(I,j,k) * CS%dt_in_T) * (G%dy_Cu(I,j) * G%IareaT(i+1,j)) - else - CFL_trans = (u(I,j,k) * CS%dt_in_T) * (G%dy_Cu(I,j) * G%IareaT(i,j)) - endif + CFL_Iarea = G%IareaT(i,j) + if (u(I,j,k) < 0.0) & + CFL_Iarea = G%IareaT(i+1,j) + + CFL_trans = abs(u(I,j,k) * CS%dt_in_T) * (G%dy_Cu(I,j) * CFL_Iarea) CFL_lin = abs(u(I,j,k) * CS%dt_in_T) * G%IdxCu(I,j) max_CFL(1) = max(max_CFL(1), CFL_trans) max_CFL(2) = max(max_CFL(2), CFL_lin) enddo ; enddo ; enddo do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie - if (v(i,J,k) < 0.0) then - CFL_trans = (-v(i,J,k) * CS%dt_in_T) * (G%dx_Cv(i,J) * G%IareaT(i,j+1)) - else - CFL_trans = (v(i,J,k) * CS%dt_in_T) * (G%dx_Cv(i,J) * G%IareaT(i,j)) - endif + CFL_Iarea = G%IareaT(i,j) + if (v(i,J,k) < 0.0) & + CFL_Iarea = G%IareaT(i,j+1) + + CFL_trans = abs(v(i,J,k) * CS%dt_in_T) * (G%dx_Cv(i,J) * CFL_Iarea) CFL_lin = abs(v(i,J,k) * CS%dt_in_T) * G%IdyCv(i,J) max_CFL(1) = max(max_CFL(1), CFL_trans) max_CFL(2) = max(max_CFL(2), CFL_lin) diff --git a/src/framework/MOM_array_transform.F90 b/src/framework/MOM_array_transform.F90 new file mode 100644 index 0000000000..09d55ad50b --- /dev/null +++ b/src/framework/MOM_array_transform.F90 @@ -0,0 +1,358 @@ +!> Module for supporting the rotation of a field's index map. +!! The implementation of each angle is described below. +!! +!! +90deg: B(i,j) = A(n-j,i) +!! = transpose, then row reverse +!! 180deg: B(i,j) = A(m-i,n-j) +!! = row reversal + column reversal +!! -90deg: B(i,j) = A(j,m-i) +!! = row reverse, then transpose +!! +!! 90 degree rotations change the shape of the field, and are handled +!! separately from 180 degree rotations. + +module MOM_array_transform + +implicit none + +private +public rotate_array +public rotate_array_pair +public rotate_vector +public allocate_rotated_array + + +!> Rotate the elements of an array to the rotated set of indices. +!! Rotation is applied across the first and second axes of the array. +interface rotate_array + module procedure rotate_array_real_2d + module procedure rotate_array_real_3d + module procedure rotate_array_real_4d + module procedure rotate_array_integer + module procedure rotate_array_logical +end interface rotate_array + + +!> Rotate a pair of arrays which map to a rotated set of indices. +!! Rotation is applied across the first and second axes of the array. +!! This rotation should be applied when one field is mapped onto the other. +!! For example, a tracer indexed along u or v face points will map from one +!! to the other after a quarter turn, and back onto itself after a half turn. +interface rotate_array_pair + module procedure rotate_array_pair_real_2d + module procedure rotate_array_pair_real_3d + module procedure rotate_array_pair_integer +end interface rotate_array_pair + + +!> Rotate an array pair representing the components of a vector. +!! Rotation is applied across the first and second axes of the array. +!! This rotation should be applied when the fields satisfy vector +!! transformation rules. For example, the u and v components of a velocity +!! will map from one to the other for quarter turns, with a sign change in one +!! component. A half turn will map elements onto themselves with sign changes +!! in both components. +interface rotate_vector + module procedure rotate_vector_real_2d + module procedure rotate_vector_real_3d + module procedure rotate_vector_real_4d +end interface rotate_vector + + +!> Allocate an array based on the rotated index map of an unrotated reference +!! array. +interface allocate_rotated_array + module procedure allocate_rotated_array_real_2d + module procedure allocate_rotated_array_real_3d + module procedure allocate_rotated_array_real_4d + module procedure allocate_rotated_array_integer +end interface allocate_rotated_array + +contains + +!> Rotate the elements of a 2d real array along first and second axes. +subroutine rotate_array_real_2d(A_in, turns, A) + real, intent(in) :: A_in(:,:) !< Unrotated array + integer, intent(in) :: turns !< Number of quarter turns + real, intent(out) :: A(:,:) !< Rotated array + + integer :: m, n + + m = size(A_in, 1) + n = size(A_in, 2) + + select case (modulo(turns, 4)) + case(0) + A = A_in + case(1) + A = transpose(A_in) + A = A(n:1:-1, :) + case(2) + A = A_in(m:1:-1, n:1:-1) + case(3) + A = transpose(A_in(m:1:-1, :)) + end select +end subroutine rotate_array_real_2d + + +!> Rotate the elements of a 3d real array along first and second axes. +subroutine rotate_array_real_3d(A_in, turns, A) + real, intent(in) :: A_in(:,:,:) !< Unrotated array + integer, intent(in) :: turns !< Number of quarter turns + real, intent(out) :: A(:,:,:) !< Rotated array + + integer :: k + + do k = lbound(A_in, 3), ubound(A_in, 3) + call rotate_array(A_in(:,:,k), turns, A(:,:,k)) + enddo +end subroutine rotate_array_real_3d + + +!> Rotate the elements of a 4d real array along first and second axes. +subroutine rotate_array_real_4d(A_in, turns, A) + real, intent(in) :: A_in(:,:,:,:) !< Unrotated array + integer, intent(in) :: turns !< Number of quarter turns + real, intent(out) :: A(:,:,:,:) !< Rotated array + + integer :: n + + do n = lbound(A_in, 4), ubound(A_in, 4) + call rotate_array(A_in(:,:,:,n), turns, A(:,:,:,n)) + enddo +end subroutine rotate_array_real_4d + + +!> Rotate the elements of a 2d integer array along first and second axes. +subroutine rotate_array_integer(A_in, turns, A) + integer, intent(in) :: A_in(:,:) !< Unrotated array + integer, intent(in) :: turns !< Number of quarter turns + integer, intent(out) :: A(:,:) !< Rotated array + + integer :: m, n + + m = size(A_in, 1) + n = size(A_in, 2) + + select case (modulo(turns, 4)) + case(0) + A = A_in + case(1) + A = transpose(A_in) + A = A(n:1:-1, :) + case(2) + A = A_in(m:1:-1, n:1:-1) + case(3) + A = transpose(A_in(m:1:-1, :)) + end select +end subroutine rotate_array_integer + + +!> Rotate the elements of a 2d logical array along first and second axes. +subroutine rotate_array_logical(A_in, turns, A) + logical, intent(in) :: A_in(:,:) !< Unrotated array + integer, intent(in) :: turns !< Number of quarter turns + logical, intent(out) :: A(:,:) !< Rotated array + + integer :: m, n + + m = size(A_in, 1) + n = size(A_in, 2) + + select case (modulo(turns, 4)) + case(0) + A = A_in + case(1) + A = transpose(A_in) + A = A(n:1:-1, :) + case(2) + A = A_in(m:1:-1, n:1:-1) + case(3) + A = transpose(A_in(m:1:-1, :)) + end select +end subroutine rotate_array_logical + + +!> Rotate the elements of a 2d real array pair along first and second axes. +subroutine rotate_array_pair_real_2d(A_in, B_in, turns, A, B) + real, intent(in) :: A_in(:,:) !< Unrotated scalar array pair + real, intent(in) :: B_in(:,:) !< Unrotated scalar array pair + integer, intent(in) :: turns !< Number of quarter turns + real, intent(out) :: A(:,:) !< Rotated scalar array pair + real, intent(out) :: B(:,:) !< Rotated scalar array pair + + if (modulo(turns, 2) /= 0) then + call rotate_array(B_in, turns, A) + call rotate_array(A_in, turns, B) + else + call rotate_array(A_in, turns, A) + call rotate_array(B_in, turns, B) + endif +end subroutine rotate_array_pair_real_2d + + +!> Rotate the elements of a 3d real array pair along first and second axes. +subroutine rotate_array_pair_real_3d(A_in, B_in, turns, A, B) + real, intent(in) :: A_in(:,:,:) !< Unrotated scalar array pair + real, intent(in) :: B_in(:,:,:) !< Unrotated scalar array pair + integer, intent(in) :: turns !< Number of quarter turns + real, intent(out) :: A(:,:,:) !< Rotated scalar array pair + real, intent(out) :: B(:,:,:) !< Rotated scalar array pair + + integer :: k + + do k = lbound(A_in, 3), ubound(A_in, 3) + call rotate_array_pair(A_in(:,:,k), B_in(:,:,k), turns, & + A(:,:,k), B(:,:,k)) + enddo +end subroutine rotate_array_pair_real_3d + + +!> Rotate the elements of a 4d real array pair along first and second axes. +subroutine rotate_array_pair_integer(A_in, B_in, turns, A, B) + integer, intent(in) :: A_in(:,:) !< Unrotated scalar array pair + integer, intent(in) :: B_in(:,:) !< Unrotated scalar array pair + integer, intent(in) :: turns !< Number of quarter turns + integer, intent(out) :: A(:,:) !< Rotated scalar array pair + integer, intent(out) :: B(:,:) !< Rotated scalar array pair + + if (modulo(turns, 2) /= 0) then + call rotate_array(B_in, turns, A) + call rotate_array(A_in, turns, B) + else + call rotate_array(A_in, turns, A) + call rotate_array(B_in, turns, B) + endif +end subroutine rotate_array_pair_integer + + +!> Rotate the elements of a 2d real vector along first and second axes. +subroutine rotate_vector_real_2d(A_in, B_in, turns, A, B) + real, intent(in) :: A_in(:,:) !< First component of unrotated vector + real, intent(in) :: B_in(:,:) !< Second component of unrotated vector + integer, intent(in) :: turns !< Number of quarter turns + real, intent(out) :: A(:,:) !< First component of rotated vector + real, intent(out) :: B(:,:) !< Second component of unrotated vector + + call rotate_array_pair(A_in, B_in, turns, A, B) + + if (modulo(turns, 4) == 1 .or. modulo(turns, 4) == 2) & + A = -A + + if (modulo(turns, 4) == 2 .or. modulo(turns, 4) == 3) & + B = -B +end subroutine rotate_vector_real_2d + + +!> Rotate the elements of a 3d real vector along first and second axes. +subroutine rotate_vector_real_3d(A_in, B_in, turns, A, B) + real, intent(in) :: A_in(:,:,:) !< First component of unrotated vector + real, intent(in) :: B_in(:,:,:) !< Second component of unrotated vector + integer, intent(in) :: turns !< Number of quarter turns + real, intent(out) :: A(:,:,:) !< First component of rotated vector + real, intent(out) :: B(:,:,:) !< Second component of unrotated vector + + integer :: k + + do k = lbound(A_in, 3), ubound(A_in, 3) + call rotate_vector(A_in(:,:,k), B_in(:,:,k), turns, A(:,:,k), B(:,:,k)) + enddo +end subroutine rotate_vector_real_3d + + +!> Rotate the elements of a 4d real vector along first and second axes. +subroutine rotate_vector_real_4d(A_in, B_in, turns, A, B) + real, intent(in) :: A_in(:,:,:,:) !< First component of unrotated vector + real, intent(in) :: B_in(:,:,:,:) !< Second component of unrotated vector + integer, intent(in) :: turns !< Number of quarter turns + real, intent(out) :: A(:,:,:,:) !< First component of rotated vector + real, intent(out) :: B(:,:,:,:) !< Second component of unrotated vector + + integer :: n + + do n = lbound(A_in, 4), ubound(A_in, 4) + call rotate_vector(A_in(:,:,:,n), B_in(:,:,:,n), turns, & + A(:,:,:,n), B(:,:,:,n)) + enddo +end subroutine rotate_vector_real_4d + + +!> Allocate a 2d real array on the rotated index map of a reference array. +subroutine allocate_rotated_array_real_2d(A_in, lb, turns, A) + ! NOTE: lb must be declared before A_in + integer, intent(in) :: lb(2) !< Lower index bounds of A_in + real, intent(in) :: A_in(lb(1):, lb(2):) !< Reference array + integer, intent(in) :: turns !< Number of quarter turns + real, allocatable, intent(inout) :: A(:,:) !< Array on rotated index + + integer :: ub(2) + + ub(:) = ubound(A_in) + + if (modulo(turns, 2) /= 0) then + allocate(A(lb(2):ub(2), lb(1):ub(1))) + else + allocate(A(lb(1):ub(1), lb(2):ub(2))) + endif +end subroutine allocate_rotated_array_real_2d + + +!> Allocate a 3d real array on the rotated index map of a reference array. +subroutine allocate_rotated_array_real_3d(A_in, lb, turns, A) + ! NOTE: lb must be declared before A_in + integer, intent(in) :: lb(3) !< Lower index bounds of A_in + real, intent(in) :: A_in(lb(1):, lb(2):, lb(3):) !< Reference array + integer, intent(in) :: turns !< Number of quarter turns + real, allocatable, intent(inout) :: A(:,:,:) !< Array on rotated index + + integer :: ub(3) + + ub(:) = ubound(A_in) + + if (modulo(turns, 2) /= 0) then + allocate(A(lb(2):ub(2), lb(1):ub(1), lb(3):ub(3))) + else + allocate(A(lb(1):ub(1), lb(2):ub(2), lb(3):ub(3))) + endif +end subroutine allocate_rotated_array_real_3d + + +!> Allocate a 4d real array on the rotated index map of a reference array. +subroutine allocate_rotated_array_real_4d(A_in, lb, turns, A) + ! NOTE: lb must be declared before A_in + integer, intent(in) :: lb(4) !< Lower index bounds of A_in + real, intent(in) :: A_in(lb(1):,lb(2):,lb(3):,lb(4):) !< Reference array + integer, intent(in) :: turns !< Number of quarter turns + real, allocatable, intent(inout) :: A(:,:,:,:) !< Array on rotated index + + integer:: ub(4) + + ub(:) = ubound(A_in) + + if (modulo(turns, 2) /= 0) then + allocate(A(lb(2):ub(2), lb(1):ub(1), lb(3):ub(3), lb(4):ub(4))) + else + allocate(A(lb(1):ub(1), lb(2):ub(2), lb(3):ub(3), lb(4):ub(4))) + endif +end subroutine allocate_rotated_array_real_4d + + +!> Allocate a 2d integer array on the rotated index map of a reference array. +subroutine allocate_rotated_array_integer(A_in, lb, turns, A) + integer, intent(in) :: lb(2) !< Lower index bounds of A_in + integer, intent(in) :: A_in(lb(1):,lb(2):) !< Reference array + integer, intent(in) :: turns !< Number of quarter turns + integer, allocatable, intent(inout) :: A(:,:) !< Array on rotated index + + integer :: ub(2) + + ub(:) = ubound(A_in) + + if (modulo(turns, 2) /= 0) then + allocate(A(lb(2):ub(2), lb(1):ub(1))) + else + allocate(A(lb(1):ub(1), lb(2):ub(2))) + endif +end subroutine allocate_rotated_array_integer + +end module MOM_array_transform diff --git a/src/framework/MOM_checksums.F90 b/src/framework/MOM_checksums.F90 index ad269f3530..3cc1f316e2 100644 --- a/src/framework/MOM_checksums.F90 +++ b/src/framework/MOM_checksums.F90 @@ -3,12 +3,13 @@ module MOM_checksums ! This file is part of MOM6. See LICENSE.md for the license. +use MOM_array_transform, only: rotate_array, rotate_array_pair, rotate_vector use MOM_coms, only : PE_here, root_PE, num_PEs, sum_across_PEs use MOM_coms, only : min_across_PEs, max_across_PEs use MOM_coms, only : reproducing_sum use MOM_error_handler, only : MOM_error, FATAL, is_root_pe use MOM_file_parser, only : log_version, param_file_type -use MOM_hor_index, only : hor_index_type +use MOM_hor_index, only : hor_index_type, rotate_hor_index use iso_fortran_env, only: error_unit @@ -191,68 +192,126 @@ subroutine subStats(array, aMean, aMin, aMax) enddo aMean = sum(array(:)) / real(n) end subroutine subStats - end subroutine zchksum !> Checksums on a pair of 2d arrays staggered at tracer points. subroutine chksum_pair_h_2d(mesg, arrayA, arrayB, HI, haloshift, omit_corners, & - scale, logunit) + scale, logunit, scalar_pair) character(len=*), intent(in) :: mesg !< Identifying messages - type(hor_index_type), intent(in) :: HI !< A horizontal index type - real, dimension(HI%isd:,HI%jsd:), intent(in) :: arrayA !< The first array to be checksummed - real, dimension(HI%isd:,HI%jsd:), intent(in) :: arrayB !< The second array to be checksummed + type(hor_index_type), target, intent(in) :: HI !< A horizontal index type + real, dimension(HI%isd:,HI%jsd:), target, intent(in) :: arrayA !< The first array to be checksummed + real, dimension(HI%isd:,HI%jsd:), target, intent(in) :: arrayB !< The second array to be checksummed integer, optional, intent(in) :: haloshift !< The width of halos to check (default 0) logical, optional, intent(in) :: omit_corners !< If true, avoid checking diagonal shifts real, optional, intent(in) :: scale !< A scaling factor for this array. integer, optional, intent(in) :: logunit !< IO unit for checksum logging + logical, optional, intent(in) :: scalar_pair !< If true, then the arrays describe + !! a scalar, rather than vector + logical :: vector_pair + integer :: turns + type(hor_index_type), pointer :: HI_in + real, dimension(:,:), pointer :: arrayA_in, arrayB_in + + vector_pair = .true. + if (present(scalar_pair)) vector_pair = .not. scalar_pair + + turns = HI%turns + if (modulo(turns, 4) /= 0) then + ! Rotate field back to the input grid + allocate(HI_in) + call rotate_hor_index(HI, -turns, HI_in) + allocate(arrayA_in(HI_in%isd:HI_in%ied, HI_in%jsd:HI_in%jed)) + allocate(arrayB_in(HI_in%isd:HI_in%ied, HI_in%jsd:HI_in%jed)) + + if (vector_pair) then + call rotate_vector(arrayA, arrayB, -turns, arrayA_in, arrayB_in) + else + call rotate_array_pair(arrayA, arrayB, -turns, arrayA_in, arrayB_in) + endif + else + HI_in => HI + arrayA_in => arrayA + arrayB_in => arrayB + endif if (present(haloshift)) then - call chksum_h_2d(arrayA, 'x '//mesg, HI, haloshift, omit_corners, & + call chksum_h_2d(arrayA_in, 'x '//mesg, HI_in, haloshift, omit_corners, & scale=scale, logunit=logunit) - call chksum_h_2d(arrayB, 'y '//mesg, HI, haloshift, omit_corners, & + call chksum_h_2d(arrayB_in, 'y '//mesg, HI_in, haloshift, omit_corners, & scale=scale, logunit=logunit) else - call chksum_h_2d(arrayA, 'x '//mesg, HI, scale=scale, logunit=logunit) - call chksum_h_2d(arrayB, 'y '//mesg, HI, scale=scale, logunit=logunit) + call chksum_h_2d(arrayA_in, 'x '//mesg, HI_in, scale=scale, logunit=logunit) + call chksum_h_2d(arrayB_in, 'y '//mesg, HI_in, scale=scale, logunit=logunit) endif - end subroutine chksum_pair_h_2d !> Checksums on a pair of 3d arrays staggered at tracer points. subroutine chksum_pair_h_3d(mesg, arrayA, arrayB, HI, haloshift, omit_corners, & - scale, logunit) + scale, logunit, scalar_pair) character(len=*), intent(in) :: mesg !< Identifying messages - type(hor_index_type), intent(in) :: HI !< A horizontal index type - real, dimension(HI%isd:,HI%jsd:, :), intent(in) :: arrayA !< The first array to be checksummed - real, dimension(HI%isd:,HI%jsd:, :), intent(in) :: arrayB !< The second array to be checksummed + type(hor_index_type), target, intent(in) :: HI !< A horizontal index type + real, dimension(HI%isd:,HI%jsd:, :), target, intent(in) :: arrayA !< The first array to be checksummed + real, dimension(HI%isd:,HI%jsd:, :), target, intent(in) :: arrayB !< The second array to be checksummed integer, optional, intent(in) :: haloshift !< The width of halos to check (default 0) logical, optional, intent(in) :: omit_corners !< If true, avoid checking diagonal shifts real, optional, intent(in) :: scale !< A scaling factor for this array. integer, optional, intent(in) :: logunit !< IO unit for checksum logging + logical, optional, intent(in) :: scalar_pair !< If true, then the arrays describe + !! a scalar, rather than vector + logical :: vector_pair + integer :: turns + type(hor_index_type), pointer :: HI_in + real, dimension(:,:,:), pointer :: arrayA_in, arrayB_in + + vector_pair = .true. + if (present(scalar_pair)) vector_pair = .not. scalar_pair + + turns = HI%turns + if (modulo(turns, 4) /= 0) then + ! Rotate field back to the input grid + allocate(HI_in) + call rotate_hor_index(HI, -turns, HI_in) + allocate(arrayA_in(HI_in%isd:HI_in%ied, HI_in%jsd:HI_in%jed, size(arrayA, 3))) + allocate(arrayB_in(HI_in%isd:HI_in%ied, HI_in%jsd:HI_in%jed, size(arrayB, 3))) + + if (vector_pair) then + call rotate_vector(arrayA, arrayB, -turns, arrayA_in, arrayB_in) + else + call rotate_array_pair(arrayA, arrayB, -turns, arrayA_in, arrayB_in) + endif + else + HI_in => HI + arrayA_in => arrayA + arrayB_in => arrayB + endif + if (present(haloshift)) then - call chksum_h_3d(arrayA, 'x '//mesg, HI, haloshift, omit_corners, & + call chksum_h_3d(arrayA_in, 'x '//mesg, HI_in, haloshift, omit_corners, & scale=scale, logunit=logunit) - call chksum_h_3d(arrayB, 'y '//mesg, HI, haloshift, omit_corners, & + call chksum_h_3d(arrayB_in, 'y '//mesg, HI_in, haloshift, omit_corners, & scale=scale, logunit=logunit) else - call chksum_h_3d(arrayA, 'x '//mesg, HI, scale=scale, logunit=logunit) - call chksum_h_3d(arrayB, 'y '//mesg, HI, scale=scale, logunit=logunit) + call chksum_h_3d(arrayA_in, 'x '//mesg, HI_in, scale=scale, logunit=logunit) + call chksum_h_3d(arrayB_in, 'y '//mesg, HI_in, scale=scale, logunit=logunit) endif + ! NOTE: automatic deallocation of array[AB]_in end subroutine chksum_pair_h_3d !> Checksums a 2d array staggered at tracer points. -subroutine chksum_h_2d(array, mesg, HI, haloshift, omit_corners, scale, logunit) - type(hor_index_type), intent(in) :: HI !< A horizontal index type - real, dimension(HI%isd:,HI%jsd:), intent(in) :: array !< The array to be checksummed +subroutine chksum_h_2d(array_m, mesg, HI_m, haloshift, omit_corners, scale, logunit) + type(hor_index_type), target, intent(in) :: HI_m !< Horizontal index bounds of the model grid + real, dimension(HI_m%isd:,HI_m%jsd:), target, intent(in) :: array_m !< Field array on the model grid character(len=*), intent(in) :: mesg !< An identifying message integer, optional, intent(in) :: haloshift !< The width of halos to check (default 0) logical, optional, intent(in) :: omit_corners !< If true, avoid checking diagonal shifts real, optional, intent(in) :: scale !< A scaling factor for this array. integer, optional, intent(in) :: logunit !< IO unit for checksum logging + real, pointer :: array(:,:) ! Field array on the input grid real, allocatable, dimension(:,:) :: rescaled_array + type(hor_index_type), pointer :: HI ! Horizontal index bounds of the input grid real :: scaling integer :: iounit !< Log IO unit integer :: i, j @@ -260,6 +319,19 @@ subroutine chksum_h_2d(array, mesg, HI, haloshift, omit_corners, scale, logunit) integer :: bc0, bcSW, bcSE, bcNW, bcNE, hshift integer :: bcN, bcS, bcE, bcW logical :: do_corners + integer :: turns ! Quarter turns from input to model grid + + ! Rotate array to the input grid + turns = HI_m%turns + if (modulo(turns, 4) /= 0) then + allocate(HI) + call rotate_hor_index(HI_m, -turns, HI) + allocate(array(HI%isd:HI%ied, HI%jsd:HI%jed)) + call rotate_array(array_m, -turns, array) + else + HI => HI_m + array => array_m + endif if (checkForNaNs) then if (is_NaN(array(HI%isc:HI%iec,HI%jsc:HI%jec))) & @@ -373,31 +445,59 @@ end subroutine chksum_h_2d !> Checksums on a pair of 2d arrays staggered at q-points. subroutine chksum_pair_B_2d(mesg, arrayA, arrayB, HI, haloshift, symmetric, & - omit_corners, scale, logunit) + omit_corners, scale, logunit, scalar_pair) character(len=*), intent(in) :: mesg !< Identifying messages - type(hor_index_type), intent(in) :: HI !< A horizontal index type - real, dimension(HI%isd:,HI%jsd:), intent(in) :: arrayA !< The first array to be checksummed - real, dimension(HI%isd:,HI%jsd:), intent(in) :: arrayB !< The second array to be checksummed + type(hor_index_type), target, intent(in) :: HI !< A horizontal index type + real, dimension(HI%isd:,HI%jsd:), target, intent(in) :: arrayA !< The first array to be checksummed + real, dimension(HI%isd:,HI%jsd:), target, intent(in) :: arrayB !< The second array to be checksummed logical, optional, intent(in) :: symmetric !< If true, do the checksums on the full !! symmetric computational domain. integer, optional, intent(in) :: haloshift !< The width of halos to check (default 0) logical, optional, intent(in) :: omit_corners !< If true, avoid checking diagonal shifts real, optional, intent(in) :: scale !< A scaling factor for this array. integer, optional, intent(in) :: logunit !< IO unit for checksum logging + logical, optional, intent(in) :: scalar_pair !< If true, then the arrays describe + !! a scalar, rather than vector logical :: sym + logical :: vector_pair + integer :: turns + type(hor_index_type), pointer :: HI_in + real, dimension(:,:), pointer :: arrayA_in, arrayB_in + + vector_pair = .true. + if (present(scalar_pair)) vector_pair = .not. scalar_pair + + turns = HI%turns + if (modulo(turns, 4) /= 0) then + ! Rotate field back to the input grid + allocate(HI_in) + call rotate_hor_index(HI, -turns, HI_in) + allocate(arrayA_in(HI_in%IsdB:HI_in%IedB, HI_in%JsdB:HI_in%JedB)) + allocate(arrayB_in(HI_in%IsdB:HI_in%IedB, HI_in%JsdB:HI_in%JedB)) + + if (vector_pair) then + call rotate_vector(arrayA, arrayB, -turns, arrayA_in, arrayB_in) + else + call rotate_array_pair(arrayA, arrayB, -turns, arrayA_in, arrayB_in) + endif + else + HI_in => HI + arrayA_in => arrayA + arrayB_in => arrayB + endif sym = .false. ; if (present(symmetric)) sym = symmetric if (present(haloshift)) then - call chksum_B_2d(arrayA, 'x '//mesg, HI, haloshift, symmetric=sym, & + call chksum_B_2d(arrayA_in, 'x '//mesg, HI_in, haloshift, symmetric=sym, & omit_corners=omit_corners, scale=scale, logunit=logunit) - call chksum_B_2d(arrayB, 'y '//mesg, HI, haloshift, symmetric=sym, & + call chksum_B_2d(arrayB_in, 'y '//mesg, HI_in, haloshift, symmetric=sym, & omit_corners=omit_corners, scale=scale, logunit=logunit) else - call chksum_B_2d(arrayA, 'x '//mesg, HI, symmetric=sym, scale=scale, & + call chksum_B_2d(arrayA_in, 'x '//mesg, HI_in, symmetric=sym, scale=scale, & logunit=logunit) - call chksum_B_2d(arrayB, 'y '//mesg, HI, symmetric=sym, scale=scale, & + call chksum_B_2d(arrayB_in, 'y '//mesg, HI_in, symmetric=sym, scale=scale, & logunit=logunit) endif @@ -405,40 +505,67 @@ end subroutine chksum_pair_B_2d !> Checksums on a pair of 3d arrays staggered at q-points. subroutine chksum_pair_B_3d(mesg, arrayA, arrayB, HI, haloshift, symmetric, & - omit_corners, scale, logunit) + omit_corners, scale, logunit, scalar_pair) character(len=*), intent(in) :: mesg !< Identifying messages - type(hor_index_type), intent(in) :: HI !< A horizontal index type - real, dimension(HI%IsdB:,HI%JsdB:, :), intent(in) :: arrayA !< The first array to be checksummed - real, dimension(HI%IsdB:,HI%JsdB:, :), intent(in) :: arrayB !< The second array to be checksummed + type(hor_index_type), target, intent(in) :: HI !< A horizontal index type + real, dimension(HI%IsdB:,HI%JsdB:, :), target, intent(in) :: arrayA !< The first array to be checksummed + real, dimension(HI%IsdB:,HI%JsdB:, :), target, intent(in) :: arrayB !< The second array to be checksummed integer, optional, intent(in) :: haloshift !< The width of halos to check (default 0) logical, optional, intent(in) :: symmetric !< If true, do the checksums on the full !! symmetric computational domain. logical, optional, intent(in) :: omit_corners !< If true, avoid checking diagonal shifts real, optional, intent(in) :: scale !< A scaling factor for this array. integer, optional, intent(in) :: logunit !< IO unit for checksum logging + logical, optional, intent(in) :: scalar_pair !< If true, then the arrays describe + !! a scalar, rather than vector logical :: sym + logical :: vector_pair + integer :: turns + type(hor_index_type), pointer :: HI_in + real, dimension(:,:,:), pointer :: arrayA_in, arrayB_in + + vector_pair = .true. + if (present(scalar_pair)) vector_pair = .not. scalar_pair + + turns = HI%turns + if (modulo(turns, 4) /= 0) then + ! Rotate field back to the input grid + allocate(HI_in) + call rotate_hor_index(HI, -turns, HI_in) + allocate(arrayA_in(HI_in%IsdB:HI_in%IedB, HI_in%JsdB:HI_in%JedB, size(arrayA, 3))) + allocate(arrayB_in(HI_in%IsdB:HI_in%IedB, HI_in%JsdB:HI_in%JedB, size(arrayB, 3))) + + if (vector_pair) then + call rotate_vector(arrayA, arrayB, -turns, arrayA_in, arrayB_in) + else + call rotate_array_pair(arrayA, arrayB, -turns, arrayA_in, arrayB_in) + endif + else + HI_in => HI + arrayA_in => arrayA + arrayB_in => arrayB + endif if (present(haloshift)) then - call chksum_B_3d(arrayA, 'x '//mesg, HI, haloshift, symmetric, & + call chksum_B_3d(arrayA_in, 'x '//mesg, HI_in, haloshift, symmetric, & omit_corners, scale=scale, logunit=logunit) - call chksum_B_3d(arrayB, 'y '//mesg, HI, haloshift, symmetric, & + call chksum_B_3d(arrayB_in, 'y '//mesg, HI_in, haloshift, symmetric, & omit_corners, scale=scale, logunit=logunit) else - call chksum_B_3d(arrayA, 'x '//mesg, HI, symmetric=symmetric, scale=scale, & + call chksum_B_3d(arrayA_in, 'x '//mesg, HI_in, symmetric=symmetric, scale=scale, & logunit=logunit) - call chksum_B_3d(arrayB, 'y '//mesg, HI, symmetric=symmetric, scale=scale, & + call chksum_B_3d(arrayB_in, 'y '//mesg, HI_in, symmetric=symmetric, scale=scale, & logunit=logunit) endif - end subroutine chksum_pair_B_3d !> Checksums a 2d array staggered at corner points. -subroutine chksum_B_2d(array, mesg, HI, haloshift, symmetric, omit_corners, & +subroutine chksum_B_2d(array_m, mesg, HI_m, haloshift, symmetric, omit_corners, & scale, logunit) - type(hor_index_type), intent(in) :: HI !< A horizontal index type - real, dimension(HI%IsdB:,HI%JsdB:), & - intent(in) :: array !< The array to be checksummed + type(hor_index_type), target, intent(in) :: HI_m !< A horizontal index type + real, dimension(HI_m%IsdB:,HI_m%JsdB:), & + target, intent(in) :: array_m !< The array to be checksummed character(len=*), intent(in) :: mesg !< An identifying message integer, optional, intent(in) :: haloshift !< The width of halos to check (default 0) logical, optional, intent(in) :: symmetric !< If true, do the checksums on the @@ -447,7 +574,9 @@ subroutine chksum_B_2d(array, mesg, HI, haloshift, symmetric, omit_corners, & real, optional, intent(in) :: scale !< A scaling factor for this array. integer, optional, intent(in) :: logunit !< IO unit for checksum logging + real, pointer :: array(:,:) ! Field array on the input grid real, allocatable, dimension(:,:) :: rescaled_array + type(hor_index_type), pointer :: HI ! Horizontal index bounds of the input grid real :: scaling integer :: iounit !< Log IO unit integer :: i, j, Is, Js @@ -455,6 +584,19 @@ subroutine chksum_B_2d(array, mesg, HI, haloshift, symmetric, omit_corners, & integer :: bc0, bcSW, bcSE, bcNW, bcNE, hshift integer :: bcN, bcS, bcE, bcW logical :: do_corners, sym, sym_stats + integer :: turns ! Quarter turns from input to model grid + + ! Rotate array to the input grid + turns = HI_m%turns + if (modulo(turns, 4) /= 0) then + allocate(HI) + call rotate_hor_index(HI_m, -turns, HI) + allocate(array(HI%IsdB:HI%IedB, HI%JsdB:HI%JedB)) + call rotate_array(array_m, -turns, array) + else + HI => HI_m + array => array_m + endif if (checkForNaNs) then if (is_NaN(array(HI%IscB:HI%IecB,HI%JscB:HI%JecB))) & @@ -585,65 +727,119 @@ end subroutine chksum_B_2d !> Checksums a pair of 2d velocity arrays staggered at C-grid locations subroutine chksum_uv_2d(mesg, arrayU, arrayV, HI, haloshift, symmetric, & - omit_corners, scale, logunit) + omit_corners, scale, logunit, scalar_pair) character(len=*), intent(in) :: mesg !< Identifying messages - type(hor_index_type), intent(in) :: HI !< A horizontal index type - real, dimension(HI%IsdB:,HI%jsd:), intent(in) :: arrayU !< The u-component array to be checksummed - real, dimension(HI%isd:,HI%JsdB:), intent(in) :: arrayV !< The v-component array to be checksummed + type(hor_index_type), target, intent(in) :: HI !< A horizontal index type + real, dimension(HI%IsdB:,HI%jsd:), target, intent(in) :: arrayU !< The u-component array to be checksummed + real, dimension(HI%isd:,HI%JsdB:), target, intent(in) :: arrayV !< The v-component array to be checksummed integer, optional, intent(in) :: haloshift !< The width of halos to check (default 0) logical, optional, intent(in) :: symmetric !< If true, do the checksums on the full !! symmetric computational domain. logical, optional, intent(in) :: omit_corners !< If true, avoid checking diagonal shifts real, optional, intent(in) :: scale !< A scaling factor for these arrays. integer, optional, intent(in) :: logunit !< IO unit for checksum logging + logical, optional, intent(in) :: scalar_pair !< If true, then the arrays describe a + !! a scalar, rather than vector + logical :: vector_pair + integer :: turns + type(hor_index_type), pointer :: HI_in + real, dimension(:,:), pointer :: arrayU_in, arrayV_in + + vector_pair = .true. + if (present(scalar_pair)) vector_pair = .not. scalar_pair + + turns = HI%turns + if (modulo(turns, 4) /= 0) then + ! Rotate field back to the input grid + allocate(HI_in) + call rotate_hor_index(HI, -turns, HI_in) + allocate(arrayU_in(HI_in%IsdB:HI_in%IedB, HI_in%jsd:HI_in%jed)) + allocate(arrayV_in(HI_in%isd:HI_in%ied, HI_in%JsdB:HI_in%JedB)) + + if (vector_pair) then + call rotate_vector(arrayU, arrayV, -turns, arrayU_in, arrayV_in) + else + call rotate_array_pair(arrayU, arrayV, -turns, arrayU_in, arrayV_in) + endif + else + HI_in => HI + arrayU_in => arrayU + arrayV_in => arrayV + endif if (present(haloshift)) then - call chksum_u_2d(arrayU, 'u '//mesg, HI, haloshift, symmetric, & - omit_corners, scale, logunit=logunit) - call chksum_v_2d(arrayV, 'v '//mesg, HI, haloshift, symmetric, & - omit_corners, scale, logunit=logunit) + call chksum_u_2d(arrayU_in, 'u '//mesg, HI_in, haloshift, symmetric, & + omit_corners, scale=scale, logunit=logunit) + call chksum_v_2d(arrayV_in, 'v '//mesg, HI_in, haloshift, symmetric, & + omit_corners, scale=scale, logunit=logunit) else - call chksum_u_2d(arrayU, 'u '//mesg, HI, symmetric=symmetric, & - logunit=logunit) - call chksum_v_2d(arrayV, 'v '//mesg, HI, symmetric=symmetric, & - logunit=logunit) + call chksum_u_2d(arrayU_in, 'u '//mesg, HI_in, symmetric=symmetric, & + scale=scale, logunit=logunit) + call chksum_v_2d(arrayV_in, 'v '//mesg, HI_in, symmetric=symmetric, & + scale=scale, logunit=logunit) endif - end subroutine chksum_uv_2d !> Checksums a pair of 3d velocity arrays staggered at C-grid locations subroutine chksum_uv_3d(mesg, arrayU, arrayV, HI, haloshift, symmetric, & - omit_corners, scale, logunit) + omit_corners, scale, logunit, scalar_pair) character(len=*), intent(in) :: mesg !< Identifying messages - type(hor_index_type), intent(in) :: HI !< A horizontal index type - real, dimension(HI%IsdB:,HI%jsd:,:), intent(in) :: arrayU !< The u-component array to be checksummed - real, dimension(HI%isd:,HI%JsdB:,:), intent(in) :: arrayV !< The v-component array to be checksummed + type(hor_index_type), target, intent(in) :: HI !< A horizontal index type + real, dimension(HI%IsdB:,HI%jsd:,:), target, intent(in) :: arrayU !< The u-component array to be checksummed + real, dimension(HI%isd:,HI%JsdB:,:), target, intent(in) :: arrayV !< The v-component array to be checksummed integer, optional, intent(in) :: haloshift !< The width of halos to check (default 0) logical, optional, intent(in) :: symmetric !< If true, do the checksums on the full !! symmetric computational domain. logical, optional, intent(in) :: omit_corners !< If true, avoid checking diagonal shifts real, optional, intent(in) :: scale !< A scaling factor for these arrays. integer, optional, intent(in) :: logunit !< IO unit for checksum logging + logical, optional, intent(in) :: scalar_pair !< If true, then the arrays describe a + !! a scalar, rather than vector + logical :: vector_pair + integer :: turns + type(hor_index_type), pointer :: HI_in + real, dimension(:,:,:), pointer :: arrayU_in, arrayV_in + + vector_pair = .true. + if (present(scalar_pair)) vector_pair = .not. scalar_pair + + turns = HI%turns + if (modulo(turns, 4) /= 0) then + ! Rotate field back to the input grid + allocate(HI_in) + call rotate_hor_index(HI, -turns, HI_in) + allocate(arrayU_in(HI_in%IsdB:HI_in%IedB, HI_in%jsd:HI_in%jed, size(arrayU, 3))) + allocate(arrayV_in(HI_in%isd:HI_in%ied, HI_in%JsdB:HI_in%JedB, size(arrayV, 3))) + + if (vector_pair) then + call rotate_vector(arrayU, arrayV, -turns, arrayU_in, arrayV_in) + else + call rotate_array_pair(arrayU, arrayV, -turns, arrayU_in, arrayV_in) + endif + else + HI_in => HI + arrayU_in => arrayU + arrayV_in => arrayV + endif if (present(haloshift)) then - call chksum_u_3d(arrayU, 'u '//mesg, HI, haloshift, symmetric, & - omit_corners, scale, logunit=logunit) - call chksum_v_3d(arrayV, 'v '//mesg, HI, haloshift, symmetric, & - omit_corners, scale, logunit=logunit) + call chksum_u_3d(arrayU_in, 'u '//mesg, HI_in, haloshift, symmetric, & + omit_corners, scale=scale, logunit=logunit) + call chksum_v_3d(arrayV_in, 'v '//mesg, HI_in, haloshift, symmetric, & + omit_corners, scale=scale, logunit=logunit) else - call chksum_u_3d(arrayU, 'u '//mesg, HI, symmetric=symmetric, & - logunit=logunit) - call chksum_v_3d(arrayV, 'v '//mesg, HI, symmetric=symmetric, & - logunit=logunit) + call chksum_u_3d(arrayU_in, 'u '//mesg, HI_in, symmetric=symmetric, & + scale=scale, logunit=logunit) + call chksum_v_3d(arrayV_in, 'v '//mesg, HI_in, symmetric=symmetric, & + scale=scale, logunit=logunit) endif - end subroutine chksum_uv_3d !> Checksums a 2d array staggered at C-grid u points. -subroutine chksum_u_2d(array, mesg, HI, haloshift, symmetric, omit_corners, & +subroutine chksum_u_2d(array_m, mesg, HI_m, haloshift, symmetric, omit_corners, & scale, logunit) - type(hor_index_type), intent(in) :: HI !< A horizontal index type - real, dimension(HI%IsdB:,HI%jsd:), intent(in) :: array !< The array to be checksummed + type(hor_index_type), target, intent(in) :: HI_m !< A horizontal index type + real, dimension(HI_m%IsdB:,HI_m%jsd:), target, intent(in) :: array_m !< The array to be checksummed character(len=*), intent(in) :: mesg !< An identifying message integer, optional, intent(in) :: haloshift !< The width of halos to check (default 0) logical, optional, intent(in) :: symmetric !< If true, do the checksums on the full @@ -652,7 +848,9 @@ subroutine chksum_u_2d(array, mesg, HI, haloshift, symmetric, omit_corners, & real, optional, intent(in) :: scale !< A scaling factor for this array. integer, optional, intent(in) :: logunit !< IO unit for checksum logging + real, pointer :: array(:,:) ! Field array on the input grid real, allocatable, dimension(:,:) :: rescaled_array + type(hor_index_type), pointer :: HI ! Horizontal index bounds of the input grid real :: scaling integer :: iounit !< Log IO unit integer :: i, j, Is @@ -660,6 +858,27 @@ subroutine chksum_u_2d(array, mesg, HI, haloshift, symmetric, omit_corners, & integer :: bc0, bcSW, bcSE, bcNW, bcNE, hshift integer :: bcN, bcS, bcE, bcW logical :: do_corners, sym, sym_stats + integer :: turns ! Quarter turns from input to model grid + + ! Rotate array to the input grid + turns = HI_m%turns + if (modulo(turns, 4) /= 0) then + allocate(HI) + call rotate_hor_index(HI_m, -turns, HI) + if (modulo(turns, 2) /= 0) then + ! Arrays originating from v-points must be handled by vchksum + allocate(array(HI%isd:HI%ied, HI%JsdB:HI%JedB)) + call rotate_array(array_m, -turns, array) + call vchksum(array, mesg, HI, haloshift, symmetric, omit_corners, scale, logunit) + return + else + allocate(array(HI%IsdB:HI%IedB, HI%jsd:HI%jed)) + call rotate_array(array_m, -turns, array) + endif + else + HI => HI_m + array => array_m + endif if (checkForNaNs) then if (is_NaN(array(HI%IscB:HI%IecB,HI%jsc:HI%jec))) & @@ -794,10 +1013,10 @@ end subroutine subStats end subroutine chksum_u_2d !> Checksums a 2d array staggered at C-grid v points. -subroutine chksum_v_2d(array, mesg, HI, haloshift, symmetric, omit_corners, & +subroutine chksum_v_2d(array_m, mesg, HI_m, haloshift, symmetric, omit_corners, & scale, logunit) - type(hor_index_type), intent(in) :: HI !< A horizontal index type - real, dimension(HI%isd:,HI%JsdB:), intent(in) :: array !< The array to be checksummed + type(hor_index_type), target, intent(in) :: HI_m !< A horizontal index type + real, dimension(HI_m%isd:,HI_m%JsdB:), target, intent(in) :: array_m !< The array to be checksummed character(len=*), intent(in) :: mesg !< An identifying message integer, optional, intent(in) :: haloshift !< The width of halos to check (default 0) logical, optional, intent(in) :: symmetric !< If true, do the checksums on the full @@ -806,7 +1025,9 @@ subroutine chksum_v_2d(array, mesg, HI, haloshift, symmetric, omit_corners, & real, optional, intent(in) :: scale !< A scaling factor for this array. integer, optional, intent(in) :: logunit !< IO unit for checksum logging + real, pointer :: array(:,:) ! Field array on the input grid real, allocatable, dimension(:,:) :: rescaled_array + type(hor_index_type), pointer :: HI ! Horizontal index bounds of the input grid real :: scaling integer :: iounit !< Log IO unit integer :: i, j, Js @@ -814,6 +1035,27 @@ subroutine chksum_v_2d(array, mesg, HI, haloshift, symmetric, omit_corners, & integer :: bc0, bcSW, bcSE, bcNW, bcNE, hshift integer :: bcN, bcS, bcE, bcW logical :: do_corners, sym, sym_stats + integer :: turns ! Quarter turns from input to model grid + + ! Rotate array to the input grid + turns = HI_m%turns + if (modulo(turns, 4) /= 0) then + allocate(HI) + call rotate_hor_index(HI_m, -turns, HI) + if (modulo(turns, 2) /= 0) then + ! Arrays originating from u-points must be handled by uchksum + allocate(array(HI%IsdB:HI%IedB, HI%jsd:HI%jed)) + call rotate_array(array_m, -turns, array) + call uchksum(array, mesg, HI, haloshift, symmetric, omit_corners, scale, logunit) + return + else + allocate(array(HI%isd:HI%ied, HI%JsdB:HI%JedB)) + call rotate_array(array_m, -turns, array) + endif + else + HI => HI_m + array => array_m + endif if (checkForNaNs) then if (is_NaN(array(HI%isc:HI%iec,HI%JscB:HI%JecB))) & @@ -948,16 +1190,18 @@ end subroutine subStats end subroutine chksum_v_2d !> Checksums a 3d array staggered at tracer points. -subroutine chksum_h_3d(array, mesg, HI, haloshift, omit_corners, scale, logunit) - type(hor_index_type), intent(in) :: HI !< A horizontal index type - real, dimension(HI%isd:,HI%jsd:,:), intent(in) :: array !< The array to be checksummed +subroutine chksum_h_3d(array_m, mesg, HI_m, haloshift, omit_corners, scale, logunit) + type(hor_index_type), target, intent(in) :: HI_m !< A horizontal index type + real, dimension(HI_m%isd:,HI_m%jsd:,:), target, intent(in) :: array_m !< The array to be checksummed character(len=*), intent(in) :: mesg !< An identifying message integer, optional, intent(in) :: haloshift !< The width of halos to check (default 0) logical, optional, intent(in) :: omit_corners !< If true, avoid checking diagonal shifts real, optional, intent(in) :: scale !< A scaling factor for this array. integer, optional, intent(in) :: logunit !< IO unit for checksum logging + real, pointer :: array(:,:,:) ! Field array on the input grid real, allocatable, dimension(:,:,:) :: rescaled_array + type(hor_index_type), pointer :: HI ! Horizontal index bounds of the input grid real :: scaling integer :: iounit !< Log IO unit integer :: i, j, k @@ -965,6 +1209,19 @@ subroutine chksum_h_3d(array, mesg, HI, haloshift, omit_corners, scale, logunit) integer :: bc0, bcSW, bcSE, bcNW, bcNE, hshift integer :: bcN, bcS, bcE, bcW logical :: do_corners + integer :: turns ! Quarter turns from input to model grid + + ! Rotate array to the input grid + turns = HI_m%turns + if (modulo(turns, 4) /= 0) then + allocate(HI) + call rotate_hor_index(HI_m, -turns, HI) + allocate(array(HI%isd:HI%ied, HI%jsd:HI%jed, size(array_m, 3))) + call rotate_array(array_m, -turns, array) + else + HI => HI_m + array => array_m + endif if (checkForNaNs) then if (is_NaN(array(HI%isc:HI%iec,HI%jsc:HI%jec,:))) & @@ -1080,10 +1337,10 @@ end subroutine subStats end subroutine chksum_h_3d !> Checksums a 3d array staggered at corner points. -subroutine chksum_B_3d(array, mesg, HI, haloshift, symmetric, omit_corners, & +subroutine chksum_B_3d(array_m, mesg, HI_m, haloshift, symmetric, omit_corners, & scale, logunit) - type(hor_index_type), intent(in) :: HI !< A horizontal index type - real, dimension(HI%IsdB:,HI%JsdB:,:), intent(in) :: array !< The array to be checksummed + type(hor_index_type), target, intent(in) :: HI_m !< A horizontal index type + real, dimension(HI_m%IsdB:,HI_m%JsdB:,:), target, intent(in) :: array_m !< The array to be checksummed character(len=*), intent(in) :: mesg !< An identifying message integer, optional, intent(in) :: haloshift !< The width of halos to check (default 0) logical, optional, intent(in) :: symmetric !< If true, do the checksums on the full @@ -1092,7 +1349,9 @@ subroutine chksum_B_3d(array, mesg, HI, haloshift, symmetric, omit_corners, & real, optional, intent(in) :: scale !< A scaling factor for this array. integer, optional, intent(in) :: logunit !< IO unit for checksum logging + real, pointer :: array(:,:,:) ! Field array on the input grid real, allocatable, dimension(:,:,:) :: rescaled_array + type(hor_index_type), pointer :: HI ! Horizontal index bounds of the input grid real :: scaling integer :: iounit !< Log IO unit integer :: i, j, k, Is, Js @@ -1100,6 +1359,19 @@ subroutine chksum_B_3d(array, mesg, HI, haloshift, symmetric, omit_corners, & integer :: bc0, bcSW, bcSE, bcNW, bcNE, hshift integer :: bcN, bcS, bcE, bcW logical :: do_corners, sym, sym_stats + integer :: turns ! Quarter turns from input to model grid + + ! Rotate array to the input grid + turns = HI_m%turns + if (modulo(turns, 4) /= 0) then + allocate(HI) + call rotate_hor_index(HI_m, -turns, HI) + allocate(array(HI%IsdB:HI%IedB, HI%JsdB:HI%JedB, size(array_m, 3))) + call rotate_array(array_m, -turns, array) + else + HI => HI_m + array => array_m + endif if (checkForNaNs) then if (is_NaN(array(HI%IscB:HI%IecB,HI%JscB:HI%JecB,:))) & @@ -1235,10 +1507,10 @@ end subroutine subStats end subroutine chksum_B_3d !> Checksums a 3d array staggered at C-grid u points. -subroutine chksum_u_3d(array, mesg, HI, haloshift, symmetric, omit_corners, & +subroutine chksum_u_3d(array_m, mesg, HI_m, haloshift, symmetric, omit_corners, & scale, logunit) - type(hor_index_type), intent(in) :: HI !< A horizontal index type - real, dimension(HI%isdB:,HI%Jsd:,:), intent(in) :: array !< The array to be checksummed + type(hor_index_type), target, intent(in) :: HI_m !< A horizontal index type + real, dimension(HI_m%isdB:,HI_m%Jsd:,:), target, intent(in) :: array_m !< The array to be checksummed character(len=*), intent(in) :: mesg !< An identifying message integer, optional, intent(in) :: haloshift !< The width of halos to check (default 0) logical, optional, intent(in) :: symmetric !< If true, do the checksums on the full @@ -1247,7 +1519,9 @@ subroutine chksum_u_3d(array, mesg, HI, haloshift, symmetric, omit_corners, & real, optional, intent(in) :: scale !< A scaling factor for this array. integer, optional, intent(in) :: logunit !< IO unit for checksum logging + real, pointer :: array(:,:,:) ! Field array on the input grid real, allocatable, dimension(:,:,:) :: rescaled_array + type(hor_index_type), pointer :: HI ! Horizontal index bounds of the input grid real :: scaling integer :: iounit !< Log IO unit integer :: i, j, k, Is @@ -1255,6 +1529,27 @@ subroutine chksum_u_3d(array, mesg, HI, haloshift, symmetric, omit_corners, & integer :: bc0, bcSW, bcSE, bcNW, bcNE, hshift integer :: bcN, bcS, bcE, bcW logical :: do_corners, sym, sym_stats + integer :: turns ! Quarter turns from input to model grid + + ! Rotate array to the input grid + turns = HI_m%turns + if (modulo(turns, 4) /= 0) then + allocate(HI) + call rotate_hor_index(HI_m, -turns, HI) + if (modulo(turns, 2) /= 0) then + ! Arrays originating from v-points must be handled by vchksum + allocate(array(HI%isd:HI%ied, HI%JsdB:HI%JedB, size(array_m, 3))) + call rotate_array(array_m, -turns, array) + call vchksum(array, mesg, HI, haloshift, symmetric, omit_corners, scale, logunit) + return + else + allocate(array(HI%IsdB:HI%IedB, HI%jsd:HI%jed, size(array_m, 3))) + call rotate_array(array_m, -turns, array) + endif + else + HI => HI_m + array => array_m + endif if (checkForNaNs) then if (is_NaN(array(HI%IscB:HI%IecB,HI%jsc:HI%jec,:))) & @@ -1389,10 +1684,10 @@ end subroutine subStats end subroutine chksum_u_3d !> Checksums a 3d array staggered at C-grid v points. -subroutine chksum_v_3d(array, mesg, HI, haloshift, symmetric, omit_corners, & +subroutine chksum_v_3d(array_m, mesg, HI_m, haloshift, symmetric, omit_corners, & scale, logunit) - type(hor_index_type), intent(in) :: HI !< A horizontal index type - real, dimension(HI%isd:,HI%JsdB:,:), intent(in) :: array !< The array to be checksummed + type(hor_index_type), target, intent(in) :: HI_m !< A horizontal index type + real, dimension(HI_m%isd:,HI_m%JsdB:,:), target, intent(in) :: array_m !< The array to be checksummed character(len=*), intent(in) :: mesg !< An identifying message integer, optional, intent(in) :: haloshift !< The width of halos to check (default 0) logical, optional, intent(in) :: symmetric !< If true, do the checksums on the full @@ -1401,7 +1696,9 @@ subroutine chksum_v_3d(array, mesg, HI, haloshift, symmetric, omit_corners, & real, optional, intent(in) :: scale !< A scaling factor for this array. integer, optional, intent(in) :: logunit !< IO unit for checksum logging + real, pointer :: array(:,:,:) ! Field array on the input grid real, allocatable, dimension(:,:,:) :: rescaled_array + type(hor_index_type), pointer :: HI ! Horizontal index bounds of the input grid real :: scaling integer :: iounit !< Log IO unit integer :: i, j, k, Js @@ -1409,6 +1706,27 @@ subroutine chksum_v_3d(array, mesg, HI, haloshift, symmetric, omit_corners, & integer :: bcN, bcS, bcE, bcW real :: aMean, aMin, aMax logical :: do_corners, sym, sym_stats + integer :: turns ! Quarter turns from input to model grid + + ! Rotate array to the input grid + turns = HI_m%turns + if (modulo(turns, 4) /= 0) then + allocate(HI) + call rotate_hor_index(HI_m, -turns, HI) + if (modulo(turns, 2) /= 0) then + ! Arrays originating from u-points must be handled by uchksum + allocate(array(HI%IsdB:HI%IedB, HI%jsd:HI%jed, size(array_m, 3))) + call rotate_array(array_m, -turns, array) + call uchksum(array, mesg, HI, haloshift, symmetric, omit_corners, scale, logunit) + return + else + allocate(array(HI%isd:HI%ied, HI%JsdB:HI%JedB, size(array_m, 3))) + call rotate_array(array_m, -turns, array) + endif + else + HI => HI_m + array => array_m + endif if (checkForNaNs) then if (is_NaN(array(HI%isc:HI%iec,HI%JscB:HI%JecB,:))) & diff --git a/src/framework/MOM_domains.F90 b/src/framework/MOM_domains.F90 index 64fddfe7fc..477ebd70df 100644 --- a/src/framework/MOM_domains.F90 +++ b/src/framework/MOM_domains.F90 @@ -3,6 +3,7 @@ module MOM_domains ! This file is part of MOM6. See LICENSE.md for the license. +use MOM_array_transform, only : rotate_array use MOM_coms, only : PE_here, root_PE, num_PEs, MOM_infra_init, MOM_infra_end use MOM_coms, only : broadcast, sum_across_PEs, min_across_PEs, max_across_PEs use MOM_cpu_clock, only : cpu_clock_begin, cpu_clock_end @@ -1599,7 +1600,7 @@ end subroutine MOM_domains_init !> clone_MD_to_MD copies one MOM_domain_type into another, while allowing !! some properties of the new type to differ from the original one. subroutine clone_MD_to_MD(MD_in, MOM_dom, min_halo, halo_size, symmetric, & - domain_name) + domain_name, turns) type(MOM_domain_type), intent(in) :: MD_in !< An existing MOM_domain type(MOM_domain_type), pointer :: MOM_dom !< A pointer to a MOM_domain that will be !! allocated if it is unassociated, and will have data @@ -1617,10 +1618,15 @@ subroutine clone_MD_to_MD(MD_in, MOM_dom, min_halo, halo_size, symmetric, & character(len=*), & optional, intent(in) :: domain_name !< A name for the new domain, "MOM" !! if missing. + integer, optional, intent(in) :: turns !< Number of quarter turns integer :: global_indices(4) logical :: mask_table_exists character(len=64) :: dom_name + integer :: qturns + + qturns = 0 + if (present(turns)) qturns = turns if (.not.associated(MOM_dom)) then allocate(MOM_dom) @@ -1629,19 +1635,37 @@ subroutine clone_MD_to_MD(MD_in, MOM_dom, min_halo, halo_size, symmetric, & endif ! Save the extra data for creating other domains of different resolution that overlay this domain - MOM_dom%niglobal = MD_in%niglobal ; MOM_dom%njglobal = MD_in%njglobal - MOM_dom%nihalo = MD_in%nihalo ; MOM_dom%njhalo = MD_in%njhalo - MOM_dom%symmetric = MD_in%symmetric MOM_dom%nonblocking_updates = MD_in%nonblocking_updates + MOM_dom%thin_halo_updates = MD_in%thin_halo_updates + + if (modulo(qturns, 2) /= 0) then + MOM_dom%niglobal = MD_in%njglobal ; MOM_dom%njglobal = MD_in%niglobal + MOM_dom%nihalo = MD_in%njhalo ; MOM_dom%njhalo = MD_in%nihalo - MOM_dom%X_FLAGS = MD_in%X_FLAGS ; MOM_dom%Y_FLAGS = MD_in%Y_FLAGS - MOM_dom%layout(:) = MD_in%layout(:) ; MOM_dom%io_layout(:) = MD_in%io_layout(:) + MOM_dom%X_FLAGS = MD_in%Y_FLAGS ; MOM_dom%Y_FLAGS = MD_in%X_FLAGS + MOM_dom%layout(:) = MD_in%layout(2:1:-1) + MOM_dom%io_layout(:) = MD_in%io_layout(2:1:-1) + else + MOM_dom%niglobal = MD_in%niglobal ; MOM_dom%njglobal = MD_in%njglobal + MOM_dom%nihalo = MD_in%nihalo ; MOM_dom%njhalo = MD_in%njhalo + + MOM_dom%X_FLAGS = MD_in%X_FLAGS ; MOM_dom%Y_FLAGS = MD_in%Y_FLAGS + MOM_dom%layout(:) = MD_in%layout(:) + MOM_dom%io_layout(:) = MD_in%io_layout(:) + endif + + global_indices(1) = 1 ; global_indices(2) = MOM_dom%niglobal + global_indices(3) = 1 ; global_indices(4) = MOM_dom%njglobal if (associated(MD_in%maskmap)) then mask_table_exists = .true. allocate(MOM_dom%maskmap(MOM_dom%layout(1), MOM_dom%layout(2))) - MOM_dom%maskmap(:,:) = MD_in%maskmap(:,:) + if (qturns /= 0) then + call rotate_array(MD_in%maskmap(:,:), qturns, MOM_dom%maskmap(:,:)) + else + MOM_dom%maskmap(:,:) = MD_in%maskmap(:,:) + endif else mask_table_exists = .false. endif @@ -1665,19 +1689,34 @@ subroutine clone_MD_to_MD(MD_in, MOM_dom, min_halo, halo_size, symmetric, & dom_name = "MOM" if (present(domain_name)) dom_name = trim(domain_name) - global_indices(1) = 1 ; global_indices(2) = MOM_dom%niglobal - global_indices(3) = 1 ; global_indices(4) = MOM_dom%njglobal if (mask_table_exists) then - call MOM_define_domain( global_indices, MOM_dom%layout, MOM_dom%mpp_domain, & + call MOM_define_domain(global_indices, MOM_dom%layout, MOM_dom%mpp_domain, & xflags=MOM_dom%X_FLAGS, yflags=MOM_dom%Y_FLAGS, & xhalo=MOM_dom%nihalo, yhalo=MOM_dom%njhalo, & - symmetry = MOM_dom%symmetric, name=dom_name, & - maskmap=MOM_dom%maskmap ) + symmetry=MOM_dom%symmetric, name=dom_name, & + maskmap=MOM_dom%maskmap) + + global_indices(2) = global_indices(2) / 2 + global_indices(4) = global_indices(4) / 2 + call MOM_define_domain(global_indices, MOM_dom%layout, & + MOM_dom%mpp_domain_d2, & + xflags=MOM_dom%X_FLAGS, yflags=MOM_dom%Y_FLAGS, & + xhalo=(MOM_dom%nihalo/2), yhalo=(MOM_dom%njhalo/2), & + symmetry=MOM_dom%symmetric, name=dom_name, & + maskmap=MOM_dom%maskmap) else - call MOM_define_domain( global_indices, MOM_dom%layout, MOM_dom%mpp_domain, & + call MOM_define_domain(global_indices, MOM_dom%layout, MOM_dom%mpp_domain, & xflags=MOM_dom%X_FLAGS, yflags=MOM_dom%Y_FLAGS, & xhalo=MOM_dom%nihalo, yhalo=MOM_dom%njhalo, & - symmetry = MOM_dom%symmetric, name=dom_name) + symmetry=MOM_dom%symmetric, name=dom_name) + + global_indices(2) = global_indices(2) / 2 + global_indices(4) = global_indices(4) / 2 + call MOM_define_domain(global_indices, MOM_dom%layout, & + MOM_dom%mpp_domain_d2, & + xflags=MOM_dom%X_FLAGS, yflags=MOM_dom%Y_FLAGS, & + xhalo=(MOM_dom%nihalo/2), yhalo=(MOM_dom%njhalo/2), & + symmetry=MOM_dom%symmetric, name=dom_name) endif if ((MOM_dom%io_layout(1) + MOM_dom%io_layout(2) > 0) .and. & @@ -1691,7 +1730,7 @@ end subroutine clone_MD_to_MD !! domain2d type, while allowing some properties of the new type to differ from !! the original one. subroutine clone_MD_to_d2D(MD_in, mpp_domain, min_halo, halo_size, symmetric, & - domain_name) + domain_name, turns) type(MOM_domain_type), intent(in) :: MD_in !< An existing MOM_domain to be cloned type(domain2d), intent(inout) :: mpp_domain !< The new mpp_domain to be set up integer, dimension(2), & @@ -1707,12 +1746,16 @@ subroutine clone_MD_to_d2D(MD_in, mpp_domain, min_halo, halo_size, symmetric, & character(len=*), & optional, intent(in) :: domain_name !< A name for the new domain, "MOM" !! if missing. + integer, optional, intent(in) :: turns !< If true, swap X and Y axes integer :: global_indices(4), layout(2), io_layout(2) integer :: X_FLAGS, Y_FLAGS, niglobal, njglobal, nihalo, njhalo logical :: symmetric_dom character(len=64) :: dom_name + if (present(turns)) & + call MOM_error(FATAL, "Rotation not supported for MOM_domain to domain2d") + ! Save the extra data for creating other domains of different resolution that overlay this domain niglobal = MD_in%niglobal ; njglobal = MD_in%njglobal nihalo = MD_in%nihalo ; njhalo = MD_in%njhalo diff --git a/src/framework/MOM_dyn_horgrid.F90 b/src/framework/MOM_dyn_horgrid.F90 index ef74a12c9d..141340047d 100644 --- a/src/framework/MOM_dyn_horgrid.F90 +++ b/src/framework/MOM_dyn_horgrid.F90 @@ -173,7 +173,7 @@ module MOM_dyn_horgrid !--------------------------------------------------------------------- !> Allocate memory used by the dyn_horgrid_type and related structures. subroutine create_dyn_horgrid(G, HI, bathymetry_at_vel) - type(dyn_horgrid_type), pointer :: G !< A pointer to the dynamic horizontal grid type + type(dyn_horgrid_type), pointer, intent(inout) :: G !< A pointer to the dynamic horizontal grid type type(hor_index_type), intent(in) :: HI !< A hor_index_type for array extents logical, optional, intent(in) :: bathymetry_at_vel !< If true, there are !! separate values for the basin depths at velocity diff --git a/src/framework/MOM_hor_index.F90 b/src/framework/MOM_hor_index.F90 index db52afcdd8..fc833eeea9 100644 --- a/src/framework/MOM_hor_index.F90 +++ b/src/framework/MOM_hor_index.F90 @@ -10,6 +10,7 @@ module MOM_hor_index implicit none ; private public :: hor_index_init, assignment(=) +public :: rotate_hor_index !> Container for horizontal index ranges for data, computational and global domains type, public :: hor_index_type @@ -49,6 +50,8 @@ module MOM_hor_index integer :: niglobal !< The global number of h-cells in the i-direction integer :: njglobal !< The global number of h-cells in the j-direction + + integer :: turns !< Number of quarter-turn rotations from input to model end type hor_index_type !> Copy the contents of one horizontal index type into another @@ -92,6 +95,7 @@ subroutine hor_index_init(Domain, HI, param_file, local_indexing, index_offset) HI%IedB = HI%ied ; HI%JedB = HI%jed HI%IegB = HI%ieg ; HI%JegB = HI%jeg + HI%turns = 0 end subroutine hor_index_init !> HIT_assign copies one hor_index_type into another. It is accessed via an @@ -110,12 +114,57 @@ subroutine HIT_assign(HI1, HI2) HI1%IsdB = HI2%IsdB ; HI1%IedB = HI2%IedB ; HI1%JsdB = HI2%JsdB ; HI1%JedB = HI2%JedB HI1%IsgB = HI2%IsgB ; HI1%IegB = HI2%IegB ; HI1%JsgB = HI2%JsgB ; HI1%JegB = HI2%JegB + HI1%niglobal = HI2%niglobal ; HI1%njglobal = HI2%njglobal HI1%idg_offset = HI2%idg_offset ; HI1%jdg_offset = HI2%jdg_offset HI1%symmetric = HI2%symmetric - HI1%niglobal = HI2%niglobal ; HI1%njglobal = HI2%njglobal - + HI1%turns = HI2%turns end subroutine HIT_assign +!> Rotate the horizontal index ranges from the input to the output map. +subroutine rotate_hor_index(HI_in, turns, HI) + type(hor_index_type), intent(in) :: HI_in !< Unrotated horizontal indices + integer, intent(in) :: turns !< Number of quarter turns + type(hor_index_type), intent(inout) :: HI !< Rotated horizontal indices + + if (modulo(turns, 2) /= 0) then + HI%isc = HI_in%jsc + HI%iec = HI_in%jec + HI%jsc = HI_in%isc + HI%jec = HI_in%iec + HI%isd = HI_in%jsd + HI%ied = HI_in%jed + HI%jsd = HI_in%isd + HI%jed = HI_in%ied + HI%isg = HI_in%jsg + HI%ieg = HI_in%jeg + HI%jsg = HI_in%isg + HI%jeg = HI_in%ieg + + HI%IscB = HI_in%JscB + HI%IecB = HI_in%JecB + HI%JscB = HI_in%IscB + HI%JecB = HI_in%IecB + HI%IsdB = HI_in%JsdB + HI%IedB = HI_in%JedB + HI%JsdB = HI_in%IsdB + HI%JedB = HI_in%IedB + HI%IsgB = HI_in%JsgB + HI%IegB = HI_in%JegB + HI%JsgB = HI_in%IsgB + HI%JegB = HI_in%IegB + + HI%niglobal = HI_in%njglobal + HI%njglobal = HI_in%niglobal + HI%idg_offset = HI_in%jdg_offset + HI%jdg_offset = HI_in%idg_offset + + HI%symmetric = HI_in%symmetric + else + HI = HI_in + endif + HI%turns = HI_in%turns + turns +end subroutine rotate_hor_index + !> \namespace mom_hor_index !! !! The hor_index_type provides the declarations and loop ranges for almost all data with horizontal extent. diff --git a/src/framework/MOM_horizontal_regridding.F90 b/src/framework/MOM_horizontal_regridding.F90 index 7c19d715db..cd8a04f2fb 100644 --- a/src/framework/MOM_horizontal_regridding.F90 +++ b/src/framework/MOM_horizontal_regridding.F90 @@ -20,8 +20,9 @@ module MOM_horizontal_regridding use MOM_io, only : slasher, vardesc, write_field use MOM_string_functions, only : uppercase use MOM_time_manager, only : time_type, get_external_field_size -use MOM_time_manager, only : init_external_field, time_interp_external +use MOM_time_manager, only : init_external_field use MOM_time_manager, only : get_external_field_axes, get_external_field_missing +use MOM_transform_FMS, only : time_interp_external => rotated_time_interp_external use MOM_variables, only : thermo_var_ptrs use mpp_io_mod, only : axistype use mpp_domains_mod, only : mpp_global_field, mpp_get_compute_domain @@ -658,6 +659,9 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, conversion, G, t real, dimension(SZI_(G),SZJ_(G)) :: tr_outf,tr_prev real, dimension(SZI_(G),SZJ_(G)) :: good2,fill2 real, dimension(SZI_(G),SZJ_(G)) :: nlevs + integer :: turns + + turns = G%HI%turns is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed @@ -753,7 +757,7 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, conversion, G, t if (.not.spongeDataOngrid) then if (is_root_pe()) & - call time_interp_external(fms_id, Time, data_in, verbose=.true.) + call time_interp_external(fms_id, Time, data_in, verbose=.true., turns=turns) ! loop through each data level and interpolate to model grid. ! after interpolating, fill in points which will be needed ! to define the layers @@ -873,7 +877,7 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, conversion, G, t enddo ! kd else - call time_interp_external(fms_id, Time, data_in, verbose=.true.) + call time_interp_external(fms_id, Time, data_in, verbose=.true., turns=turns) do k=1,kd do j=js,je do i=is,ie diff --git a/src/framework/MOM_restart.F90 b/src/framework/MOM_restart.F90 index ec9789c20b..c918f3a9ee 100644 --- a/src/framework/MOM_restart.F90 +++ b/src/framework/MOM_restart.F90 @@ -9,16 +9,18 @@ module MOM_restart use MOM_string_functions, only : lowercase use MOM_grid, only : ocean_grid_type use MOM_io, only : create_file, fieldtype, file_exists, open_file, close_file -use MOM_io, only : write_field, MOM_read_data, read_data, get_filename_appendix +use MOM_io, only : MOM_read_data, read_data, get_filename_appendix use MOM_io, only : get_file_info, get_file_atts, get_file_fields, get_file_times use MOM_io, only : vardesc, var_desc, query_vardesc, modify_vardesc use MOM_io, only : MULTIPLE, NETCDF_FILE, READONLY_FILE, SINGLE_FILE use MOM_io, only : CENTER, CORNER, NORTH_FACE, EAST_FACE use MOM_time_manager, only : time_type, time_type_to_real, real_to_time use MOM_time_manager, only : days_in_month, get_date, set_date +use MOM_transform_FMS, only : mpp_chksum => rotated_mpp_chksum +use MOM_transform_FMS, only : write_field => rotated_write_field use MOM_verticalGrid, only : verticalGrid_type -use mpp_mod, only: mpp_chksum,mpp_pe -use mpp_io_mod, only: mpp_attribute_exist, mpp_get_atts +use mpp_io_mod, only : mpp_attribute_exist, mpp_get_atts +use mpp_mod, only : mpp_pe implicit none ; private @@ -26,6 +28,7 @@ module MOM_restart public save_restart, query_initialized, restart_init_end, vardesc public restart_files_exist, determine_is_new_run, is_new_run public register_restart_field_as_obsolete +public register_restart_pair !> A type for making arrays of pointers to 4-d arrays type p4d @@ -86,6 +89,7 @@ module MOM_restart !! made from a run with a different mask_table than the current run, !! in which case the checksums will not match and cause crash. character(len=240) :: restartfile !< The name or name root for MOM restart files. + integer :: turns !< Number of quarter turns from input to model domain !> An array of descriptions of the registered fields type(field_restart), pointer :: restart_field(:) => NULL() @@ -112,6 +116,13 @@ module MOM_restart module procedure register_restart_field_ptr0d, register_restart_field_0d end interface +!> Register a pair of restart fieilds whose rotations map onto each other +interface register_restart_pair + module procedure register_restart_pair_ptr2d + module procedure register_restart_pair_ptr3d + module procedure register_restart_pair_ptr4d +end interface register_restart_pair + !> Indicate whether a field has been read from a restart file interface query_initialized module procedure query_initialized_name @@ -287,6 +298,67 @@ subroutine register_restart_field_ptr0d(f_ptr, var_desc, mandatory, CS) end subroutine register_restart_field_ptr0d + +!> Register a pair of rotationally equivalent 2d restart fields +subroutine register_restart_pair_ptr2d(a_ptr, b_ptr, a_desc, b_desc, & + mandatory, CS) + real, dimension(:,:), target, intent(in) :: a_ptr !< First field pointer + real, dimension(:,:), target, intent(in) :: b_ptr !< Second field pointer + type(vardesc), intent(in) :: a_desc !< First field descriptor + type(vardesc), intent(in) :: b_desc !< Second field descriptor + logical, intent(in) :: mandatory !< If true, abort if field is missing + type(MOM_restart_CS), pointer :: CS !< MOM restart control structure + + if (modulo(CS%turns, 2) /= 0) then + call register_restart_field(b_ptr, a_desc, mandatory, CS) + call register_restart_field(a_ptr, b_desc, mandatory, CS) + else + call register_restart_field(a_ptr, a_desc, mandatory, CS) + call register_restart_field(b_ptr, b_desc, mandatory, CS) + endif +end subroutine register_restart_pair_ptr2d + + +!> Register a pair of rotationally equivalent 3d restart fields +subroutine register_restart_pair_ptr3d(a_ptr, b_ptr, a_desc, b_desc, & + mandatory, CS) + real, dimension(:,:,:), target, intent(in) :: a_ptr !< First field pointer + real, dimension(:,:,:), target, intent(in) :: b_ptr !< Second field pointer + type(vardesc), intent(in) :: a_desc !< First field descriptor + type(vardesc), intent(in) :: b_desc !< Second field descriptor + logical, intent(in) :: mandatory !< If true, abort if field is missing + type(MOM_restart_CS), pointer :: CS !< MOM restart control structure + + if (modulo(CS%turns, 2) /= 0) then + call register_restart_field(b_ptr, a_desc, mandatory, CS) + call register_restart_field(a_ptr, b_desc, mandatory, CS) + else + call register_restart_field(a_ptr, a_desc, mandatory, CS) + call register_restart_field(b_ptr, b_desc, mandatory, CS) + endif +end subroutine register_restart_pair_ptr3d + + +!> Register a pair of rotationally equivalent 2d restart fields +subroutine register_restart_pair_ptr4d(a_ptr, b_ptr, a_desc, b_desc, & + mandatory, CS) + real, dimension(:,:,:,:), target, intent(in) :: a_ptr !< First field pointer + real, dimension(:,:,:,:), target, intent(in) :: b_ptr !< Second field pointer + type(vardesc), intent(in) :: a_desc !< First field descriptor + type(vardesc), intent(in) :: b_desc !< Second field descriptor + logical, intent(in) :: mandatory !< If true, abort if field is missing + type(MOM_restart_CS), pointer :: CS !< MOM restart control structure + + if (modulo(CS%turns, 2) /= 0) then + call register_restart_field(b_ptr, a_desc, mandatory, CS) + call register_restart_field(a_ptr, b_desc, mandatory, CS) + else + call register_restart_field(a_ptr, a_desc, mandatory, CS) + call register_restart_field(b_ptr, b_desc, mandatory, CS) + endif +end subroutine register_restart_pair_ptr4d + + ! The following provide alternate interfaces to register restarts. !> Register a 4-d field for restarts, providing the metadata as individual arguments @@ -815,6 +887,9 @@ subroutine save_restart(directory, time, G, CS, time_stamped, filename, GV) integer :: length integer(kind=8) :: check_val(CS%max_fields,1) integer :: isL, ieL, jsL, jeL, pos + integer :: turns + + turns = CS%turns if (.not.associated(CS)) call MOM_error(FATAL, "MOM_restart " // & "save_restart: Module must be initialized before it is used.") @@ -927,14 +1002,21 @@ subroutine save_restart(directory, time, G, CS, time_stamped, filename, GV) end select !Prepare the checksum of the restart fields to be written to restart files - call get_checksum_loop_ranges(G, pos, isL, ieL, jsL, jeL) + if (modulo(turns, 2) /= 0) then + call get_checksum_loop_ranges(G, pos, jsL, jeL, isL, ieL) + else + call get_checksum_loop_ranges(G, pos, isL, ieL, jsL, jeL) + endif do m=start_var,next_var-1 if (associated(CS%var_ptr3d(m)%p)) then - check_val(m-start_var+1,1) = mpp_chksum(CS%var_ptr3d(m)%p(isL:ieL,jsL:jeL,:)) + check_val(m-start_var+1,1) = & + mpp_chksum(CS%var_ptr3d(m)%p(isL:ieL,jsL:jeL,:), turns=-turns) elseif (associated(CS%var_ptr2d(m)%p)) then - check_val(m-start_var+1,1) = mpp_chksum(CS%var_ptr2d(m)%p(isL:ieL,jsL:jeL)) + check_val(m-start_var+1,1) = & + mpp_chksum(CS%var_ptr2d(m)%p(isL:ieL,jsL:jeL), turns=-turns) elseif (associated(CS%var_ptr4d(m)%p)) then - check_val(m-start_var+1,1) = mpp_chksum(CS%var_ptr4d(m)%p(isL:ieL,jsL:jeL,:,:)) + check_val(m-start_var+1,1) = & + mpp_chksum(CS%var_ptr4d(m)%p(isL:ieL,jsL:jeL,:,:), turns=-turns) elseif (associated(CS%var_ptr1d(m)%p)) then check_val(m-start_var+1,1) = mpp_chksum(CS%var_ptr1d(m)%p) elseif (associated(CS%var_ptr0d(m)%p)) then @@ -951,16 +1033,15 @@ subroutine save_restart(directory, time, G, CS, time_stamped, filename, GV) endif do m=start_var,next_var-1 - if (associated(CS%var_ptr3d(m)%p)) then call write_field(unit,fields(m-start_var+1), G%Domain%mpp_domain, & - CS%var_ptr3d(m)%p, restart_time) + CS%var_ptr3d(m)%p, restart_time, turns=-turns) elseif (associated(CS%var_ptr2d(m)%p)) then call write_field(unit,fields(m-start_var+1), G%Domain%mpp_domain, & - CS%var_ptr2d(m)%p, restart_time) + CS%var_ptr2d(m)%p, restart_time, turns=-turns) elseif (associated(CS%var_ptr4d(m)%p)) then call write_field(unit,fields(m-start_var+1), G%Domain%mpp_domain, & - CS%var_ptr4d(m)%p, restart_time) + CS%var_ptr4d(m)%p, restart_time, turns=-turns) elseif (associated(CS%var_ptr1d(m)%p)) then call write_field(unit, fields(m-start_var+1), CS%var_ptr1d(m)%p, & restart_time) @@ -1425,6 +1506,8 @@ subroutine restart_init(param_file, CS, restart_root) !! set by RESTARTFILE to enable the use of this module by !! other components than MOM. + logical :: rotate_index + ! This include declares and sets the variable "version". #include "version_variable.h" character(len=40) :: mdl = "MOM_restart" ! This module's name. @@ -1464,6 +1547,16 @@ subroutine restart_init(param_file, CS, restart_root) "in which case the checksums will not match and cause crash.",& default=.true.) + ! Maybe not the best place to do this? + call get_param(param_file, mdl, "ROTATE_INDEX", rotate_index, & + default=.false., do_not_log=.true.) + + CS%turns = 0 + if (rotate_index) then + call get_param(param_file, mdl, "INDEX_TURNS", CS%turns, & + default=1, do_not_log=.true.) + endif + allocate(CS%restart_field(CS%max_fields)) allocate(CS%restart_obsolete(CS%max_fields)) allocate(CS%var_ptr0d(CS%max_fields)) diff --git a/src/framework/MOM_transform_FMS.F90 b/src/framework/MOM_transform_FMS.F90 new file mode 100644 index 0000000000..2af6088c90 --- /dev/null +++ b/src/framework/MOM_transform_FMS.F90 @@ -0,0 +1,399 @@ +!> Support functions and interfaces to permit transformed model domains to +!! interact with FMS operations registered on the non-transformed domains. + +module MOM_transform_FMS + +use horiz_interp_mod, only : horiz_interp_type +use MOM_error_handler, only : MOM_error, FATAL +use MOM_io, only : fieldtype, write_field +use mpp_domains_mod, only : domain2D +use fms_mod, only : mpp_chksum +use time_manager_mod, only : time_type +use time_interp_external_mod, only : time_interp_external + +use MOM_array_transform, only : allocate_rotated_array, rotate_array + +implicit none + +private +public rotated_mpp_chksum +public rotated_write_field +public rotated_time_interp_external + +!> Rotate and compute the FMS (mpp) checksum of a field +interface rotated_mpp_chksum + module procedure rotated_mpp_chksum_real_0d + module procedure rotated_mpp_chksum_real_1d + module procedure rotated_mpp_chksum_real_2d + module procedure rotated_mpp_chksum_real_3d + module procedure rotated_mpp_chksum_real_4d +end interface rotated_mpp_chksum + +!> Rotate and write a registered field to an FMS output file +interface rotated_write_field + module procedure rotated_write_field_real_0d + module procedure rotated_write_field_real_1d + module procedure rotated_write_field_real_2d + module procedure rotated_write_field_real_3d + module procedure rotated_write_field_real_4d +end interface rotated_write_field + +!> Read a field based on model time, and rotate to the model domain +interface rotated_time_interp_external + module procedure rotated_time_interp_external_0d + module procedure rotated_time_interp_external_2d + module procedure rotated_time_interp_external_3d +end interface rotated_time_interp_external + +contains + +! NOTE: No transformations are applied to the 0d and 1d field implementations, +! but are provided to maintain compatibility with the FMS interfaces. + + +!> Compute the FMS (mpp) checksum of a scalar. +!! This function is provided to support the full FMS mpp_chksum interface. +function rotated_mpp_chksum_real_0d(field, pelist, mask_val, turns) & + result(chksum) + real, intent(in) :: field !> Input scalar + integer, optional, intent(in) :: pelist(:) !> PE list of ranks to checksum + real, optional, intent(in) :: mask_val !> FMS mask value + integer, optional, intent(in) :: turns !> Number of quarter turns + integer :: chksum !> FMS checksum of scalar + + if (present(turns)) & + call MOM_error(FATAL, "Rotation not supported for 0d fields.") + + chksum = mpp_chksum(field, pelist=pelist, mask_val=mask_val) +end function rotated_mpp_chksum_real_0d + + +!> Compute the FMS (mpp) checksum of a 1d field. +!! This function is provided to support the full FMS mpp_chksum interface. +function rotated_mpp_chksum_real_1d(field, pelist, mask_val, turns) & + result(chksum) + real, intent(in) :: field(:) !> Input field + integer, optional, intent(in) :: pelist(:) !> PE list of ranks to checksum + real, optional, intent(in) :: mask_val !> FMS mask value + integer, optional, intent(in) :: turns !> Number of quarter-turns + integer :: chksum !> FMS checksum of field + + if (present(turns)) & + call MOM_error(FATAL, "Rotation not supported for 1d fields.") + + chksum = mpp_chksum(field, pelist=pelist, mask_val=mask_val) +end function rotated_mpp_chksum_real_1d + + +!> Compute the FMS (mpp) checksum of a rotated 2d field. +function rotated_mpp_chksum_real_2d(field, pelist, mask_val, turns) & + result(chksum) + real, intent(in) :: field(:,:) !> Unrotated input field + integer, optional, intent(in) :: pelist(:) !> PE list of ranks to checksum + real, optional, intent(in) :: mask_val !> FMS mask value + integer, optional, intent(in) :: turns !> Number of quarter-turns + integer :: chksum !> FMS checksum of field + + real, allocatable :: field_rot(:,:) + integer :: qturns + + qturns = 0 + if (present(turns)) & + qturns = turns + + if (qturns == 0) then + chksum = mpp_chksum(field, pelist=pelist, mask_val=mask_val) + else + call allocate_rotated_array(field, [1,1], qturns, field_rot) + call rotate_array(field, qturns, field_rot) + chksum = mpp_chksum(field_rot, pelist=pelist, mask_val=mask_val) + deallocate(field_rot) + endif +end function rotated_mpp_chksum_real_2d + + +!> Compute the FMS (mpp) checksum of a rotated 3d field. +function rotated_mpp_chksum_real_3d(field, pelist, mask_val, turns) & + result(chksum) + real, intent(in) :: field(:,:,:) !> Unrotated input field + integer, optional, intent(in) :: pelist(:) !> PE list of ranks to checksum + real, optional, intent(in) :: mask_val !> FMS mask value + integer, optional, intent(in) :: turns !> Number of quarter-turns + integer :: chksum !> FMS checksum of field + + real, allocatable :: field_rot(:,:,:) + integer :: qturns + + qturns = 0 + if (present(turns)) & + qturns = turns + + if (qturns == 0) then + chksum = mpp_chksum(field, pelist=pelist, mask_val=mask_val) + else + call allocate_rotated_array(field, [1,1,1], qturns, field_rot) + call rotate_array(field, qturns, field_rot) + chksum = mpp_chksum(field_rot, pelist=pelist, mask_val=mask_val) + deallocate(field_rot) + endif +end function rotated_mpp_chksum_real_3d + + +!> Compute the FMS (mpp) checksum of a rotated 4d field. +function rotated_mpp_chksum_real_4d(field, pelist, mask_val, turns) & + result(chksum) + real, intent(in) :: field(:,:,:,:) !> Unrotated input field + integer, optional, intent(in) :: pelist(:) !> PE list of ranks to checksum + real, optional, intent(in) :: mask_val !> FMS mask value + integer, optional, intent(in) :: turns !> Number of quarter-turns + integer :: chksum !> FMS checksum of field + + real, allocatable :: field_rot(:,:,:,:) + integer :: qturns + + qturns = 0 + if (present(turns)) & + qturns = turns + + if (qturns == 0) then + chksum = mpp_chksum(field, pelist=pelist, mask_val=mask_val) + else + call allocate_rotated_array(field, [1,1,1,1], qturns, field_rot) + call rotate_array(field, qturns, field_rot) + chksum = mpp_chksum(field_rot, pelist=pelist, mask_val=mask_val) + deallocate(field_rot) + endif +end function rotated_mpp_chksum_real_4d + + +! NOTE: In MOM_io, write_field points to mpp_write, which supports a very broad +! range of interfaces. Here, we only support the much more narrow family of +! mpp_write_2ddecomp functions used to write tiled data. + + +!> Write the rotation of a 1d field to an FMS output file +!! This function is provided to support the full FMS write_field interface. +subroutine rotated_write_field_real_0d(io_unit, field_md, field, tstamp, turns) + integer, intent(in) :: io_unit !> File I/O unit handle + type(fieldtype), intent(in) :: field_md !> FMS field metadata + real, intent(inout) :: field !> Unrotated field array + real, optional, intent(in) :: tstamp !> Model timestamp + integer, optional, intent(in) :: turns !> Number of quarter-turns + + if (present(turns)) & + call MOM_error(FATAL, "Rotation not supported for 0d fields.") + + call write_field(io_unit, field_md, field, tstamp=tstamp) +end subroutine rotated_write_field_real_0d + + +!> Write the rotation of a 1d field to an FMS output file +!! This function is provided to support the full FMS write_field interface. +subroutine rotated_write_field_real_1d(io_unit, field_md, field, tstamp, turns) + integer, intent(in) :: io_unit !> File I/O unit handle + type(fieldtype), intent(in) :: field_md !> FMS field metadata + real, intent(inout) :: field(:) !> Unrotated field array + real, optional, intent(in) :: tstamp !> Model timestamp + integer, optional, intent(in) :: turns !> Number of quarter-turns + + if (present(turns)) & + call MOM_error(FATAL, "Rotation not supported for 0d fields.") + + call write_field(io_unit, field_md, field, tstamp=tstamp) +end subroutine rotated_write_field_real_1d + + +!> Write the rotation of a 2d field to an FMS output file +subroutine rotated_write_field_real_2d(io_unit, field_md, domain, field, & + tstamp, tile_count, default_data, turns) + integer, intent(in) :: io_unit !> File I/O unit handle + type(fieldtype), intent(in) :: field_md !> FMS field metadata + type(domain2D), intent(inout) :: domain !> FMS MPP domain + real, intent(inout) :: field(:,:) !> Unrotated field array + real, optional, intent(in) :: tstamp !> Model timestamp + integer, optional, intent(in) :: tile_count !> PEs per tile (default: 1) + real, optional, intent(in) :: default_data !> Default fill value + integer, optional, intent(in) :: turns !> Number of quarter-turns + + real, allocatable :: field_rot(:,:) + integer :: qturns + + qturns = 0 + if (present(turns)) & + qturns = turns + + if (qturns == 0) then + call write_field(io_unit, field_md, domain, field, tstamp=tstamp, & + tile_count=tile_count, default_data=default_data) + else + call allocate_rotated_array(field, [1,1], qturns, field_rot) + call rotate_array(field, qturns, field_rot) + call write_field(io_unit, field_md, domain, field_rot, tstamp=tstamp, & + tile_count=tile_count, default_data=default_data) + deallocate(field_rot) + endif +end subroutine rotated_write_field_real_2d + + +!> Write the rotation of a 3d field to an FMS output file +subroutine rotated_write_field_real_3d(io_unit, field_md, domain, field, & + tstamp, tile_count, default_data, turns) + integer, intent(in) :: io_unit !> File I/O unit handle + type(fieldtype), intent(in) :: field_md !> FMS field metadata + type(domain2D), intent(inout) :: domain !> FMS MPP domain + real, intent(inout) :: field(:,:,:) !> Unrotated field array + real, optional, intent(in) :: tstamp !> Model timestamp + integer, optional, intent(in) :: tile_count !> PEs per tile (default: 1) + real, optional, intent(in) :: default_data !> Default fill value + integer, optional, intent(in) :: turns !> Number of quarter-turns + + real, allocatable :: field_rot(:,:,:) + integer :: qturns + + qturns = 0 + if (present(turns)) & + qturns = turns + + if (qturns == 0) then + call write_field(io_unit, field_md, domain, field, tstamp=tstamp, & + tile_count=tile_count, default_data=default_data) + else + call allocate_rotated_array(field, [1,1,1], qturns, field_rot) + call rotate_array(field, qturns, field_rot) + call write_field(io_unit, field_md, domain, field_rot, tstamp=tstamp, & + tile_count=tile_count, default_data=default_data) + deallocate(field_rot) + endif +end subroutine rotated_write_field_real_3d + + +!> Write the rotation of a 4d field to an FMS output file +subroutine rotated_write_field_real_4d(io_unit, field_md, domain, field, & + tstamp, tile_count, default_data, turns) + integer, intent(in) :: io_unit !> File I/O unit handle + type(fieldtype), intent(in) :: field_md !> FMS field metadata + type(domain2D), intent(inout) :: domain !> FMS MPP domain + real, intent(inout) :: field(:,:,:,:) !> Unrotated field array + real, optional, intent(in) :: tstamp !> Model timestamp + integer, optional, intent(in) :: tile_count !> PEs per tile (default: 1) + real, optional, intent(in) :: default_data !> Default fill value + integer, optional, intent(in) :: turns !> Number of quarter-turns + + real, allocatable :: field_rot(:,:,:,:) + integer :: qturns + + qturns = 0 + if (present(turns)) qturns = turns + + if (qturns == 0) then + call write_field(io_unit, field_md, domain, field, tstamp=tstamp, & + tile_count=tile_count, default_data=default_data) + else + call allocate_rotated_array(field, [1,1,1,1], qturns, field_rot) + call rotate_array(field, qturns, field_rot) + call write_field(io_unit, field_md, domain, field_rot, tstamp=tstamp, & + tile_count=tile_count, default_data=default_data) + deallocate(field_rot) + endif +end subroutine rotated_write_field_real_4d + + +!> Read a scalar field based on model time +!! This function is provided to support the full FMS time_interp_external +!! interface. +subroutine rotated_time_interp_external_0d(fms_id, time, data_in, verbose, & + turns) + integer, intent(in) :: fms_id !< FMS field ID + type(time_type), intent(in) :: time !< Model time + real, intent(inout) :: data_in !< field to write data + logical, intent(in), optional :: verbose !< Verbose output + integer, intent(in), optional :: turns !< Number of quarter turns + + if (present(turns)) & + call MOM_error(FATAL, "Rotation not supported for 0d fields.") + + call time_interp_external(fms_id, time, data_in, verbose=verbose) +end subroutine rotated_time_interp_external_0d + +!> Read a 2d field based on model time, and rotate to the model grid +subroutine rotated_time_interp_external_2d(fms_id, time, data_in, interp, & + verbose, horz_interp, mask_out, is_in, ie_in, js_in, je_in, window_id, & + turns) + integer, intent(in) :: fms_id + type(time_type), intent(in) :: time + real, dimension(:,:), intent(inout) :: data_in + integer, intent(in), optional :: interp + logical, intent(in), optional :: verbose + type(horiz_interp_type),intent(in), optional :: horz_interp + logical, dimension(:,:), intent(out), optional :: mask_out + integer, intent(in), optional :: is_in, ie_in, js_in, je_in + integer, intent(in), optional :: window_id + integer, intent(in), optional :: turns + + real, allocatable :: data_pre(:,:) + integer :: qturns + + ! TODO: Mask rotation requires logical array rotation support + if (present(mask_out)) & + call MOM_error(FATAL, "Rotation of masked output not yet support") + + qturns = 0 + if (present(turns)) qturns = turns + + if (qturns == 0) then + call time_interp_external(fms_id, time, data_in, interp=interp, & + verbose=verbose, horz_interp=horz_interp, mask_out=mask_out, & + is_in=is_in, ie_in=ie_in, js_in=js_in, je_in=je_in, & + window_id=window_id) + else + call allocate_rotated_array(data_in, [1,1], -qturns, data_pre) + call time_interp_external(fms_id, time, data_pre, interp=interp, & + verbose=verbose, horz_interp=horz_interp, mask_out=mask_out, & + is_in=is_in, ie_in=ie_in, js_in=js_in, je_in=je_in, & + window_id=window_id) + call rotate_array(data_pre, turns, data_in) + endif +end subroutine rotated_time_interp_external_2d + + +!> Read a 3d field based on model time, and rotate to the model grid +subroutine rotated_time_interp_external_3d(fms_id, time, data_in, interp, & + verbose, horz_interp, mask_out, is_in, ie_in, js_in, je_in, window_id, & + turns) + integer, intent(in) :: fms_id + type(time_type), intent(in) :: time + real, dimension(:,:,:), intent(inout) :: data_in + integer, intent(in), optional :: interp + logical, intent(in), optional :: verbose + type(horiz_interp_type),intent(in), optional :: horz_interp + logical, dimension(:,:,:), intent(out), optional :: mask_out + integer, intent(in), optional :: is_in, ie_in, js_in, je_in + integer, intent(in), optional :: window_id + integer, intent(in), optional :: turns + + real, allocatable :: data_pre(:,:,:) + integer :: qturns + + ! TODO: Mask rotation requires logical array rotation support + if (present(mask_out)) & + call MOM_error(FATAL, "Rotation of masked output not yet support") + + qturns = 0 + if (present(turns)) qturns = turns + + if (qturns == 0) then + call time_interp_external(fms_id, time, data_in, interp=interp, & + verbose=verbose, horz_interp=horz_interp, mask_out=mask_out, & + is_in=is_in, ie_in=ie_in, js_in=js_in, je_in=je_in, & + window_id=window_id) + else + call allocate_rotated_array(data_in, [1,1,1], -qturns, data_pre) + call time_interp_external(fms_id, time, data_pre, interp=interp, & + verbose=verbose, horz_interp=horz_interp, mask_out=mask_out, & + is_in=is_in, ie_in=ie_in, js_in=js_in, je_in=je_in, & + window_id=window_id) + call rotate_array(data_pre, turns, data_in) + endif +end subroutine rotated_time_interp_external_3d + +end module MOM_transform_FMS diff --git a/src/initialization/MOM_grid_initialize.F90 b/src/initialization/MOM_grid_initialize.F90 index 1c594f45c1..45c903f4ff 100644 --- a/src/initialization/MOM_grid_initialize.F90 +++ b/src/initialization/MOM_grid_initialize.F90 @@ -125,7 +125,8 @@ subroutine grid_metrics_chksum(parent, G, US) halo = min(G%ied-G%iec, G%jed-G%jec, 1) - call hchksum_pair(trim(parent)//': d[xy]T', G%dxT, G%dyT, G%HI, haloshift=halo, scale=L_to_m) + call hchksum_pair(trim(parent)//': d[xy]T', G%dxT, G%dyT, G%HI, & + haloshift=halo, scale=L_to_m, scalar_pair=.true.) call uvchksum(trim(parent)//': dxC[uv]', G%dxCu, G%dyCv, G%HI, haloshift=halo, scale=L_to_m) @@ -133,7 +134,8 @@ subroutine grid_metrics_chksum(parent, G, US) call Bchksum_pair(trim(parent)//': dxB[uv]', G%dxBu, G%dyBu, G%HI, haloshift=halo, scale=L_to_m) - call hchksum_pair(trim(parent)//': Id[xy]T', G%IdxT, G%IdyT, G%HI, haloshift=halo, scale=m_to_L) + call hchksum_pair(trim(parent)//': Id[xy]T', G%IdxT, G%IdyT, G%HI, & + haloshift=halo, scale=m_to_L, scalar_pair=.true.) call uvchksum(trim(parent)//': Id[xy]C[uv]', G%IdxCu, G%IdyCv, G%HI, haloshift=halo, scale=m_to_L) diff --git a/src/initialization/MOM_shared_initialization.F90 b/src/initialization/MOM_shared_initialization.F90 index 3338f1fedb..9311003863 100644 --- a/src/initialization/MOM_shared_initialization.F90 +++ b/src/initialization/MOM_shared_initialization.F90 @@ -1202,6 +1202,8 @@ subroutine write_ocean_geometry_file(G, param_file, directory, geom_file, US) real, dimension(G%IsdB:G%IedB,G%jsd :G%jed ) :: out_u real, dimension(G%isd :G%ied ,G%JsdB:G%JedB) :: out_v + call callTree_enter('write_ocean_geometry_file()') + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed @@ -1331,6 +1333,7 @@ subroutine write_ocean_geometry_file(G, param_file, directory, geom_file, US) call close_file(unit) + call callTree_leave('write_ocean_geometry_file()') end subroutine write_ocean_geometry_file end module MOM_shared_initialization diff --git a/src/parameterizations/lateral/MOM_MEKE.F90 b/src/parameterizations/lateral/MOM_MEKE.F90 index 4033d64f3c..eedd9e9268 100644 --- a/src/parameterizations/lateral/MOM_MEKE.F90 +++ b/src/parameterizations/lateral/MOM_MEKE.F90 @@ -196,8 +196,10 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h if (associated(MEKE%GM_src)) & call hchksum(MEKE%GM_src, 'MEKE GM_src', G%HI, scale=US%RZ3_T3_to_W_m2*US%L_to_Z**2) if (associated(MEKE%MEKE)) call hchksum(MEKE%MEKE, 'MEKE MEKE', G%HI, scale=US%L_T_to_m_s**2) - call uvchksum("MEKE SN_[uv]", SN_u, SN_v, G%HI, scale=US%s_to_T) - call uvchksum("MEKE h[uv]", hu, hv, G%HI, haloshift=1, scale=GV%H_to_m*US%L_to_m**2) + call uvchksum("MEKE SN_[uv]", SN_u, SN_v, G%HI, scale=US%s_to_T, & + scalar_pair=.true.) + call uvchksum("MEKE h[uv]", hu, hv, G%HI, haloshift=1, & + scale=GV%H_to_m*(US%L_to_m**2)) endif sdt = dt*CS%MEKE_dtScale ! Scaled dt to use for time-stepping @@ -287,7 +289,8 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h call MEKE_lengthScales(CS, MEKE, G, GV, US, SN_u, SN_v, MEKE%MEKE, bottomFac2, barotrFac2, LmixScale) if (CS%debug) then if (CS%visc_drag) & - call uvchksum("MEKE drag_vel_[uv]", drag_vel_u, drag_vel_v, G%HI, scale=US%Z_to_m*US%s_to_T) + call uvchksum("MEKE drag_vel_[uv]", drag_vel_u, drag_vel_v, G%HI, & + scale=US%Z_to_m*US%s_to_T, scalar_pair=.true.) call hchksum(mass, 'MEKE mass',G%HI,haloshift=1, scale=US%RZ_to_kg_m2) call hchksum(drag_rate_visc, 'MEKE drag_rate_visc', G%HI, scale=US%L_T_to_m_s) call hchksum(bottomFac2, 'MEKE bottomFac2', G%HI) diff --git a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 index 0e237fac55..e5e699ebee 100644 --- a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 +++ b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 @@ -602,9 +602,12 @@ subroutine calc_Visbeck_coeffs(h, slope_x, slope_y, N2_u, N2_v, G, GV, US, CS) endif if (CS%debug) then - call uvchksum("calc_Visbeck_coeffs slope_[xy]", slope_x, slope_y, G%HI, haloshift=1) - call uvchksum("calc_Visbeck_coeffs N2_u, N2_v", N2_u, N2_v, G%HI, scale=US%s_to_T**2) - call uvchksum("calc_Visbeck_coeffs SN_[uv]", CS%SN_u, CS%SN_v, G%HI, scale=US%s_to_T) + call uvchksum("calc_Visbeck_coeffs slope_[xy]", slope_x, slope_y, G%HI, & + haloshift=1) + call uvchksum("calc_Visbeck_coeffs N2_u, N2_v", N2_u, N2_v, G%HI, & + scale=US%s_to_T**2, scalar_pair=.true.) + call uvchksum("calc_Visbeck_coeffs SN_[uv]", CS%SN_u, CS%SN_v, G%HI, & + scale=US%s_to_T, scalar_pair=.true.) endif end subroutine calc_Visbeck_coeffs diff --git a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 index 4da62ed5df..6796da5b57 100644 --- a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 +++ b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 @@ -422,7 +422,8 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp endif if (CS%debug) then - call uvchksum("Kh_[uv]", Kh_u, Kh_v, G%HI, haloshift=0, scale=US%L_to_m**2*US%s_to_T) + call uvchksum("Kh_[uv]", Kh_u, Kh_v, G%HI, haloshift=0, & + scale=(US%L_to_m**2)*US%s_to_T, scalar_pair=.true.) call uvchksum("int_slope_[uv]", int_slope_u, int_slope_v, G%HI, haloshift=0) call hchksum(h, "thickness_diffuse_1 h", G%HI, haloshift=1, scale=GV%H_to_m) call hchksum(e, "thickness_diffuse_1 e", G%HI, haloshift=1, scale=US%Z_to_m) diff --git a/src/parameterizations/vertical/MOM_ALE_sponge.F90 b/src/parameterizations/vertical/MOM_ALE_sponge.F90 index ccd85280f5..b791535ed1 100644 --- a/src/parameterizations/vertical/MOM_ALE_sponge.F90 +++ b/src/parameterizations/vertical/MOM_ALE_sponge.F90 @@ -12,6 +12,7 @@ module MOM_ALE_sponge ! This file is part of MOM6. See LICENSE.md for the license. +use MOM_array_transform, only: rotate_array use MOM_coms, only : sum_across_PEs use MOM_diag_mediator, only : post_data, query_averaging_enabled, register_diag_field use MOM_diag_mediator, only : diag_ctrl @@ -54,6 +55,7 @@ module MOM_ALE_sponge public set_up_ALE_sponge_field, set_up_ALE_sponge_vel_field public get_ALE_sponge_thicknesses, get_ALE_sponge_nz_data public initialize_ALE_sponge, apply_ALE_sponge, ALE_sponge_end, init_ALE_sponge_diags +public rotate_ALE_sponge, update_ALE_sponge_field ! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional ! consistency testing. These are noted in comments with units like Z, H, L, and T, along with @@ -999,6 +1001,163 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) end subroutine apply_ALE_sponge +!> Rotate the ALE sponge fields from the input to the model index map. +subroutine rotate_ALE_sponge(sponge_in, G_in, sponge, G, turns, param_file) + type(ALE_sponge_CS), intent(in) :: sponge_in + type(ocean_grid_type), intent(in) :: G_in + type(ALE_sponge_CS), pointer :: sponge + type(ocean_grid_type), intent(in) :: G + integer, intent(in) :: turns + type(param_file_type), intent(in) :: param_file + + ! First part: Index construction + ! 1. Reconstruct Iresttime(:,:) from sponge_in + ! 2. rotate Iresttime(:,:) + ! 3. Call initialize_sponge using new grid and rotated Iresttime(:,:) + ! All the index adjustment should follow from the Iresttime rotation + + real, dimension(:,:), allocatable :: Iresttime_in, Iresttime + real, dimension(:,:,:), allocatable :: data_h_in, data_h + real, dimension(:,:,:), allocatable :: sp_val_in, sp_val + real, dimension(:,:,:), pointer :: sp_ptr => NULL() + integer :: c, c_i, c_j + integer :: k, nz_data + integer :: n + logical :: fixed_sponge + + fixed_sponge = .not. sponge_in%time_varying_sponges + ! NOTE: nz_data is only conditionally set when fixed_sponge is true. + + allocate(Iresttime_in(G_in%isd:G_in%ied, G_in%jsd:G_in%jed)) + allocate(Iresttime(G%isd:G%ied, G%jsd:G%jed)) + Iresttime_in(:,:) = 0.0 + + if (fixed_sponge) then + nz_data = sponge_in%nz_data + allocate(data_h_in(G_in%isd:G_in%ied, G_in%jsd:G_in%jed, nz_data)) + allocate(data_h(G%isd:G%ied, G%jsd:G%jed, nz_data)) + data_h_in(:,:,:) = 0. + endif + + ! Re-populate the 2D Iresttime and data_h arrays on the original grid + do c = 1, sponge_in%num_col + c_i = sponge_in%col_i(c) + c_j = sponge_in%col_j(c) + Iresttime_in(c_i, c_j) = sponge_in%Iresttime_col(c) + if (fixed_sponge) then + do k = 1, nz_data + data_h(c_i, c_j, k) = sponge_in%Ref_h%p(k,c) + enddo + endif + enddo + + call rotate_array(Iresttime_in, turns, Iresttime) + if (fixed_sponge) then + call rotate_array(data_h_in, turns, data_h) + call initialize_ALE_sponge_fixed(Iresttime, G, param_file, sponge, & + data_h, nz_data) + else + call initialize_ALE_sponge_varying(Iresttime, G, param_file, sponge) + endif + + deallocate(Iresttime_in) + deallocate(Iresttime) + if (fixed_sponge) then + deallocate(data_h_in) + deallocate(data_h) + endif + + ! Second part: Provide rotated fields for which relaxation is applied + + sponge%fldno = sponge_in%fldno + + if (fixed_sponge) then + allocate(sp_val_in(G_in%isd:G_in%ied, G_in%jsd:G_in%jed, nz_data)) + allocate(sp_val(G%isd:G%ied, G%jsd:G%jed, nz_data)) + endif + + do n = 1, sponge_in%fldno + ! Assume that tracers are pointers and are remapped in other functions(?) + sp_ptr => sponge_in%var(n)%p + sp_val_in(:,:,:) = 0.0 + do c = 1, sponge_in%num_col + c_i = sponge_in%col_i(c) + c_j = sponge_in%col_j(c) + if (fixed_sponge) then + do k = 1, nz_data + sp_val_in(c_i, c_j, k) = sponge_in%Ref_val(n)%p(k,c) + enddo + endif + enddo + + call rotate_array(sp_val_in, turns, sp_val) + if (fixed_sponge) then + ! NOTE: This points sp_val with the unrotated field. See note below. + call set_up_ALE_sponge_field(sp_val, G, sp_ptr, sponge) + else + ! We don't want to repeat FMS init in set_up_ALE_sponge_field_varying() + ! (time_interp_external_init, init_external_field, etc), so we manually + ! do a portion of this function below. + sponge%Ref_val(n)%id = sponge_in%Ref_val(n)%id + sponge%Ref_val(n)%num_tlevs = sponge_in%Ref_val(n)%num_tlevs + + nz_data = sponge_in%Ref_val(n)%nz_data + sponge%Ref_val(n)%nz_data = nz_data + + allocate(sponge%Ref_val(n)%p(nz_data, sponge_in%num_col)) + allocate(sponge%Ref_val(n)%h(nz_data, sponge_in%num_col)) + sponge%Ref_val(n)%p(:,:) = 0.0 + sponge%Ref_val(n)%h(:,:) = 0.0 + + ! TODO: There is currently no way to associate a generic field pointer to + ! its rotated equivalent without introducing a new data structure which + ! explicitly tracks the pairing. + ! + ! As a temporary fix, we store the pointer to the unrotated field in + ! the rotated sponge, and use this reference to replace the pointer + ! to the rotated field update_ALE_sponge field. + ! + ! This makes a lot of unverifiable assumptions, and should not be + ! considered the final solution. + sponge%var(n)%p => sp_ptr + endif + enddo + + if (fixed_sponge) then + deallocate(sp_val_in) + deallocate(sp_val) + endif + + ! TODO: var_u and var_v sponge dampling is not yet supported. + if (associated(sponge_in%var_u%p) .or. associated(sponge_in%var_v%p)) & + call MOM_error(FATAL, "Rotation of ALE sponge velocities is not yet " & + // "implemented.") + + ! Transfer any existing diag_CS reference pointer + sponge%diag => sponge_in%diag + + ! NOTE: initialize_ALE_sponge_* resolves remap_cs +end subroutine rotate_ALE_sponge + + +!> Scan the ALE sponge variables and replace a prescribed pointer to a new value. +! TODO: This function solely exists to replace field pointers in the sponge +! after rotation. This function is part of a temporary solution until +! something more robust is developed. +subroutine update_ALE_sponge_field(sponge, p_old, p_new) + type(ALE_sponge_CS), pointer :: sponge + real, dimension(:,:,:), target, intent(in) :: p_old + real, dimension(:,:,:), target, intent(in) :: p_new + + integer :: n + + do n = 1, sponge%fldno + if (associated(sponge%var(n)%p, p_old)) & + sponge%var(n)%p => p_new + enddo +end subroutine update_ALE_sponge_field + + ! GMM: I could not find where sponge_end is being called, but I am keeping ! ALE_sponge_end here so we can add that if needed. !> This subroutine deallocates any memory associated with the ALE_sponge module. diff --git a/src/parameterizations/vertical/MOM_set_diffusivity.F90 b/src/parameterizations/vertical/MOM_set_diffusivity.F90 index 3045639232..d7dfcea090 100644 --- a/src/parameterizations/vertical/MOM_set_diffusivity.F90 +++ b/src/parameterizations/vertical/MOM_set_diffusivity.F90 @@ -3,12 +3,11 @@ module MOM_set_diffusivity ! This file is part of MOM6. See LICENSE.md for the license. -use MOM_checksums, only : hchksum_pair use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end use MOM_cpu_clock, only : CLOCK_MODULE_DRIVER, CLOCK_MODULE, CLOCK_ROUTINE use MOM_diag_mediator, only : diag_ctrl, time_type use MOM_diag_mediator, only : post_data, register_diag_field -use MOM_debugging, only : hchksum, uvchksum, Bchksum +use MOM_debugging, only : hchksum, uvchksum, Bchksum, hchksum_pair use MOM_EOS, only : calculate_density, calculate_density_derivs use MOM_error_handler, only : MOM_error, is_root_pe, FATAL, WARNING, NOTE use MOM_error_handler, only : callTree_showQuery @@ -536,13 +535,15 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & endif if (associated(visc%kv_bbl_u) .and. associated(visc%kv_bbl_v)) then - call uvchksum("BBL Kv_bbl_[uv]", visc%kv_bbl_u, visc%kv_bbl_v, & - G%HI, 0, symmetric=.true., scale=US%Z2_T_to_m2_s) + call uvchksum("BBL Kv_bbl_[uv]", visc%kv_bbl_u, visc%kv_bbl_v, G%HI, & + haloshift=0, symmetric=.true., scale=US%Z2_T_to_m2_s, & + scalar_pair=.true.) endif if (associated(visc%bbl_thick_u) .and. associated(visc%bbl_thick_v)) then - call uvchksum("BBL bbl_thick_[uv]", visc%bbl_thick_u, & - visc%bbl_thick_v, G%HI, 0, symmetric=.true., scale=US%Z_to_m) + call uvchksum("BBL bbl_thick_[uv]", visc%bbl_thick_u, visc%bbl_thick_v, & + G%HI, haloshift=0, symmetric=.true., scale=US%Z_to_m, & + scalar_pair=.true.) endif if (associated(visc%Ray_u) .and. associated(visc%Ray_v)) then diff --git a/src/parameterizations/vertical/MOM_set_viscosity.F90 b/src/parameterizations/vertical/MOM_set_viscosity.F90 index f5412facae..062642c3da 100644 --- a/src/parameterizations/vertical/MOM_set_viscosity.F90 +++ b/src/parameterizations/vertical/MOM_set_viscosity.F90 @@ -943,10 +943,11 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, symmetrize) if (associated(visc%Ray_u) .and. associated(visc%Ray_v)) & call uvchksum("Ray [uv]", visc%Ray_u, visc%Ray_v, G%HI, haloshift=0, scale=US%Z_to_m) if (associated(visc%kv_bbl_u) .and. associated(visc%kv_bbl_v)) & - call uvchksum("kv_bbl_[uv]", visc%kv_bbl_u, visc%kv_bbl_v, G%HI, haloshift=0, scale=US%Z2_T_to_m2_s) + call uvchksum("kv_bbl_[uv]", visc%kv_bbl_u, visc%kv_bbl_v, G%HI, & + haloshift=0, scale=US%Z2_T_to_m2_s, scalar_pair=.true.) if (associated(visc%bbl_thick_u) .and. associated(visc%bbl_thick_v)) & - call uvchksum("bbl_thick_[uv]", visc%bbl_thick_u, & - visc%bbl_thick_v, G%HI, haloshift=0, scale=US%Z_to_m) + call uvchksum("bbl_thick_[uv]", visc%bbl_thick_u, visc%bbl_thick_v, & + G%HI, haloshift=0, scale=US%Z_to_m, scalar_pair=.true.) endif end subroutine set_viscous_BBL @@ -1710,8 +1711,8 @@ subroutine set_viscous_ML(u, v, h, tv, forces, visc, dt, G, GV, US, CS, symmetri if (CS%debug) then if (associated(visc%nkml_visc_u) .and. associated(visc%nkml_visc_v)) & - call uvchksum("nkml_visc_[uv]", visc%nkml_visc_u, & - visc%nkml_visc_v, G%HI,haloshift=0) + call uvchksum("nkml_visc_[uv]", visc%nkml_visc_u, visc%nkml_visc_v, & + G%HI, haloshift=0, scalar_pair=.true.) endif if (CS%id_nkml_visc_u > 0) & call post_data(CS%id_nkml_visc_u, visc%nkml_visc_u, CS%diag) diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90 index 5a610095ce..b4e2e302c8 100644 --- a/src/parameterizations/vertical/MOM_vert_friction.F90 +++ b/src/parameterizations/vertical/MOM_vert_friction.F90 @@ -558,7 +558,8 @@ subroutine vertvisc_remnant(visc, visc_rem_u, visc_rem_v, dt, G, GV, US, CS) enddo ! end of v-component J loop if (CS%debug) then - call uvchksum("visc_rem_[uv]", visc_rem_u, visc_rem_v, G%HI,haloshift=0) + call uvchksum("visc_rem_[uv]", visc_rem_u, visc_rem_v, G%HI, haloshift=0, & + scalar_pair=.true.) endif end subroutine vertvisc_remnant @@ -1008,10 +1009,13 @@ subroutine vertvisc_coef(u, v, h, forces, visc, dt, G, GV, US, CS, OBC) enddo ! end of v-point j loop if (CS%debug) then - call uvchksum("vertvisc_coef h_[uv]", CS%h_u, CS%h_v, G%HI, haloshift=0, scale=GV%H_to_m) - call uvchksum("vertvisc_coef a_[uv]", CS%a_u, CS%a_v, G%HI, haloshift=0, scale=US%Z_to_m*US%s_to_T) + call uvchksum("vertvisc_coef h_[uv]", CS%h_u, CS%h_v, G%HI, haloshift=0, & + scale=GV%H_to_m, scalar_pair=.true.) + call uvchksum("vertvisc_coef a_[uv]", CS%a_u, CS%a_v, G%HI, haloshift=0, & + scale=US%Z_to_m*US%s_to_T, scalar_pair=.true.) if (allocated(hML_u) .and. allocated(hML_v)) & - call uvchksum("vertvisc_coef hML_[uv]", hML_u, hML_v, G%HI, haloshift=0, scale=GV%H_to_m) + call uvchksum("vertvisc_coef hML_[uv]", hML_u, hML_v, G%HI, & + haloshift=0, scale=GV%H_to_m, scalar_pair=.true.) endif ! Offer diagnostic fields for averaging. diff --git a/src/tracer/MOM_tracer_hor_diff.F90 b/src/tracer/MOM_tracer_hor_diff.F90 index d18bb3e330..aec9c8ccf2 100644 --- a/src/tracer/MOM_tracer_hor_diff.F90 +++ b/src/tracer/MOM_tracer_hor_diff.F90 @@ -556,10 +556,12 @@ subroutine tracer_hordiff(h, dt, MEKE, VarMix, G, GV, US, CS, Reg, tv, do_online if (CS%debug) then call uvchksum("After tracer diffusion khdt_[xy]", khdt_x, khdt_y, & - G%HI, haloshift=0, symmetric=.true., scale=US%L_to_m**2) + G%HI, haloshift=0, symmetric=.true., scale=US%L_to_m**2, & + scalar_pair=.true.) if (CS%use_neutral_diffusion) then call uvchksum("After tracer diffusion Coef_[xy]", Coef_x, Coef_y, & - G%HI, haloshift=0, symmetric=.true., scale=US%L_to_m**2) + G%HI, haloshift=0, symmetric=.true., scale=US%L_to_m**2, & + scalar_pair=.true.) endif endif