Skip to content

Commit

Permalink
+Add Hybgen regridding
Browse files Browse the repository at this point in the history 
  This commit adds the ability to use Hybgen regridding, derived from Hycom
code.  There are two new files, MOM_hybgen_regrid.F90 and MOM_hybgen_unmix.F90,
which in turn add 15 new runtime parameters (with names like HYBGEN_...) to
control the hybgen code.  This new regridding option is specified via
REGRIDDING_COORDINATE_MODE="HYBGEN", and this new option is listed in the
MOM_parameter_doc files for cases that have USE_REGRIDDING=True.  There is also
a new publicly visible parameter, REGRIDDING_HYBGEN, in regrid_consts.F90.

  In addition, the new routine regridding_preadjust_reqs is
provided in MOM_regridding to specify whether convective adjustment or
hybgen_unmixing should be done before regridding.  This is used along with the
optional argument conv_adjust to regridding_main to move these pre-adjustment
steps out of regridding_main.

  The unused routine ALE_build_grid was eliminated, and comments were added
describing a few undocumented internal real variables.

  All answers are bitwise identical, but there are several new public interfaces
and there will be new lines in the comments in some MOM_parameter_doc files.
  • Loading branch information
@Hallberg-NOAA @marshallward
Hallberg-NOAA authored and marshallward committed Apr 9, 2022
1 parent 78c574b commit fd2060c
Show file tree
Hide file tree
Showing 6 changed files with 1,646 additions and 83 deletions.
103 changes: 52 additions & 51 deletions src/ALE/MOM_ALE.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,13 +20,16 @@ module MOM_ALE
use MOM_error_handler, only : MOM_error, FATAL, WARNING
use MOM_error_handler, only : callTree_showQuery
use MOM_error_handler, only : callTree_enter, callTree_leave, callTree_waypoint
use MOM_hybgen_unmix, only : hybgen_unmix, init_hybgen_unmix, end_hybgen_unmix, hybgen_unmix_CS
use MOM_hybgen_regrid, only : hybgen_regrid_CS
use MOM_file_parser, only : get_param, param_file_type, log_param
use MOM_interface_heights,only : find_eta
use MOM_open_boundary, only : ocean_OBC_type, OBC_DIRECTION_E, OBC_DIRECTION_W
use MOM_open_boundary, only : OBC_DIRECTION_N, OBC_DIRECTION_S
use MOM_regridding, only : initialize_regridding, regridding_main, end_regridding
use MOM_regridding, only : uniformResolution
use MOM_regridding, only : inflate_vanished_layers_old
use MOM_regridding, only : regridding_preadjust_reqs, convective_adjustment
use MOM_regridding, only : set_target_densities_from_GV, set_target_densities
use MOM_regridding, only : regriddingCoordinateModeDoc, DEFAULT_COORDINATE_MODE
use MOM_regridding, only : regriddingInterpSchemeDoc, regriddingDefaultInterpScheme
Expand Down Expand Up @@ -73,6 +76,11 @@ module MOM_ALE
type(remapping_CS) :: remapCS !< Remapping parameters and work arrays
type(remapping_CS) :: vel_remapCS !< Remapping parameters for velocities and work arrays

type(hybgen_unmix_CS), pointer :: hybgen_unmixCS => NULL() !< Parameters for hybgen remapping

logical :: use_hybgen_unmix !< If true, use the hybgen unmixing code before regridding
logical :: do_conv_adj !< If true, do convective adjustment before regridding

integer :: nk !< Used only for queries, not directly by this module
real :: BBL_h_vel_mask !< The thickness of a bottom boundary layer within which velocities in
!! thin layers are zeroed out after remapping, following practice with
Expand Down Expand Up @@ -119,7 +127,6 @@ module MOM_ALE
public ALE_main_offline
public ALE_offline_inputs
public ALE_offline_tracer_final
public ALE_build_grid
public ALE_regrid_accelerated
public ALE_remap_scalar
public ALE_PLM_edge_values
Expand Down Expand Up @@ -165,6 +172,8 @@ subroutine ALE_init( param_file, GV, US, max_depth, CS)
logical :: force_bounds_in_subcell
logical :: local_logical
logical :: remap_boundary_extrap
type(hybgen_regrid_CS), pointer :: hybgen_regridCS => NULL() ! Control structure for hybgen regridding
! for sharing parameters.

