From 987a0ef3e84e7b4f08f4034030bb88cefc29ba63 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Fri, 18 Nov 2022 05:51:14 -0500 Subject: [PATCH] +Add REMAP_AUXILIARY_VARS to remap accelerations Added REMAP_AUXILIARY_VARS to control whether to remap the accelerations that are used in the predictor step of the split RK2 time stepping scheme. Also added the new routines remap_dyn_split_rk2_aux_vars, remap_OBC_fields and remap_vertvisc_aux_vars to do the remapping, and code to call these routines when REMAP_AUXILIARY_VARS is true. By default, REMAP_AUXILIARY_VARS is false, and all answers are bitwise identical, but the entire MOM6-examples regression suite has been run with this set to true, and they do appear to give physically plausible answers in all cases, partially addressing the issue noted at github.com/NOAA-GFDL/MOM6/issues/203. New entries are added to the MOM_parameter_doc files, and there are three new publicly visible routines, but by default answers do not change. --- src/core/MOM.F90 | 31 ++- src/core/MOM_dynamics_split_RK2.F90 | 39 +++- src/core/MOM_open_boundary.F90 | 179 ++++++++++++++++++ .../vertical/MOM_set_viscosity.F90 | 72 +++++-- 4 files changed, 293 insertions(+), 28 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index fef71ab4d5..b20940aeba 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -74,7 +74,7 @@ module MOM use MOM_dynamics_unsplit, only : MOM_dyn_unsplit_CS use MOM_dynamics_split_RK2, only : step_MOM_dyn_split_RK2, register_restarts_dyn_split_RK2 use MOM_dynamics_split_RK2, only : initialize_dyn_split_RK2, end_dyn_split_RK2 -use MOM_dynamics_split_RK2, only : MOM_dyn_split_RK2_CS +use MOM_dynamics_split_RK2, only : MOM_dyn_split_RK2_CS, remap_dyn_split_rk2_aux_vars use MOM_dynamics_unsplit_RK2, only : step_MOM_dyn_unsplit_RK2, register_restarts_dyn_unsplit_RK2 use MOM_dynamics_unsplit_RK2, only : initialize_dyn_unsplit_RK2, end_dyn_unsplit_RK2 use MOM_dynamics_unsplit_RK2, only : MOM_dyn_unsplit_RK2_CS @@ -103,14 +103,13 @@ module MOM use MOM_mixed_layer_restrat, only : mixedlayer_restrat_register_restarts use MOM_obsolete_diagnostics, only : register_obsolete_diagnostics use MOM_open_boundary, only : ocean_OBC_type, OBC_registry_type -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 : register_temp_salt_segments, update_segment_tracer_reservoirs +use MOM_open_boundary, only : open_boundary_register_restarts, remap_OBC_fields use MOM_open_boundary, only : rotate_OBC_config, rotate_OBC_init use MOM_porous_barriers, only : porous_widths_layer, porous_widths_interface, porous_barriers_init use MOM_porous_barriers, only : porous_barrier_CS -use MOM_set_visc, only : set_viscous_BBL, set_viscous_ML -use MOM_set_visc, only : set_visc_register_restarts, set_visc_CS +use MOM_set_visc, only : set_viscous_BBL, set_viscous_ML, set_visc_CS +use MOM_set_visc, only : set_visc_register_restarts, remap_vertvisc_aux_vars use MOM_set_visc, only : set_visc_init, set_visc_end use MOM_shared_initialization, only : write_ocean_geometry_file use MOM_sponge, only : init_sponge_diags, sponge_CS @@ -250,6 +249,10 @@ module MOM logical :: use_ALE_algorithm !< If true, use the ALE algorithm rather than layered !! isopycnal/stacked shallow water mode. This logical is set by calling the !! function useRegridding() from the MOM_regridding module. + logical :: remap_aux_vars !< If true, apply ALE remapping to all of the auxiliary 3-D + !! variables that are needed to reproduce across restarts, + !! similarly to what is done with the primary state variables. + type(MOM_stoch_eos_CS) :: stoch_eos_CS !< structure containing random pattern for stoch EOS logical :: alternate_first_direction !< If true, alternate whether the x- or y-direction !! updates occur first in directionally split parts of the calculation. @@ -1547,6 +1550,16 @@ subroutine step_MOM_thermo(CS, G, GV, US, u, v, h, tv, fluxes, dtdia, & call ALE_remap_tracers(CS%ALE_CSp, G, GV, h, h_new, CS%tracer_Reg, showCallTree, dtdia, PCM_cell) call ALE_remap_velocities(CS%ALE_CSp, G, GV, h, h_new, u, v, CS%OBC, dzRegrid, showCallTree, dtdia) + if (CS%remap_aux_vars) then + if (CS%split) & + call remap_dyn_split_RK2_aux_vars(G, GV, CS%dyn_split_RK2_CSp, h, h_new, CS%ALE_CSp, CS%OBC, dzRegrid) + + if (associated(CS%OBC)) & + call remap_OBC_fields(G, GV, h, h_new, CS%OBC, PCM_cell=PCM_cell) + + call remap_vertvisc_aux_vars(G, GV, CS%visc, h, h_new, CS%ALE_CSp, CS%OBC) + endif + ! Replace the old grid with new one. All remapping must be done by this point in the code. !$OMP parallel do default(shared) do k=1,nz ; do j=js-1,je+1 ; do i=is-1,ie+1 @@ -2072,6 +2085,12 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & call get_param(param_file, "MOM", "USE_REGRIDDING", CS%use_ALE_algorithm, & "If True, use the ALE algorithm (regridding/remapping). "//& "If False, use the layered isopycnal algorithm.", default=.false. ) + call get_param(param_file, "MOM", "REMAP_AUXILIARY_VARS", CS%remap_aux_vars, & + "If true, apply ALE remapping to all of the auxiliary 3-dimensional "//& + "variables that are needed to reproduce across restarts, similarly to "//& + "what is already being done with the primary state variables. "//& + "The default should be changed to true.", default=.false., & + do_not_log=.not.CS%use_ALE_algorithm) call get_param(param_file, "MOM", "BULKMIXEDLAYER", bulkmixedlayer, & "If true, use a Kraus-Turner-like bulk mixed layer "//& "with transitional buffer layers. Layers 1 through "//& diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index 68f8c97669..f438c14a05 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -36,7 +36,7 @@ module MOM_dynamics_split_RK2 use MOM_time_manager, only : time_type, time_type_to_real, operator(+) use MOM_time_manager, only : operator(-), operator(>), operator(*), operator(/) -use MOM_ALE, only : ALE_CS +use MOM_ALE, only : ALE_CS, ALE_remap_velocities use MOM_barotropic, only : barotropic_init, btstep, btcalc, bt_mass_source use MOM_barotropic, only : register_barotropic_restarts, set_dtbt, barotropic_CS use MOM_barotropic, only : barotropic_end @@ -160,6 +160,9 @@ module MOM_dynamics_split_RK2 !! predictor step. This is used to accomodate various generations !! of restart files. logical :: use_tides !< If true, tidal forcing is enabled. + logical :: remap_aux !< If true, apply ALE remapping to all of the auxiliary 3-D + !! variables that are needed to reproduce across restarts, + !! similarly to what is done with the primary state variables. real :: be !< A nondimensional number from 0.5 to 1 that controls !! the backward weighting of the time stepping scheme [nondim] @@ -256,6 +259,7 @@ module MOM_dynamics_split_RK2 public step_MOM_dyn_split_RK2 public register_restarts_dyn_split_RK2 public initialize_dyn_split_RK2 +public remap_dyn_split_RK2_aux_vars public end_dyn_split_RK2 !>@{ CPU time clock IDs @@ -1160,6 +1164,32 @@ subroutine register_restarts_dyn_split_RK2(HI, GV, US, param_file, CS, restart_C end subroutine register_restarts_dyn_split_RK2 +!> This subroutine does remapping for the auxiliary restart variables that are used +!! with the split RK2 time stepping scheme. +subroutine remap_dyn_split_RK2_aux_vars(G, GV, CS, h_old, h_new, ALE_CSp, OBC, dzRegrid) + type(ocean_grid_type), intent(inout) :: G !< ocean grid structure + type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure + type(MOM_dyn_split_RK2_CS), pointer :: CS !< module control structure + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + intent(in) :: h_old !< Thickness of source grid [H ~> m or kg m-2] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + intent(in) :: h_new !< Thickness of destination grid [H ~> m or kg m-2] + type(ALE_CS), pointer :: ALE_CSp !< ALE control structure to use when remapping + type(ocean_OBC_type), pointer :: OBC !< OBC control structure to use when remapping + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), & + optional, intent(in) :: dzRegrid !< Change in interface position [H ~> m or kg m-2] + + if (.not.CS%remap_aux) return + + if (CS%store_CAu) then + call ALE_remap_velocities(ALE_CSp, G, GV, h_old, h_new, CS%u_av, CS%v_av, OBC, dzRegrid) + call ALE_remap_velocities(ALE_CSp, G, GV, h_old, h_new, CS%CAu_pred, CS%CAv_pred, OBC, dzRegrid) + endif + + call ALE_remap_velocities(ALE_CSp, G, GV, h_old, h_new, CS%diffu, CS%diffv, OBC, dzRegrid) + +end subroutine remap_dyn_split_RK2_aux_vars + !> This subroutine initializes all of the variables that are used by this !! dynamic core, including diagnostics and the cpu clocks. subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param_file, & @@ -1276,6 +1306,13 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param "If true, calculate the Coriolis accelerations at the end of each "//& "timestep for use in the predictor step of the next split RK2 timestep.", & default=.true.) + call get_param(param_file, mdl, "REMAP_AUXILIARY_VARS", CS%remap_aux, & + "If true, apply ALE remapping to all of the auxiliary 3-dimensional "//& + "variables that are needed to reproduce across restarts, similarly to "//& + "what is already being done with the primary state variables. "//& + "The default should be changed to true.", default=.false., do_not_log=.true.) + if (CS%remap_aux .and. .not.CS%store_CAu) call MOM_error(FATAL, & + "REMAP_AUXILIARY_VARS requires that STORE_CORIOLIS_ACCEL = True.") call get_param(param_file, mdl, "DEBUG", CS%debug, & "If true, write out verbose debugging data.", & default=.false., debuggingParam=.true.) diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index 7d89d97c0f..abe63a950a 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -65,6 +65,7 @@ module MOM_open_boundary public open_boundary_register_restarts public update_segment_tracer_reservoirs public update_OBC_ramp +public remap_OBC_fields public rotate_OBC_config public rotate_OBC_init public initialize_segment_data @@ -5408,6 +5409,184 @@ subroutine update_segment_tracer_reservoirs(G, GV, uhr, vhr, h, OBC, dt, Reg) end subroutine update_segment_tracer_reservoirs +!> Vertically remap the OBC tracer reservoirs and radiation rates that are filtered in time. +subroutine remap_OBC_fields(G, GV, h_old, h_new, OBC, PCM_cell) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_old !< Thickness of source grid [H ~> m or kg m-2] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_new !< Thickness of destination grid [H ~> m or kg m-2] + type(ocean_OBC_type), pointer :: OBC !< Open boundary structure + logical, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + optional, intent(in) :: PCM_cell !< Use PCM remapping in cells where true + + ! Local variables + type(OBC_segment_type), pointer :: segment => NULL() ! A pointer to the various segments, used just for shorthand. + + real :: tr_column(GV%ke) ! A column of updated tracer concentrations [CU ~> Conc] + real :: r_norm_col(GV%ke) ! A column of updated radiation rates, in grid points per timestep [nondim] + real :: rxy_col(GV%ke) ! A column of updated radiation rates for oblique OBCs [L2 T-2 ~> m2 s-2] + real :: h1(GV%ke) ! A column of source grid layer thicknesses [H ~> m or kg m-2] + real :: h2(GV%ke) ! A column of target grid layer thicknesses [H ~> m or kg m-2] + real :: I_scale ! The inverse of the scaling factor for the tracers. + ! For salinity the units would be [ppt S-1 ~> 1]. + real :: h_neglect ! Tiny thickness used in remapping [H ~> m or kg m-2] + logical :: PCM(GV%ke) ! If true, do PCM remapping from a cell. + integer :: i, j, k, m, n, ntr, nz + + if (.not.associated(OBC)) return + + nz = GV%ke + ntr = OBC%ntr + h_neglect = GV%H_subroundoff + + if (.not.present(PCM_cell)) PCM(:) = .false. + + if (associated(OBC)) then ; if (OBC%OBC_pe) then ; do n=1,OBC%number_of_segments + segment => OBC%segment(n) + if (.not.associated(segment%tr_Reg)) cycle + + if (segment%is_E_or_W) then + I = segment%HI%IsdB + do j=segment%HI%jsd,segment%HI%jed + + ! Store a column of the start and final grids + if (segment%direction == OBC_DIRECTION_W) then + if (G%mask2dT(i+1,j) == 0.0) cycle + h1(:) = h_old(i+1,j,:) + h2(:) = h_new(i+1,j,:) + if (present(PCM_cell)) then ; PCM(:) = PCM_cell(i+1,j,:) ; endif + else + if (G%mask2dT(i,j) == 0.0) cycle + h1(:) = h_old(i,j,:) + h2(:) = h_new(i,j,:) + if (present(PCM_cell)) then ; PCM(:) = PCM_cell(i,j,:) ; endif + endif + + ! Vertically remap the reservoir tracer concentrations + do m=1,ntr ; if (allocated(segment%tr_Reg%Tr(m)%tres)) then + I_scale = 1.0 ; if (segment%tr_Reg%Tr(m)%scale /= 0.0) I_scale = 1.0 / segment%tr_Reg%Tr(m)%scale + + if (present(PCM_cell)) then + call remapping_core_h(OBC%remap_CS, nz, h1, segment%tr_Reg%Tr(m)%tres(I,j,:), nz, h2, tr_column, & + h_neglect, h_neglect, PCM_cell=PCM) + else + call remapping_core_h(OBC%remap_CS, nz, h1, segment%tr_Reg%Tr(m)%tres(I,j,:), nz, h2, tr_column, & + h_neglect, h_neglect) + endif + + ! Possibly underflow any very tiny tracer concentrations to 0? + + ! Update tracer concentrations + segment%tr_Reg%Tr(m)%tres(I,j,:) = tr_column(:) + if (allocated(OBC%tres_x)) then ; do k=1,nz + OBC%tres_x(I,j,k,m) = I_scale * segment%tr_Reg%Tr(m)%tres(I,j,k) + enddo ; endif + + endif ; enddo + + if (segment%radiation .and. (OBC%gamma_uv < 1.0)) then + call remapping_core_h(OBC%remap_CS, nz, h1, segment%rx_norm_rad(I,j,:), nz, h2, r_norm_col, & + h_neglect, h_neglect, PCM_cell=PCM) + + do k=1,nz + segment%rx_norm_rad(I,j,k) = r_norm_col(k) + OBC%rx_normal(I,j,k) = segment%rx_norm_rad(I,j,k) + enddo + endif + + if (segment%oblique .and. (OBC%gamma_uv < 1.0)) then + call remapping_core_h(OBC%remap_CS, nz, h1, segment%rx_norm_obl(I,j,:), nz, h2, rxy_col, & + h_neglect, h_neglect, PCM_cell=PCM) + segment%rx_norm_obl(I,j,:) = rxy_col(:) + call remapping_core_h(OBC%remap_CS, nz, h1, segment%ry_norm_obl(I,j,:), nz, h2, rxy_col, & + h_neglect, h_neglect, PCM_cell=PCM) + segment%ry_norm_obl(I,j,:) = rxy_col(:) + call remapping_core_h(OBC%remap_CS, nz, h1, segment%cff_normal(I,j,:), nz, h2, rxy_col, & + h_neglect, h_neglect, PCM_cell=PCM) + segment%cff_normal(I,j,:) = rxy_col(:) + + do k=1,nz + OBC%rx_oblique_u(I,j,k) = segment%rx_norm_obl(I,j,k) + OBC%ry_oblique_u(I,j,k) = segment%ry_norm_obl(I,j,k) + OBC%cff_normal_u(I,j,k) = segment%cff_normal(I,j,k) + enddo + endif + + enddo + elseif (segment%is_N_or_S) then + J = segment%HI%JsdB + do i=segment%HI%isd,segment%HI%ied + + ! Store a column of the start and final grids + if (segment%direction == OBC_DIRECTION_S) then + if (G%mask2dT(i,j+1) == 0.0) cycle + h1(:) = h_old(i,j+1,:) + h2(:) = h_new(i,j+1,:) + if (present(PCM_cell)) then ; PCM(:) = PCM_cell(i,j+1,:) ; endif + else + if (G%mask2dT(i,j) == 0.0) cycle + h1(:) = h_old(i,j,:) + h2(:) = h_new(i,j,:) + if (present(PCM_cell)) then ; PCM(:) = PCM_cell(i,j,:) ; endif + endif + + ! Vertically remap the reservoir tracer concentrations + do m=1,ntr ; if (allocated(segment%tr_Reg%Tr(m)%tres)) then + I_scale = 1.0 ; if (segment%tr_Reg%Tr(m)%scale /= 0.0) I_scale = 1.0 / segment%tr_Reg%Tr(m)%scale + + if (present(PCM_cell)) then + call remapping_core_h(OBC%remap_CS, nz, h1, segment%tr_Reg%Tr(m)%tres(i,J,:), nz, h2, tr_column, & + h_neglect, h_neglect, PCM_cell=PCM) + else + call remapping_core_h(OBC%remap_CS, nz, h1, segment%tr_Reg%Tr(m)%tres(i,J,:), nz, h2, tr_column, & + h_neglect, h_neglect) + endif + + ! Possibly underflow any very tiny tracer concentrations to 0? + + ! Update tracer concentrations + segment%tr_Reg%Tr(m)%tres(i,J,:) = tr_column(:) + if (allocated(OBC%tres_y)) then ; do k=1,nz + OBC%tres_y(i,J,k,m) = I_scale * segment%tr_Reg%Tr(m)%tres(i,J,k) + enddo ; endif + + endif ; enddo + + if (segment%radiation .and. (OBC%gamma_uv < 1.0)) then + call remapping_core_h(OBC%remap_CS, nz, h1, segment%ry_norm_rad(i,J,:), nz, h2, r_norm_col, & + h_neglect, h_neglect, PCM_cell=PCM) + + do k=1,nz + segment%ry_norm_rad(i,J,k) = r_norm_col(k) + OBC%ry_normal(i,J,k) = segment%ry_norm_rad(i,J,k) + enddo + endif + + if (segment%oblique .and. (OBC%gamma_uv < 1.0)) then + call remapping_core_h(OBC%remap_CS, nz, h1, segment%rx_norm_obl(i,J,:), nz, h2, rxy_col, & + h_neglect, h_neglect, PCM_cell=PCM) + segment%rx_norm_obl(i,J,:) = rxy_col(:) + call remapping_core_h(OBC%remap_CS, nz, h1, segment%ry_norm_obl(i,J,:), nz, h2, rxy_col, & + h_neglect, h_neglect, PCM_cell=PCM) + segment%ry_norm_obl(i,J,:) = rxy_col(:) + call remapping_core_h(OBC%remap_CS, nz, h1, segment%cff_normal(i,J,:), nz, h2, rxy_col, & + h_neglect, h_neglect, PCM_cell=PCM) + segment%cff_normal(i,J,:) = rxy_col(:) + + do k=1,nz + OBC%rx_oblique_v(i,J,k) = segment%rx_norm_obl(i,J,k) + OBC%ry_oblique_v(i,J,k) = segment%ry_norm_obl(i,J,k) + OBC%cff_normal_v(i,J,k) = segment%cff_normal(i,J,k) + enddo + endif + + enddo + endif + enddo ; endif ; endif + +end subroutine remap_OBC_fields + + !> Adjust interface heights to fit the bathymetry and diagnose layer thickness. !! !! If the bottom most interface is below the topography then the bottom-most diff --git a/src/parameterizations/vertical/MOM_set_viscosity.F90 b/src/parameterizations/vertical/MOM_set_viscosity.F90 index 1d2bbbb048..e9497d0e92 100644 --- a/src/parameterizations/vertical/MOM_set_viscosity.F90 +++ b/src/parameterizations/vertical/MOM_set_viscosity.F90 @@ -4,37 +4,39 @@ module MOM_set_visc ! This file is part of MOM6. See LICENSE.md for the license. -use MOM_debugging, only : uvchksum, hchksum -use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, CLOCK_ROUTINE +use MOM_ALE, only : ALE_CS, ALE_remap_velocities, ALE_remap_interface_vals, ALE_remap_vertex_vals +use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, CLOCK_ROUTINE +use MOM_debugging, only : uvchksum, hchksum use MOM_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr use MOM_diag_mediator, only : diag_ctrl, time_type -use MOM_domains, only : pass_var, CORNER +use MOM_domains, only : pass_var, CORNER use MOM_error_handler, only : MOM_error, FATAL, WARNING -use MOM_file_parser, only : get_param, log_param, log_version, param_file_type -use MOM_forcing_type, only : forcing, mech_forcing -use MOM_grid, only : ocean_grid_type -use MOM_hor_index, only : hor_index_type -use MOM_io, only : slasher, MOM_read_data -use MOM_kappa_shear, only : kappa_shear_is_used, kappa_shear_at_vertex -use MOM_cvmix_shear, only : cvmix_shear_is_used -use MOM_cvmix_conv, only : cvmix_conv_is_used -use MOM_CVMix_ddiff, only : CVMix_ddiff_is_used -use MOM_restart, only : register_restart_field, query_initialized, MOM_restart_CS -use MOM_restart, only : register_restart_field_as_obsolete -use MOM_safe_alloc, only : safe_alloc_ptr, safe_alloc_alloc -use MOM_unit_scaling, only : unit_scale_type -use MOM_variables, only : thermo_var_ptrs, vertvisc_type, porous_barrier_type -use MOM_verticalGrid, only : verticalGrid_type -use MOM_EOS, only : calculate_density, calculate_density_derivs +use MOM_file_parser, only : get_param, log_param, log_version, param_file_type +use MOM_forcing_type, only : forcing, mech_forcing +use MOM_grid, only : ocean_grid_type +use MOM_hor_index, only : hor_index_type +use MOM_io, only : slasher, MOM_read_data +use MOM_kappa_shear, only : kappa_shear_is_used, kappa_shear_at_vertex +use MOM_cvmix_shear, only : cvmix_shear_is_used +use MOM_cvmix_conv, only : cvmix_conv_is_used +use MOM_CVMix_ddiff, only : CVMix_ddiff_is_used +use MOM_restart, only : register_restart_field, query_initialized, MOM_restart_CS +use MOM_restart, only : register_restart_field_as_obsolete +use MOM_safe_alloc, only : safe_alloc_ptr, safe_alloc_alloc +use MOM_unit_scaling, only : unit_scale_type +use MOM_variables, only : thermo_var_ptrs, vertvisc_type, porous_barrier_type +use MOM_verticalGrid, only : verticalGrid_type +use MOM_EOS, only : calculate_density, calculate_density_derivs use MOM_open_boundary, only : ocean_OBC_type, OBC_NONE, OBC_DIRECTION_E use MOM_open_boundary, only : OBC_DIRECTION_W, OBC_DIRECTION_N, OBC_DIRECTION_S use MOM_open_boundary, only : OBC_segment_type + implicit none ; private #include public set_viscous_BBL, set_viscous_ML, set_visc_init, set_visc_end -public set_visc_register_restarts +public set_visc_register_restarts, remap_vertvisc_aux_vars ! 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 @@ -1942,6 +1944,34 @@ subroutine set_visc_register_restarts(HI, GV, US, param_file, visc, restart_CS) end subroutine set_visc_register_restarts +!> This subroutine does remapping for the auxiliary restart variables in a vertvisc_type +!! that are used across timesteps +subroutine remap_vertvisc_aux_vars(G, GV, visc, h_old, h_new, ALE_CSp, OBC) + type(ocean_grid_type), intent(inout) :: G !< ocean grid structure + type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure + type(vertvisc_type), intent(inout) :: visc !< A structure containing vertical + !! viscosities and related fields. + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + intent(in) :: h_old !< Thickness of source grid [H ~> m or kg m-2] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + intent(in) :: h_new !< Thickness of destination grid [H ~> m or kg m-2] + type(ALE_CS), pointer :: ALE_CSp !< ALE control structure to use when remapping + type(ocean_OBC_type), pointer :: OBC !< Open boundary structure + + if (associated(visc%Kd_shear)) then + call ALE_remap_interface_vals(ALE_CSp, G, GV, h_old, h_new, visc%Kd_shear) + endif + + if (associated(visc%Kv_shear)) then + call ALE_remap_interface_vals(ALE_CSp, G, GV, h_old, h_new, visc%Kv_shear) + endif + + if (associated(visc%Kv_shear_Bu)) then + call ALE_remap_vertex_vals(ALE_CSp, G, GV, h_old, h_new, visc%Kv_shear_Bu) + endif + +end subroutine remap_vertvisc_aux_vars + !> Initializes the MOM_set_visc control structure subroutine set_visc_init(Time, G, GV, US, param_file, diag, visc, CS, restart_CS, OBC) type(time_type), target, intent(in) :: Time !< The current model time. @@ -1953,7 +1983,7 @@ subroutine set_visc_init(Time, G, GV, US, param_file, diag, visc, CS, restart_CS type(diag_ctrl), target, intent(inout) :: diag !< A structure that is used to regulate diagnostic !! output. type(vertvisc_type), intent(inout) :: visc !< A structure containing vertical viscosities and - !! related fields. Allocated here. + !! related fields. type(set_visc_CS), intent(inout) :: CS !< Vertical viscosity control structure type(MOM_restart_CS), intent(inout) :: restart_CS !< MOM restart control structure type(ocean_OBC_type), pointer :: OBC !< A pointer to an open boundary condition structure