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