Skip to content

Commit

Permalink
+Add REMAP_AUXILIARY_VARS to remap accelerations
Browse files Browse the repository at this point in the history 
  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/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.
  • Loading branch information
@Hallberg-NOAA
Hallberg-NOAA committed Nov 18, 2022
1 parent ba145bf commit 987a0ef
Show file tree
Hide file tree
Showing 4 changed files with 293 additions and 28 deletions.
31 changes: 25 additions & 6 deletions src/core/MOM.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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.
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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 "//&
Expand Down
39 changes: 38 additions & 1 deletion src/core/MOM_dynamics_split_RK2.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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]
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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, &
Expand Down Expand Up @@ -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.)
Expand Down
179 changes: 179 additions & 0 deletions src/core/MOM_open_boundary.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
Loading

0 comments on commit 987a0ef

Please sign in to comment.