if (associated(CS)) then
call MOM_error(WARNING, "ALE_init called with an associated "// &
Expand All @@ -184,6 +193,7 @@ subroutine ALE_init( param_file, GV, US, max_depth, CS)

! Initialize and configure regridding
call ALE_initRegridding(GV, US, max_depth, param_file, mdl, CS%regridCS)
call regridding_preadjust_reqs(CS%regridCS, CS%do_conv_adj, CS%use_hybgen_unmix, hybgen_CS=hybgen_regridCS)

! Initialize and configure remapping that is orchestrated by ALE.
call get_param(param_file, mdl, "REMAPPING_SCHEME", string, &
Expand Down Expand Up @@ -277,6 +287,9 @@ subroutine ALE_init( param_file, GV, US, max_depth, CS)
"to avoid such filtering altogether.", &
default=1.0e-6, units="m", scale=GV%m_to_H, do_not_log=(CS%BBL_h_vel_mask<=0.0))

if (CS%use_hybgen_unmix) &
call init_hybgen_unmix(CS%hybgen_unmixCS, GV, US, param_file, hybgen_regridCS)

call get_param(param_file, "MOM", "DEBUG", CS%debug, &
"If true, write out verbose debugging data.", &
default=.false., debuggingParam=.true.)
Expand Down Expand Up @@ -350,6 +363,7 @@ subroutine ALE_end(CS)
! Deallocate memory used for the regridding
call end_remapping( CS%remapCS )

if (CS%use_hybgen_unmix) call end_hybgen_unmix( CS%hybgen_unmixCS )
call end_regridding( CS%regridCS )

deallocate(CS)
Expand Down Expand Up @@ -379,9 +393,9 @@ subroutine ALE_main( G, GV, US, h, u, v, tv, Reg, CS, OBC, dt, frac_shelf_h)
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: eta_preale
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h_new ! New 3D grid obtained after last time step [H ~> m or kg m-2]
logical :: PCM_cell(SZI_(G),SZJ_(G),SZK_(GV)) !< If true, PCM remapping should be used in a cell.
integer :: nk, i, j, k, isc, iec, jsc, jec, ntr
integer :: ntr, i, j, k, isc, iec, jsc, jec, nk

nk = GV%ke; isc = G%isc; iec = G%iec; jsc = G%jsc; jec = G%jec
isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec ; nk = GV%ke

if (CS%show_call_tree) call callTree_enter("ALE_main(), MOM_ALE.F90")

Expand All @@ -401,10 +415,19 @@ subroutine ALE_main( G, GV, US, h, u, v, tv, Reg, CS, OBC, dt, frac_shelf_h)
endif
dzRegrid(:,:,:) = 0.0

! If necessary, do some preparatory work to clean up the model state before regridding.

! This adjusts the input thicknesses prior to remapping, based on the verical coordinate.
if (CS%do_conv_adj) call convective_adjustment(G, GV, h, tv)
if (CS%use_hybgen_unmix) then
ntr = 0 ; if (associated(Reg)) ntr = Reg%ntr
call hybgen_unmix(G, GV, G%US, CS%hybgen_unmixCS, tv, Reg, ntr, h)
endif

! Build new grid. The new grid is stored in h_new. The old grid is h.
! Both are needed for the subsequent remapping of variables.
call regridding_main( CS%remapCS, CS%regridCS, G, GV, h, tv, h_new, dzRegrid, &
frac_shelf_h, PCM_cell=PCM_cell)
frac_shelf_h, conv_adjust=.false., PCM_cell=PCM_cell)

call check_grid( G, GV, h, 0. )

Expand Down Expand Up @@ -459,9 +482,9 @@ subroutine ALE_main_offline( G, GV, h, tv, Reg, CS, OBC, dt)
! Local variables
real, dimension(SZI_(G), SZJ_(G), SZK_(GV)+1) :: dzRegrid ! The change in grid interface positions
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h_new ! New 3D grid obtained after last time step [H ~> m or kg m-2]
integer :: nk, i, j, k, isc, iec, jsc, jec
integer :: ntr, i, j, k, isc, iec, jsc, jec, nk

nk = GV%ke; isc = G%isc; iec = G%iec; jsc = G%jsc; jec = G%jec
isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec ; nk = GV%ke

if (CS%show_call_tree) call callTree_enter("ALE_main_offline(), MOM_ALE.F90")

Expand All @@ -470,9 +493,16 @@ subroutine ALE_main_offline( G, GV, h, tv, Reg, CS, OBC, dt)
endif
dzRegrid(:,:,:) = 0.0

! This adjusts the input state prior to remapping, depending on the verical coordinate.
if (CS%do_conv_adj) call convective_adjustment(G, GV, h, tv)
if (CS%use_hybgen_unmix) then
ntr = 0 ; if (associated(Reg)) ntr = Reg%ntr
call hybgen_unmix(G, GV, G%US, CS%hybgen_unmixCS, tv, Reg, ntr, h)
endif

! Build new grid. The new grid is stored in h_new. The old grid is h.
! Both are needed for the subsequent remapping of variables.
call regridding_main( CS%remapCS, CS%regridCS, G, GV, h, tv, h_new, dzRegrid )
call regridding_main( CS%remapCS, CS%regridCS, G, GV, h, tv, h_new, dzRegrid, conv_adjust=.false. )

call check_grid( G, GV, h, 0. )

Expand Down Expand Up @@ -519,7 +549,7 @@ subroutine ALE_offline_inputs(CS, G, GV, h, tv, Reg, uhtr, vhtr, Kd, debug, OBC)
real, dimension(SZK_(GV)) :: h_dest ! Destination grid thicknesses at velocity points [H ~> m or kg m-2]
real, dimension(SZK_(GV)) :: temp_vec ! Transports on the destination grid [H L2 ~> m3 or kg]

nk = GV%ke; isc = G%isc; iec = G%iec; jsc = G%jsc; jec = G%jec
isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec ; nk = GV%ke
dzRegrid(:,:,:) = 0.0
h_new(:,:,:) = 0.0

Expand Down Expand Up @@ -595,13 +625,18 @@ subroutine ALE_offline_tracer_final( G, GV, h, tv, h_target, Reg, CS, OBC)

real, dimension(SZI_(G), SZJ_(G), SZK_(GV)+1) :: dzRegrid !< The change in grid interface positions
real, dimension(SZI_(G), SZJ_(G), SZK_(GV)) :: h_new !< Regridded target thicknesses
integer :: nk, i, j, k, isc, iec, jsc, jec
integer :: ntr, i, j, k, isc, iec, jsc, jec, nk

nk = GV%ke; isc = G%isc; iec = G%iec; jsc = G%jsc; jec = G%jec
isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec ; nk = GV%ke

if (CS%show_call_tree) call callTree_enter("ALE_offline_tracer_final(), MOM_ALE.F90")
! Need to make sure that h_target is consistent with the current offline ALE confiuration
call regridding_main( CS%remapCS, CS%regridCS, G, GV, h_target, tv, h_new, dzRegrid )
if (CS%do_conv_adj) call convective_adjustment(G, GV, h_target, tv)
if (CS%use_hybgen_unmix) then
ntr = 0 ; if (associated(Reg)) ntr = Reg%ntr
call hybgen_unmix(G, GV, G%US, CS%hybgen_unmixCS, tv, Reg, ntr, h)
endif
call regridding_main( CS%remapCS, CS%regridCS, G, GV, h_target, tv, h_new, dzRegrid, conv_adjust=.false. )
call check_grid( G, GV, h_target, 0. )


Expand Down Expand Up @@ -651,52 +686,17 @@ subroutine check_grid( G, GV, h, threshold )

end subroutine check_grid

!### This routine does not appear to be used.
!> Generates new grid
subroutine ALE_build_grid( G, GV, regridCS, remapCS, h, tv, debug, frac_shelf_h )
type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
type(regridding_CS), intent(in) :: regridCS !< Regridding parameters and options
type(remapping_CS), intent(in) :: remapCS !< Remapping parameters and options
type(thermo_var_ptrs), intent(inout) :: tv !< Thermodynamical variable structure
real, dimension(SZI_(G),SZJ_(G), SZK_(GV)), intent(inout) :: h !< Current 3D grid obtained after the
!! last time step [H ~> m or kg m-2]
logical, optional, intent(in) :: debug !< If true, show the call tree
real, dimension(SZI_(G),SZJ_(G)), optional, intent(in):: frac_shelf_h !< Fractional ice shelf coverage [nondim]
! Local variables
integer :: nk, i, j, k
real, dimension(SZI_(G), SZJ_(G), SZK_(GV)+1) :: dzRegrid ! The change in grid interface positions
real, dimension(SZI_(G), SZJ_(G), SZK_(GV)) :: h_new ! The new grid thicknesses
logical :: show_call_tree

show_call_tree = .false.
if (present(debug)) show_call_tree = debug
if (show_call_tree) call callTree_enter("ALE_build_grid(), MOM_ALE.F90")

! Build new grid. The new grid is stored in h_new. The old grid is h.
! Both are needed for the subsequent remapping of variables.
call regridding_main( remapCS, regridCS, G, GV, h, tv, h_new, dzRegrid, frac_shelf_h )

! Override old grid with new one. The new grid 'h_new' is built in
! one of the 'build_...' routines above.
!$OMP parallel do default(none) shared(G,h,h_new)
do j = G%jsc,G%jec ; do i = G%isc,G%iec
if (G%mask2dT(i,j)>0.) h(i,j,:) = h_new(i,j,:)
enddo ; enddo

if (show_call_tree) call callTree_leave("ALE_build_grid()")
end subroutine ALE_build_grid

!> For a state-based coordinate, accelerate the process of regridding by
!! repeatedly applying the grid calculation algorithm
subroutine ALE_regrid_accelerated(CS, G, GV, h, tv, n, u, v, OBC, Reg, dt, dzRegrid, initial)
subroutine ALE_regrid_accelerated(CS, G, GV, h, tv, n_itt, u, v, OBC, Reg, dt, dzRegrid, initial)
type(ALE_CS), pointer :: CS !< ALE control structure
type(ocean_grid_type), intent(inout) :: G !< Ocean grid
type(verticalGrid_type), intent(in) :: GV !< Vertical grid
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
intent(inout) :: h !< Original thicknesses [H ~> m or kg m-2]
type(thermo_var_ptrs), intent(inout) :: tv !< Thermo vars (T/S/EOS)
integer, intent(in) :: n !< Number of times to regrid
integer, intent(in) :: n_itt !< Number of times to regrid
real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
intent(inout) :: u !< Zonal velocity [L T-1 ~> m s-1]
real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
Expand All @@ -711,7 +711,7 @@ subroutine ALE_regrid_accelerated(CS, G, GV, h, tv, n, u, v, OBC, Reg, dt, dzReg
!! routine (and expect diagnostics to work)

! Local variables
integer :: i, j, k, nz
integer :: i, j, itt, nz
type(thermo_var_ptrs) :: tv_local ! local/intermediate temp/salt
type(group_pass_type) :: pass_T_S_h ! group pass if the coordinate has a stencil
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h_loc, h_orig ! A working copy of layer thicknesses
Expand Down Expand Up @@ -744,11 +744,12 @@ subroutine ALE_regrid_accelerated(CS, G, GV, h, tv, n, u, v, OBC, Reg, dt, dzReg
if (present(dt)) &
call ALE_update_regrid_weights(dt, CS)

do k = 1, n
do itt = 1, n_itt
call do_group_pass(pass_T_S_h, G%domain)

! generate new grid
call regridding_main(CS%remapCS, CS%regridCS, G, GV, h_loc, tv_local, h, dzInterface)
if (CS%do_conv_adj) call convective_adjustment(G, GV, h_loc, tv_local)
call regridding_main(CS%remapCS, CS%regridCS, G, GV, h_loc, tv_local, h, dzInterface, conv_adjust=.false.)
dzIntTotal(:,:,:) = dzIntTotal(:,:,:) + dzInterface(:,:,:)

! remap from original grid onto new grid
Expand Down
Loading

0 comments on commit fd2060c

Please sign in to comment.