diff --git a/.testing/tc2.b/MOM_input b/.testing/tc2.b/MOM_input new file mode 120000 index 0000000000..b0cf8cd51c --- /dev/null +++ b/.testing/tc2.b/MOM_input @@ -0,0 +1 @@ +../tc2/MOM_input \ No newline at end of file diff --git a/.testing/tc2.b/MOM_override b/.testing/tc2.b/MOM_override new file mode 100644 index 0000000000..e69de29bb2 diff --git a/.testing/tc2.b/MOM_tc_variant b/.testing/tc2.b/MOM_tc_variant new file mode 100644 index 0000000000..9937095b42 --- /dev/null +++ b/.testing/tc2.b/MOM_tc_variant @@ -0,0 +1,5 @@ +#override TOPO_CONFIG = "spoon" +#override REMAPPING_SCHEME = "PPM_H4" +#override REGRIDDING_COORDINATE_MODE = "ADAPTIVE" + +ADAPT_ADJUSTMENT_SCALE = 0.0 diff --git a/.testing/tc2.b/diag_table b/.testing/tc2.b/diag_table new file mode 120000 index 0000000000..fcf2284f5f --- /dev/null +++ b/.testing/tc2.b/diag_table @@ -0,0 +1 @@ +../tc2/diag_table \ No newline at end of file diff --git a/.testing/tc2.b/input.nml b/.testing/tc2.b/input.nml new file mode 100644 index 0000000000..3c7dcf7bea --- /dev/null +++ b/.testing/tc2.b/input.nml @@ -0,0 +1,20 @@ +&mom_input_nml + output_directory = './' + input_filename = 'n' + restart_input_dir = 'INPUT/' + restart_output_dir = 'RESTART/' + parameter_filename = + 'MOM_input', + 'MOM_tc_variant', + 'MOM_override', +/ + +&diag_manager_nml +/ + +&fms_nml + clock_grain = 'ROUTINE' + clock_flags = 'SYNC' + domains_stack_size = 955296 + stack_size = 0 +/ diff --git a/src/ALE/MOM_ALE.F90 b/src/ALE/MOM_ALE.F90 index e319b71ddc..7c8bbc20a1 100644 --- a/src/ALE/MOM_ALE.F90 +++ b/src/ALE/MOM_ALE.F90 @@ -37,7 +37,8 @@ module MOM_ALE use MOM_regridding, only : regridding_CS, set_regrid_params, write_regrid_file use MOM_regridding, only : getCoordinateInterfaces use MOM_regridding, only : getCoordinateUnits, getCoordinateShortName -use MOM_regridding, only : getStaticThickness +use MOM_regridding, only : getCoordinateMode, getStaticThickness +use MOM_regridding, only : get_adapt_CS use MOM_remapping, only : initialize_remapping, end_remapping use MOM_remapping, only : remapping_core_h, remapping_core_w use MOM_remapping, only : remappingSchemesDoc, remappingDefaultScheme @@ -51,11 +52,14 @@ module MOM_ALE !use regrid_consts, only : coordinateMode, DEFAULT_COORDINATE_MODE use regrid_consts, only : coordinateUnits, coordinateMode, state_dependent +use regrid_consts, only : REGRIDDING_ADAPTIVE use regrid_edge_values, only : edge_values_implicit_h4 use PLM_functions, only : PLM_reconstruction, PLM_boundary_extrapolation use PLM_functions, only : PLM_extrapolate_slope, PLM_monotonized_slope, PLM_slope_wa use PPM_functions, only : PPM_reconstruction, PPM_boundary_extrapolation +use coord_adapt, only : adapt_CS, adapt_diag_CS, associate_adapt_diag, get_adapt_diag_CS + implicit none ; private #include @@ -155,6 +159,7 @@ module MOM_ALE public ALE_remap_init_conds public ALE_register_diags public ALE_set_extrap_boundaries +public ALE_register_coord_diags ! 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 @@ -427,6 +432,140 @@ subroutine ALE_register_diags(Time, G, GV, US, diag, CS) end subroutine ALE_register_diags + +!> Register any diagnostics from within a specific regridding routine +subroutine ALE_register_coord_diags(Time, G, GV, US, diag, CS) + type(time_type), target, intent(in) :: Time !< Time structure + type(ocean_grid_type), intent(in) :: G !< Grid structure + type(unit_scale_type), intent(in) :: US !< Dimensional unit scaling type + type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure + type(diag_ctrl), target, intent(in) :: diag !< Diagnostics control structure + type(ALE_CS), pointer :: CS !< Module control structure + + type(adapt_CS), pointer :: adapt_CS + type(adapt_diag_CS), pointer :: diag_CS + + integer :: isd, ied, jsd, jed, isdB, iedB, jsdB, jedB, nk + + isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed + isdB = G%isdB ; iedB = G%iedB ; jsdB = G%jsdB ; jedB = G%jedB + nk = GV%ke + + select case (getCoordinateMode(CS%regridCS)) + case ( REGRIDDING_ADAPTIVE ) + allocate(diag_CS) + + diag_CS%id_slope_u = register_diag_field("ocean_model", "adapt_slope_u", diag%axesCui, Time, & + "Adaptive coordinate slope on u-points in density space (used for flux calc.)") + diag_CS%id_slope_v = register_diag_field("ocean_model", "adapt_slope_v", diag%axesCvi, Time, & + "Adaptive coordinate slope on v-points in density space (used for flux calc.)") + + diag_CS%id_denom_u = register_diag_field("ocean_model", "adapt_denom_u", diag%axesCui, Time, & + "Adaptive coordinate density denominator on u-points (used for density displacement)") + diag_CS%id_denom_v = register_diag_field("ocean_model", "adapt_denom_v", diag%axesCvi, Time, & + "Adaptive coordinate density denominator on v-points (used for density displacement)") + + diag_CS%id_phys_u = register_diag_field("ocean_model", "adapt_phys_u", diag%axesCui, Time, & + "Adaptive coordinate physical slope on u-points (used as weighting)") + diag_CS%id_phys_v = register_diag_field("ocean_model", "adapt_phys_v", diag%axesCvi, Time, & + "Adaptive coordinate physical slope on v-points (used as weighting)") + + + diag_CS%id_coord_u = register_diag_field("ocean_model", "adapt_coord_u", diag%axesCui, Time, & + "Adaptive coordinate along-coordinate slope on u-points (used as weighting)") + diag_CS%id_coord_v = register_diag_field("ocean_model", "adapt_coord_v", diag%axesCvi, Time, & + "Adaptive coordinate along-coordinate slope on v-points (used as weighting)") + + diag_CS%id_limiting_density = register_diag_field("ocean_model", "adapt_limiting_density", & + diag%axesTi, Time, & + "Adaptive coordinate layer-limiting on density term (difference between "// & + "unlimited and limited flux, before weighting", & + conversion=(GV%H_to_Z * US%Z_to_m) ** 2) + diag_CS%id_limiting_smoothing = register_diag_field("ocean_model", "adapt_limiting_smoothing", & + diag%axesTi, Time, & + "Adaptive coordinate layer-limiting on smoothing term (difference between "// & + "unlimited and limited flux, before weighting", & + conversion=(GV%H_to_Z * US%Z_to_m) ** 2) + + diag_CS%id_w_adjust = register_diag_field("ocean_model", "adapt_w_adjust", diag%axesTi, Time, & + "Adaptive coordinate interface velocity due to hydrostatic adjustment") + + diag_CS%id_disp_density = register_diag_field("ocean_model", "adapt_disp_density", & + diag%axesTi, Time, & + "Adaptive coordinate interface displacement due to density adaptivity", & + conversion=GV%H_to_Z * US%Z_to_m) + diag_CS%id_disp_smoothing = register_diag_field("ocean_model", "adapt_disp_smoothing", & + diag%axesTi, Time, & + "Adaptive coordinate interface displacement due to (limited) smoothing", & + conversion=GV%H_to_Z * US%Z_to_m) + diag_CS%id_disp_unlimited = register_diag_field("ocean_model", "adapt_disp_unlimited", & + diag%axesTi, Time, & + "Adaptive coordinate interface displacement due to (barotropic) smoothing", & + conversion=US%z_to_m) + + if (diag_CS%id_slope_u > 0) allocate(diag_CS%slope_u(isdB:iedB,jsd:jed,nk+1)) + if (diag_CS%id_slope_v > 0) allocate(diag_CS%slope_v(isd:ied,jsdB:jedB,nk+1)) + if (diag_CS%id_denom_u > 0) allocate(diag_CS%denom_u(isdB:iedB,jsd:jed,nk+1)) + if (diag_CS%id_denom_v > 0) allocate(diag_CS%denom_v(isd:ied,jsdB:jedB,nk+1)) + + if (diag_CS%id_phys_u > 0) allocate(diag_CS%phys_u(isdB:iedB,jsd:jed,nk+1)) + if (diag_CS%id_phys_v > 0) allocate(diag_CS%phys_v(isd:ied,jsdB:jedB,nk+1)) + if (diag_CS%id_coord_u > 0) allocate(diag_CS%coord_u(isdB:iedB,jsd:jed,nk+1)) + if (diag_CS%id_coord_v > 0) allocate(diag_CS%coord_v(isd:ied,jsdB:jedB,nk+1)) + + if (diag_CS%id_limiting_density > 0) & + allocate(diag_CS%limiting_density(isd:ied,jsd:jed,nk+1)) + if (diag_CS%id_limiting_smoothing > 0) & + allocate(diag_CS%limiting_smoothing(isd:ied,jsd:jed,nk+1)) + + if (diag_CS%id_w_adjust > 0) & + allocate(diag_CS%w_adjust(isd:ied,jsd:jed,nk+1)) + + if (diag_CS%id_disp_density > 0) & + allocate(diag_CS%disp_density(isd:ied,jsd:jed,nk+1)) + if (diag_CS%id_disp_smoothing > 0) & + allocate(diag_CS%disp_smoothing(isd:ied,jsd:jed,nk+1)) + if (diag_CS%id_disp_unlimited > 0) & + allocate(diag_CS%disp_unlimited(isd:ied,jsd:jed,nk+1)) + + adapt_CS => get_adapt_CS(CS%regridCS) + call associate_adapt_diag(adapt_CS, diag_CS) + end select + +end subroutine ALE_register_coord_diags + +subroutine ALE_post_coord_diags(CS) + type(ALE_CS), pointer :: CS !< Module control structure + + type(adapt_CS), pointer :: adapt_CS + type(adapt_diag_CS), pointer :: diag_CS + + select case (getCoordinateMode(CS%regridCS)) + case ( REGRIDDING_ADAPTIVE ) + diag_CS => get_adapt_diag_CS(get_adapt_CS(CS%regridCS)) + if (.not. associated(diag_CS)) return + + if (diag_CS%id_slope_u > 0) call post_data(diag_CS%id_slope_u, diag_CS%slope_u, CS%diag) + if (diag_CS%id_slope_v > 0) call post_data(diag_CS%id_slope_v, diag_CS%slope_v, CS%diag) + if (diag_CS%id_denom_u > 0) call post_data(diag_CS%id_denom_u, diag_CS%denom_u, CS%diag) + if (diag_CS%id_denom_v > 0) call post_data(diag_CS%id_denom_v, diag_CS%denom_v, CS%diag) + if (diag_CS%id_phys_u > 0) call post_data(diag_CS%id_phys_u, diag_CS%phys_u, CS%diag) + if (diag_CS%id_phys_v > 0) call post_data(diag_CS%id_phys_v, diag_CS%phys_v, CS%diag) + if (diag_CS%id_coord_u > 0) call post_data(diag_CS%id_coord_u, diag_CS%coord_u, CS%diag) + if (diag_CS%id_coord_v > 0) call post_data(diag_CS%id_coord_v, diag_CS%coord_v, CS%diag) + if (diag_CS%id_limiting_density > 0) & + call post_data(diag_CS%id_limiting_density, diag_CS%limiting_density, CS%diag) + if (diag_CS%id_limiting_smoothing > 0) & + call post_data(diag_CS%id_limiting_smoothing, diag_CS%limiting_smoothing, CS%diag) + if (diag_CS%id_w_adjust > 0) call post_data(diag_CS%id_w_adjust, diag_CS%w_adjust, CS%diag) + if (diag_CS%id_disp_density > 0) call post_data(diag_CS%id_disp_density, diag_CS%disp_density, CS%diag) + if (diag_CS%id_disp_smoothing > 0) & + call post_data(diag_CS%id_disp_smoothing, diag_CS%disp_smoothing, CS%diag) + if (diag_CS%id_disp_unlimited > 0) & + call post_data(diag_CS%id_disp_unlimited, diag_CS%disp_unlimited, CS%diag) + end select +end subroutine ALE_post_coord_diags + !> Crudely adjust (initial) grid for integrity. !! This routine is typically called (from initialize_MOM in file MOM.F90) !! before the main time integration loop to initialize the regridding stuff. @@ -519,7 +658,7 @@ end subroutine pre_ALE_adjustments !> Takes care of building a new grid. The creation of the new grid can be based on z coordinates, !! target interface densities, sigma coordinates or any arbitrary coordinate system. -subroutine ALE_regrid( G, GV, US, h, h_new, dzRegrid, tv, CS, frac_shelf_h, PCM_cell) +subroutine ALE_regrid( G, GV, US, h, h_new, dzRegrid, tv, CS, frac_shelf_h, PCM_cell, dt) type(ocean_grid_type), intent(in) :: G !< Ocean grid informations type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type @@ -534,7 +673,8 @@ subroutine ALE_regrid( G, GV, US, h, h_new, dzRegrid, tv, CS, frac_shelf_h, PCM_ type(ALE_CS), pointer :: CS !< Regridding parameters and options real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: frac_shelf_h !< Fractional ice shelf coverage [nondim] logical, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & - optional, intent(out) :: PCM_cell !< If true, use PCM remapping in a cell. + optional, intent(out) :: PCM_cell !< If true, use PCM remapping in a cell. + real, optional, intent(in) :: dt !< Time step between calls to ALE_regrid [T ~> s] ! Local variables logical :: showCallTree @@ -547,7 +687,7 @@ subroutine ALE_regrid( G, GV, US, h, h_new, dzRegrid, tv, CS, frac_shelf_h, PCM_ ! Both are needed for the subsequent remapping of variables. dzRegrid(:,:,:) = 0.0 call regridding_main( CS%remapCS, CS%regridCS, G, GV, US, h, tv, h_new, dzRegrid, & - frac_shelf_h=frac_shelf_h, PCM_cell=PCM_cell) + frac_shelf_h=frac_shelf_h, PCM_cell=PCM_cell, dt=dt) if (CS%id_dzRegrid>0) then ; if (query_averaging_enabled(CS%diag)) then call post_data(CS%id_dzRegrid, dzRegrid, CS%diag, alt_h=h_new) @@ -725,7 +865,7 @@ subroutine ALE_regrid_accelerated(CS, G, GV, US, h, tv, n_itt, u, v, OBC, Reg, d ! Update the layer specific volumes if necessary if (allocated(tv_local%SpV_avg)) call calc_derived_thermo(tv_local, h, G, GV, US, halo=1) - call regridding_main(CS%remapCS, CS%regridCS, G, GV, US, h_loc, tv_local, h, dzInterface) + call regridding_main(CS%remapCS, CS%regridCS, G, GV, US, h_loc, tv_local, h, dzInterface, dt=dt) if (CS%remap_uv_using_old_alg) & dzIntTotal(:,:,:) = dzIntTotal(:,:,:) + dzInterface(:,:,:) diff --git a/src/ALE/MOM_regridding.F90 b/src/ALE/MOM_regridding.F90 index 7e24d80a21..fe489d9eb2 100644 --- a/src/ALE/MOM_regridding.F90 +++ b/src/ALE/MOM_regridding.F90 @@ -21,6 +21,7 @@ module MOM_regridding use MOM_string_functions, only : uppercase, extractWord, extract_integer, extract_real use MOM_remapping, only : remapping_CS +use filter_utils, only : filter_CS, filtered_grid_motion use regrid_consts, only : state_dependent, coordinateUnits use regrid_consts, only : coordinateMode, DEFAULT_COORDINATE_MODE use regrid_consts, only : REGRIDDING_LAYER, REGRIDDING_ZSTAR @@ -40,7 +41,7 @@ module MOM_regridding use coord_hycom, only : init_coord_hycom, set_hycom_params, build_hycom1_column, end_coord_hycom use coord_hycom, only : init_3d_coord_hycom use coord_adapt, only : adapt_CS -use coord_adapt, only : init_coord_adapt, set_adapt_params, build_adapt_column, end_coord_adapt +use coord_adapt, only : init_coord_adapt, set_adapt_params, build_adapt_grid, end_coord_adapt use MOM_hybgen_regrid, only : hybgen_regrid, hybgen_regrid_CS, init_hybgen_regrid, end_hybgen_regrid use MOM_hybgen_regrid, only : write_Hybgen_coord_file @@ -105,17 +106,8 @@ module MOM_regridding !> Reference pressure for potential density calculations [R L2 T-2 ~> Pa] real :: ref_pressure = 2.e7 - !> Weight given to old coordinate when blending between new and old grids [nondim] - !! Used only below depth_of_time_filter_shallow, with a cubic variation - !! from zero to full effect between depth_of_time_filter_shallow and - !! depth_of_time_filter_deep. - real :: old_grid_weight = 0. - - !> Depth above which no time-filtering of grid is applied [H ~> m or kg m-2] - real :: depth_of_time_filter_shallow = 0. - - !> Depth below which time-filtering of grid is applied at full effect [H ~> m or kg m-2] - real :: depth_of_time_filter_deep = 0. + !> The control structure for coordinate filtering + type(filter_CS) :: filter_CS !> Fraction (between 0 and 1) of compressibility to add to potential density !! profiles when interpolating for target grid positions [nondim] @@ -153,11 +145,11 @@ module MOM_regridding public uniformResolution, setCoordinateResolution public set_target_densities_from_GV, set_target_densities public set_regrid_max_depths, set_regrid_max_thickness -public getCoordinateResolution, getCoordinateInterfaces +public getCoordinateResolution, getCoordinateInterfaces, getCoordinateMode public getCoordinateUnits, getCoordinateShortName, getStaticThickness public DEFAULT_COORDINATE_MODE public set_h_neglect, set_dz_neglect -public get_zlike_CS, get_sigma_CS, get_rho_CS +public get_zlike_CS, get_sigma_CS, get_rho_CS, get_adapt_CS !> Documentation for coordinate options character(len=*), parameter, public :: regriddingCoordinateModeDoc = & @@ -245,6 +237,7 @@ subroutine initialize_regridding(CS, G, GV, US, max_depth, param_file, mdl, & real :: adaptZoom ! The thickness of the near-surface zooming region with the adaptive coordinate [H ~> m or kg m-2] real :: adaptDrho0 ! Reference density difference for stratification-dependent diffusion. [R ~> kg m-3] integer :: i, j, k, nzf(4) + real :: adapt_alpha_rho, adapt_alpha_p, adapt_timescale, adapt_cutoff, adapt_smooth, adapt_adjustment real, dimension(:), allocatable :: dz ! Resolution (thickness) in units of coordinate, which may be [m] ! or [Z ~> m] or [H ~> m or kg m-2] or [R ~> kg m-3] or other units. real, dimension(:,:), allocatable :: dz_2d ! 2D resolution (thickness) in units of coordinate, which may be [m] @@ -922,8 +915,8 @@ subroutine initialize_regridding(CS, G, GV, US, max_depth, param_file, mdl, & elseif (coordinateMode(coord_mode) == REGRIDDING_RHO) then call setCoordinateResolution(dz, CS, scale=US%kg_m3_to_R) elseif (coordinateMode(coord_mode) == REGRIDDING_ADAPTIVE) then - call setCoordinateResolution(dz, CS, scale=GV%m_to_H) - CS%coord_scale = GV%H_to_m + call setCoordinateResolution(dz, CS, scale=US%m_to_Z) + CS%coord_scale = US%Z_to_m else call setCoordinateResolution(dz, CS, scale=US%m_to_Z) CS%coord_scale = US%Z_to_m @@ -1002,27 +995,55 @@ subroutine initialize_regridding(CS, G, GV, US, max_depth, param_file, mdl, & endif if (coordinateMode(coord_mode) == REGRIDDING_ADAPTIVE) then - call get_param(param_file, mdl, "ADAPT_TIME_RATIO", adaptTimeRatio, & - "Ratio of ALE timestep to grid timescale.", units="nondim", default=1.0e-1) - call get_param(param_file, mdl, "ADAPT_ZOOM_DEPTH", adaptZoom, & - "Depth of near-surface zooming region.", units="m", default=200.0, scale=GV%m_to_H) - call get_param(param_file, mdl, "ADAPT_ZOOM_COEFF", adaptZoomCoeff, & - "Coefficient of near-surface zooming diffusivity.", units="nondim", default=0.2) - call get_param(param_file, mdl, "ADAPT_BUOY_COEFF", adaptBuoyCoeff, & - "Coefficient of buoyancy diffusivity.", units="nondim", default=0.8) - call get_param(param_file, mdl, "ADAPT_ALPHA", adaptAlpha, & - "Scaling on optimization tendency.", units="nondim", default=1.0) - call get_param(param_file, mdl, "ADAPT_DO_MIN_DEPTH", tmpLogical, & - "If true, make a HyCOM-like mixed layer by preventing interfaces "//& - "from being shallower than the depths specified by the regridding coordinate.", & - default=.false.) - call get_param(param_file, mdl, "ADAPT_DRHO0", adaptDrho0, & - "Reference density difference for stratification-dependent diffusion.", & - units="kg m-3", default=0.5, scale=US%kg_m3_to_R) - - call set_regrid_params(CS, adaptTimeRatio=adaptTimeRatio, adaptZoom=adaptZoom, & - adaptZoomCoeff=adaptZoomCoeff, adaptBuoyCoeff=adaptBuoyCoeff, adaptAlpha=adaptAlpha, & - adaptDoMin=tmpLogical, adaptDrho0=adaptDrho0) + call get_param(param_file, mdl, "ADAPT_ALPHA_RHO", adapt_alpha_rho, & + "Density adaptivity coefficient (use negative value for automatic)", & + units="nondim", default=-1.0) + call get_param(param_file, mdl, "ADAPT_ALPHA_P", adapt_alpha_p, & + "Pressure adaptivity coefficient (use negative value for automatic)", & + units="nondim", default=-1.0) + call get_param(param_file, mdl, "ADAPT_TIMESCALE", adapt_timescale, & + "Timescale for adaptivity diffusivity (defaults to a day)", & + units="s", default=86400.0, scale=US%s_to_T) + call get_param(param_file, mdl, "ADAPT_MEAN_H", tmpLogical, & + "Use mean rather than 'upstream' h in calculations", default=.false.) + call get_param(param_file, mdl, "ADAPT_SLOPE_CUTOFF", adapt_cutoff, & + "Slope cutoff between stratified and unstratified regions", units="nondim", default=1e-2) + call get_param(param_file, mdl, "ADAPT_SMOOTH_MIN", adapt_smooth, & + "Minimum weight toward smoothing term", units="nondim", default=0.) + + call get_param(param_file, mdl, "ADAPT_ADJUSTMENT_SCALE", adapt_adjustment, & + "Non-dimensional scale for adjusting interface positions when\n"//& + "a diagonal convective instability would occur. When set to 1,\n"//& + "perform the full adjustment permitted by the local CFL value.", & + units="nondim", default=1.0) + + call set_regrid_params(CS, adapt_alpha_rho=adapt_alpha_rho, adapt_alpha_p=adapt_alpha_p, & + adapt_timescale=adapt_timescale, adapt_mean=tmpLogical, & + adapt_cutoff=adapt_cutoff, adapt_smooth=adapt_smooth, adapt_adjustment_scale=adapt_adjustment) + + call get_param(param_file, mdl, "ADAPT_RESTORING_TIMESCALE", adapt_timescale, & + "Timescale for adaptivity restoring (default 10 days)", & + units="s", default=864000.0, scale=US%s_to_T) + call set_regrid_params(CS, adapt_restoring_timescale=adapt_timescale) + + call get_param(param_file, mdl, "ADAPT_TWIN_GRADIENT", tmpLogical, & + "Use twin gradient approach, requiring sign of gradient above\n"//& + "and below an interface to agree, avoiding a vertical null mode.", & + default=.true.) + call set_regrid_params(CS, adapt_twin=tmpLogical) + + call get_param(param_file, mdl, "ADAPT_PHYSICAL_SLOPE", tmpLogical, & + "Use a physical slope calculation for weighting of terms.", & + default=.true.) + call set_regrid_params(CS, adapt_physical_slope=tmpLogical) + + call get_param(param_file, mdl, "ADAPT_RESTORE_MEAN", tmpLogical, & + "If True, restore toward a dynamically-calculated mean interface position.\n"//& + "Otherwise, restore toward the profile given by ALE_COORDINATE_CONFIG.", & + default=.false.) + + call set_regrid_params(CS, adapt_restore_mean=tmpLogical, & + adapt_restoring_timescale=adapt_timescale) endif if (main_parameters .and. coord_is_state_dependent) then @@ -1179,7 +1200,7 @@ end subroutine end_regridding !------------------------------------------------------------------------------ !> Dispatching regridding routine for orchestrating regridding & remapping subroutine regridding_main( remapCS, CS, G, GV, US, h, tv, h_new, dzInterface, & - frac_shelf_h, PCM_cell) + frac_shelf_h, PCM_cell, dt) !------------------------------------------------------------------------------ ! This routine takes care of (1) building a new grid and (2) remapping between ! the old grid and the new grid. The creation of the new grid can be based @@ -1213,6 +1234,7 @@ subroutine regridding_main( remapCS, CS, G, GV, US, h, tv, h_new, dzInterface, & real, dimension(SZI_(G),SZJ_(G)), optional, intent(in ) :: frac_shelf_h !< Fractional ice shelf coverage [nomdim] logical, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & optional, intent(out ) :: PCM_cell !< Use PCM remapping in cells where true + real, optional, intent(in) :: dt !< Current model timestep ! Local variables real :: nom_depth_H(SZI_(G),SZJ_(G)) !< The nominal ocean depth at each point in thickness units [H ~> m or kg m-2] @@ -1278,7 +1300,7 @@ subroutine regridding_main( remapCS, CS, G, GV, US, h, tv, h_new, dzInterface, & call hybgen_regrid(G, GV, G%US, h, nom_depth_H, tv, CS%hybgen_CS, dzInterface, PCM_cell) call calc_h_new_by_dz(CS, G, GV, h, dzInterface, h_new) case ( REGRIDDING_ADAPTIVE ) - call build_grid_adaptive(G, GV, G%US, h, nom_depth_H, tv, dzInterface, remapCS, CS) + call build_grid_adaptive(G, GV, G%US, h, tv, CS, dzInterface, dt) call calc_h_new_by_dz(CS, G, GV, h, dzInterface, h_new) case ( REGRIDDING_ARBITRARY ) @@ -1432,159 +1454,6 @@ subroutine check_grid_column( nk, h, dzInterface, msg ) end subroutine check_grid_column -!> Returns the change in interface position motion after filtering and -!! assuming the top and bottom interfaces do not move. The filtering is -!! a function of depth, and is applied as the integrated average filtering -!! over the trajectory of the interface. By design, this code can not give -!! tangled interfaces provided that z_old and z_new are not already tangled. -subroutine filtered_grid_motion( CS, nk, z_old, z_new, dz_g ) - type(regridding_CS), intent(in) :: CS !< Regridding control structure - integer, intent(in) :: nk !< Number of cells in source grid - real, dimension(nk+1), intent(in) :: z_old !< Old grid position [H ~> m or kg m-2] - real, dimension(CS%nk+1), intent(in) :: z_new !< New grid position before filtering [H ~> m or kg m-2] - real, dimension(CS%nk+1), intent(inout) :: dz_g !< Change in interface positions including - !! the effects of filtering [H ~> m or kg m-2] - ! Local variables - real :: sgn ! The sign convention for downward [nondim]. - real :: dz_tgt ! The target grid movement of the unfiltered grid [H ~> m or kg m-2] - real :: zr1 ! The old grid position of an interface relative to the surface [H ~> m or kg m-2] - real :: z_old_k ! The corrected position of the old grid [H ~> m or kg m-2] - real :: Aq ! A temporary variable related to the grid weights [nondim] - real :: Bq ! A temporary variable used in the linear term in the quadratic expression for the - ! filtered grid movement [H ~> m or kg m-2] - real :: z0, dz0 ! Together these give the position of an interface relative to a reference hieght - ! that may be adjusted for numerical accuracy in a solver [H ~> m or kg m-2] - real :: F0 ! An estimated grid movement [H ~> m or kg m-2] - real :: zs ! The depth at which the shallow filtering timescale applies [H ~> m or kg m-2] - real :: zd ! The depth at which the deep filtering timescale applies [H ~> m or kg m-2] - real :: dzwt ! The depth range over which the transition in the filtering timescale occurs [H ~> m or kg m-2] - real :: Idzwt ! The Adcroft reciprocal of dzwt [H-1 ~> m-1 or m2 kg-1] - real :: wtd ! The weight given to the new grid when time filtering [nondim] - real :: Iwtd ! The inverse of wtd [nondim] - real :: Int_zs ! A depth integral of the weights in [H ~> m or kg m-2] - real :: Int_zd ! A depth integral of the weights in [H ~> m or kg m-2] - real :: dInt_zs_zd ! The depth integral of the weights between the deep and shallow depths in [H ~> m or kg m-2] -! For debugging: - real, dimension(nk+1) :: z_act ! The final grid positions after the filtered movement [H ~> m or kg m-2] -! real, dimension(nk+1) :: ddz_g_s, ddz_g_d - logical :: debug = .false. - integer :: k - - if ((z_old(nk+1) - z_old(1)) * (z_new(CS%nk+1) - z_new(1)) < 0.0) then - call MOM_error(FATAL, "filtered_grid_motion: z_old and z_new use different sign conventions.") - elseif ((z_old(nk+1) - z_old(1)) * (z_new(CS%nk+1) - z_new(1)) == 0.0) then - ! This is a massless column, so do nothing and return. - do k=1,CS%nk+1 ; dz_g(k) = 0.0 ; enddo ; return - elseif ((z_old(nk+1) - z_old(1)) + (z_new(CS%nk+1) - z_new(1)) > 0.0) then - sgn = 1.0 - else - sgn = -1.0 - endif - - if (debug) then - do k=2,CS%nk+1 - if (sgn*(z_new(k)-z_new(k-1)) < -5e-16*(abs(z_new(k))+abs(z_new(k-1))) ) & - call MOM_error(FATAL, "filtered_grid_motion: z_new is tangled.") - enddo - do k=2,nk+1 - if (sgn*(z_old(k)-z_old(k-1)) < -5e-16*(abs(z_old(k))+abs(z_old(k-1))) ) & - call MOM_error(FATAL, "filtered_grid_motion: z_old is tangled.") - enddo - ! ddz_g_s(:) = 0.0 ; ddz_g_d(:) = 0.0 - endif - - zs = CS%depth_of_time_filter_shallow - zd = CS%depth_of_time_filter_deep - wtd = 1.0 - CS%old_grid_weight - Iwtd = 1.0 / wtd - - dzwt = (zd - zs) - Idzwt = 0.0 ; if (abs(zd - zs) > 0.0) Idzwt = 1.0 / (zd - zs) - dInt_zs_zd = 0.5*(1.0 + Iwtd) * (zd - zs) - Aq = 0.5*(Iwtd - 1.0) - - dz_g(1) = 0.0 - z_old_k = z_old(1) - do k = 2,CS%nk+1 - if (k<=nk+1) z_old_k = z_old(k) ! This allows for virtual z_old interface at bottom of the model - ! zr1 is positive and increases with depth, and dz_tgt is positive downward. - dz_tgt = sgn*(z_new(k) - z_old_k) - zr1 = sgn*(z_old_k - z_old(1)) - - ! First, handle the two simple and common cases that do not pass through - ! the adjustment rate transition zone. - if ((zr1 > zd) .and. (zr1 + wtd * dz_tgt > zd)) then - dz_g(k) = sgn * wtd * dz_tgt - elseif ((zr1 < zs) .and. (zr1 + dz_tgt < zs)) then - dz_g(k) = sgn * dz_tgt - else - ! Find the new value by inverting the equation - ! integral(0 to dz_new) Iwt(z) dz = dz_tgt - ! This is trivial where Iwt is a constant, and agrees with the two limits above. - - ! Take test values at the transition points to figure out which segment - ! the new value will be found in. - if (zr1 >= zd) then - Int_zd = Iwtd*(zd - zr1) - Int_zs = Int_zd - dInt_zs_zd - elseif (zr1 <= zs) then - Int_zs = (zs - zr1) - Int_zd = dInt_zs_zd + (zs - zr1) - else -! Int_zd = (zd - zr1) * (Iwtd + 0.5*(1.0 - Iwtd) * (zd - zr1) / (zd - zs)) - Int_zd = (zd - zr1) * (Iwtd*(0.5*(zd+zr1) - zs) + 0.5*(zd - zr1)) * Idzwt - Int_zs = (zs - zr1) * (0.5*Iwtd * ((zr1 - zs)) + (zd - 0.5*(zr1+zs))) * Idzwt - ! It has been verified that Int_zs = Int_zd - dInt_zs_zd to within roundoff. - endif - - if (dz_tgt >= Int_zd) then ! The new location is in the deep, slow region. - dz_g(k) = sgn * ((zd-zr1) + wtd*(dz_tgt - Int_zd)) - elseif (dz_tgt <= Int_zs) then ! The new location is in the shallow region. - dz_g(k) = sgn * ((zs-zr1) + (dz_tgt - Int_zs)) - else ! We need to solve a quadratic equation for z_new. - ! For accuracy, do the integral from the starting depth or the nearest - ! edge of the transition region. The results with each choice are - ! mathematically equivalent, but differ in roundoff, and this choice - ! should minimize the likelihood of inadvertently overlapping interfaces. - if (zr1 <= zs) then ; dz0 = zs-zr1 ; z0 = zs ; F0 = dz_tgt - Int_zs - elseif (zr1 >= zd) then ; dz0 = zd-zr1 ; z0 = zd ; F0 = dz_tgt - Int_zd - else ; dz0 = 0.0 ; z0 = zr1 ; F0 = dz_tgt ; endif - - Bq = (dzwt + 2.0*Aq*(z0-zs)) - ! Solve the quadratic: Aq*(zn-z0)**2 + Bq*(zn-z0) - F0*dzwt = 0 - ! Note that b>=0, and the two terms in the standard form cancel for the right root. - dz_g(k) = sgn * (dz0 + 2.0*F0*dzwt / (Bq + sqrt(Bq**2 + 4.0*Aq*F0*dzwt) )) - -! if (debug) then -! dz0 = zs-zr1 ; z0 = zs ; F0 = dz_tgt - Int_zs ; Bq = (dzwt + 2.0*Aq*(z0-zs)) -! ddz_g_s(k) = sgn * (dz0 + 2.0*F0*dzwt / (Bq + sqrt(Bq**2 + 4.0*Aq*F0*dzwt) )) - dz_g(k) -! dz0 = zd-zr1 ; z0 = zd ; F0 = dz_tgt - Int_zd ; Bq = (dzwt + 2.0*Aq*(z0-zs)) -! ddz_g_d(k) = sgn * (dz0 + 2.0*F0*dzwt / (Bq + sqrt(Bq**2 + 4.0*Aq*F0*dzwt) )) - dz_g(k) -! -! if (abs(ddz_g_s(k)) > 1e-12*(abs(dz_g(k)) + abs(dz_g(k)+ddz_g_s(k)))) & -! call MOM_error(WARNING, "filtered_grid_motion: Expect z_output to be tangled (sc).") -! if (abs(ddz_g_d(k) - ddz_g_s(k)) > 1e-12*(abs(dz_g(k)+ddz_g_d(k)) + abs(dz_g(k)+ddz_g_s(k)))) & -! call MOM_error(WARNING, "filtered_grid_motion: Expect z_output to be tangled.") -! endif - endif - - endif - enddo - !dz_g(CS%nk+1) = 0.0 - - if (debug) then - z_old_k = z_old(1) - do k=1,CS%nk+1 - if (k<=nk+1) z_old_k = z_old(k) ! This allows for virtual z_old interface at bottom of the model - z_act(k) = z_old_k + dz_g(k) - enddo - do k=2,CS%nk+1 - if (sgn*((z_act(k))-z_act(k-1)) < -1e-15*(abs(z_act(k))+abs(z_act(k-1))) ) & - call MOM_error(FATAL, "filtered_grid_motion: z_output is tangled.") - enddo - endif - -end subroutine filtered_grid_motion !> Builds a z*-coordinate grid with partial steps (Adcroft and Campin, 2004). !! z* is defined as @@ -1643,7 +1512,7 @@ subroutine build_zstar_grid( CS, G, GV, h, nom_depth_H, dzInterface, frac_shelf_ ! Determine water column thickness totalThickness = 0.0 do k = 1,nz - totalThickness = totalThickness + h(i,j,k) + totalThickness = totalThickness + h(i,j,k) enddo ! if (GV%Boussinesq) then @@ -1669,7 +1538,7 @@ subroutine build_zstar_grid( CS, G, GV, h, nom_depth_H, dzInterface, frac_shelf_ endif ! Calculate the final change in grid position after blending new and old grids - call filtered_grid_motion( CS, nz, zOld, zNew, dzInterface(i,j,:) ) + call filtered_grid_motion( CS%filter_CS, CS%nk, nz, zOld, zNew, dzInterface(i,j,:) ) #ifdef __DO_SAFETY_CHECKS__ dh = max(nominalDepth,totalThickness) @@ -1768,7 +1637,7 @@ subroutine build_sigma_grid( CS, G, GV, h, nom_depth_H, dzInterface ) zOld(k) = zOld(k+1) + h(i,j,k) enddo - call filtered_grid_motion( CS, nz, zOld, zNew, dzInterface(i,j,:) ) + call filtered_grid_motion( CS%filter_CS, CS%nk, nz, zOld, zNew, dzInterface(i,j,:) ) #ifdef __DO_SAFETY_CHECKS__ dh = max(nominalDepth,totalThickness) @@ -1915,7 +1784,7 @@ subroutine build_rho_grid( G, GV, US, h, nom_depth_H, tv, dzInterface, remapCS, endif ! Calculate the final change in grid position after blending new and old grids - call filtered_grid_motion( CS, nz, zOld, zNew, dzInterface(i,j,:) ) + call filtered_grid_motion( CS%filter_CS, CS%nk, nz, zOld, zNew, dzInterface(i,j,:) ) #ifdef __DO_SAFETY_CHECKS__ do k=2,CS%nk @@ -2045,7 +1914,7 @@ subroutine build_grid_HyCOM1( G, GV, US, h, nom_depth_H, tv, h_new, dzInterface, h_neglect=h_neglect, h_neglect_edge=h_neglect_edge) ! Calculate the final change in grid position after blending new and old grids - call filtered_grid_motion( CS, GV%ke, z_col, z_col_new, dz_col ) + call filtered_grid_motion( CS%filter_CS, CS%nk, GV%ke, z_col, z_col_new, dz_col ) ! This adjusts things robust to round-off errors dz_col(:) = -dz_col(:) @@ -2063,72 +1932,33 @@ subroutine build_grid_HyCOM1( G, GV, US, h, nom_depth_H, tv, h_new, dzInterface, end subroutine build_grid_HyCOM1 -!> This subroutine builds an adaptive grid that follows density surfaces where -!! possible, subject to constraints on the smoothness of interface heights. -subroutine build_grid_adaptive(G, GV, US, h, nom_depth_H, tv, dzInterface, remapCS, CS) +!> Build a grid using the AG adaptive-density/smoothing algorithm +subroutine build_grid_adaptive(G, GV, US, h, tv, CS, dzInterface, dt) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] - real, dimension(SZI_(G),SZJ_(G)), intent(in) :: nom_depth_H !< The bathymetric depth of this column - !! relative to mean sea level or another locally - !! valid reference height, converted to thickness - !! units [H ~> m or kg m-2] type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various !! thermodynamic variables - type(regridding_CS), intent(in) :: CS !< Regridding control structure + type(regridding_CS), intent(in) :: CS !< Regridding control structure real, dimension(SZI_(G),SZJ_(G),CS%nk+1), intent(inout) :: dzInterface !< The change in interface depth !! [H ~> m or kg m-2] - type(remapping_CS), intent(in) :: remapCS !< The remapping control structure + real, optional, intent(in) :: dt !< The intended timestep over which this + !! regridding operation applies - ! local variables - integer :: i, j, k, nz ! indices and dimension lengths - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: tInt ! Temperature on interfaces [C ~> degC] - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: sInt ! Salinity on interfaces [S ~> ppt] - ! current interface positions and after tendency term is applied - ! positive downward - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: zInt ! Interface depths [H ~> m or kg m-2] - real, dimension(SZK_(GV)+1) :: zNext ! New interface depths [H ~> m or kg m-2] + integer :: i, j - nz = GV%ke + call build_adapt_grid(G, GV, US, h, tv, dzInterface, CS%adapt_CS, CS%filter_CS, CS%min_thickness, dt) call assert((GV%ke == CS%nk), "build_grid_adaptive is only written to work "//& "with the same number of input and target layers.") - ! position surface at z = 0. - zInt(:,:,1) = 0. - - ! work on interior interfaces - do K = 2, nz ; do j = G%jsc-2,G%jec+2 ; do i = G%isc-2,G%iec+2 - tInt(i,j,K) = 0.5 * (tv%T(i,j,k-1) + tv%T(i,j,k)) - sInt(i,j,K) = 0.5 * (tv%S(i,j,k-1) + tv%S(i,j,k)) - zInt(i,j,K) = zInt(i,j,K-1) + h(i,j,k-1) ! zInt in [H] - enddo ; enddo ; enddo - - ! top and bottom temp/salt interfaces are just the layer - ! average values - tInt(:,:,1) = tv%T(:,:,1) ; tInt(:,:,nz+1) = tv%T(:,:,nz) - sInt(:,:,1) = tv%S(:,:,1) ; sInt(:,:,nz+1) = tv%S(:,:,nz) - - ! set the bottom interface depth - zInt(:,:,nz+1) = zInt(:,:,nz) + h(:,:,nz) - - ! calculate horizontal density derivatives (alpha/beta) - ! between cells in a 5-point stencil, columnwise - do j = G%jsc-1,G%jec+1 ; do i = G%isc-1,G%iec+1 - if (G%mask2dT(i,j) < 0.5) then - dzInterface(i,j,:) = 0. ! land point, don't move interfaces, and skip - cycle - endif - - call build_adapt_column(CS%adapt_CS, G, GV, US, tv, i, j, zInt, tInt, sInt, h, & - nom_depth_H, zNext) + do j = G%jsc-1,G%jec+1 + do i = G%isc-1,G%iec+1 + call check_grid_column(GV%ke, h(i,j,:), dzInterface(i,j,:), 'in build_grid_adaptive') + enddo + enddo - call filtered_grid_motion(CS, nz, zInt(i,j,:), zNext, dzInterface(i,j,:)) - ! convert from depth to z - do K = 1, nz+1 ; dzInterface(i,j,K) = -dzInterface(i,j,K) ; enddo - call adjust_interface_motion(CS, nz, h(i,j,:), dzInterface(i,j,:)) - enddo ; enddo end subroutine build_grid_adaptive !> Adjust dz_Interface to ensure non-negative future thicknesses @@ -2385,7 +2215,7 @@ subroutine initCoord(CS, G, GV, US, coord_mode, param_file) case (REGRIDDING_HYBGEN) call init_hybgen_regrid(CS%hybgen_CS, GV, US, param_file) case (REGRIDDING_ADAPTIVE) - call init_coord_adapt(CS%adapt_CS, CS%nk, CS%coordinateResolution, GV%m_to_H, US%kg_m3_to_R) + call init_coord_adapt(CS%adapt_CS, CS%nk, CS%coordinateResolution) end select end subroutine initCoord @@ -2699,6 +2529,14 @@ function getCoordinateInterfaces( CS, undo_scaling ) end function getCoordinateInterfaces +!> Query the regrdding scheme (coordinate mode) +function getCoordinateMode( CS ) + type(regridding_CS), intent(in) :: CS !< Regridding control structure + integer :: getCoordinateMode + + getCoordinateMode = CS%regridding_scheme +end function getCoordinateMode + !------------------------------------------------------------------------------ !> Query the target coordinate units function getCoordinateUnits( CS ) @@ -2761,8 +2599,9 @@ subroutine set_regrid_params( CS, boundary_extrapolation, min_thickness, old_gri interp_scheme, depth_of_time_filter_shallow, depth_of_time_filter_deep, & compress_fraction, ref_pressure, & integrate_downward_for_e, remap_answers_2018, remap_answer_date, regrid_answer_date, & - adaptTimeRatio, adaptZoom, adaptZoomCoeff, adaptBuoyCoeff, & - adaptAlpha, adaptDoMin, adaptDrho0) + adapt_alpha_rho, adapt_alpha_p, adapt_timescale, & + adapt_mean, adapt_twin, adapt_cutoff, adapt_smooth, adapt_physical_slope, & + adapt_restoring_timescale, adapt_restore_mean, adapt_adjustment_scale) type(regridding_CS), intent(inout) :: CS !< Regridding control structure logical, optional, intent(in) :: boundary_extrapolation !< Extrapolate in boundary cells real, optional, intent(in) :: min_thickness !< Minimum thickness allowed when building the @@ -2781,16 +2620,18 @@ subroutine set_regrid_params( CS, boundary_extrapolation, min_thickness, old_gri !! use more robust but mathematically equivalent expressions. integer, optional, intent(in) :: remap_answer_date !< The vintage of the expressions to use for remapping integer, optional, intent(in) :: regrid_answer_date !< The vintage of the expressions to use for regridding - real, optional, intent(in) :: adaptTimeRatio !< Ratio of the ALE timestep to the grid timescale [nondim]. - real, optional, intent(in) :: adaptZoom !< Depth of near-surface zooming region [H ~> m or kg m-2]. - real, optional, intent(in) :: adaptZoomCoeff !< Coefficient of near-surface zooming diffusivity [nondim]. - real, optional, intent(in) :: adaptBuoyCoeff !< Coefficient of buoyancy diffusivity [nondim]. - real, optional, intent(in) :: adaptAlpha !< Scaling factor on optimization tendency [nondim]. - logical, optional, intent(in) :: adaptDoMin !< If true, make a HyCOM-like mixed layer by - !! preventing interfaces from being shallower than - !! the depths specified by the regridding coordinate. - real, optional, intent(in) :: adaptDrho0 !< Reference density difference for stratification-dependent - !! diffusion. [R ~> kg m-3] + real, optional, intent(in) :: adapt_alpha_rho !< Manually-specified weight for density adaptivity + real, optional, intent(in) :: adapt_alpha_p !< Manually-specified weight for pressure adaptivity + real, optional, intent(in) :: adapt_timescale !< Timescale over which to apply adaptivity terms + real, optional, intent(in) :: adapt_restoring_timescale !< Timescale for coordinate restoration + real, optional, intent(in) :: adapt_cutoff !< Interface slope cutoff defining stratified/unstratified regions + real, optional, intent(in) :: adapt_smooth !< Minimum weight for smoothing term + real, optional, intent(in) :: adapt_adjustment_scale !< Non-dimensional scale for diagonal convective instability + logical, optional, intent(in) :: adapt_mean !< Use mean rather than "upstream" thickness + logical, optional, intent(in) :: adapt_twin !< Calculate sign of density gradient above and below interfaces + logical, optional, intent(in) :: adapt_physical_slope !< Use along-coordinate or physical-space slope? + logical, optional, intent(in) :: adapt_restore_mean !< Restore towards dynamically-calculated interface mean, + !! or specified coordinate if (present(interp_scheme)) call set_interp_scheme(CS%interp_CS, interp_scheme) if (present(boundary_extrapolation)) call set_interp_extrap(CS%interp_CS, boundary_extrapolation) @@ -2799,15 +2640,15 @@ subroutine set_regrid_params( CS, boundary_extrapolation, min_thickness, old_gri if (present(old_grid_weight)) then if (old_grid_weight<0. .or. old_grid_weight>1.) & call MOM_error(FATAL,'MOM_regridding, set_regrid_params: Weight is out side the range 0..1!') - CS%old_grid_weight = old_grid_weight + CS%filter_CS%old_grid_weight = old_grid_weight endif - if (present(depth_of_time_filter_shallow)) CS%depth_of_time_filter_shallow = & + if (present(depth_of_time_filter_shallow)) CS%filter_CS%depth_of_time_filter_shallow = & depth_of_time_filter_shallow - if (present(depth_of_time_filter_deep)) CS%depth_of_time_filter_deep = & + if (present(depth_of_time_filter_deep)) CS%filter_CS%depth_of_time_filter_deep = & depth_of_time_filter_deep if (present(depth_of_time_filter_shallow) .or. present(depth_of_time_filter_deep)) then - if (CS%depth_of_time_filter_deep This returns a pointer to the adapt_CS stored in the regridding control structure. +function get_adapt_CS(CS) + type(regridding_CS), intent(in) :: CS !< Regridding control structure + type(adapt_CS), pointer :: get_adapt_CS + + get_adapt_CS => CS%adapt_CS +end function get_adapt_CS + !------------------------------------------------------------------------------ !> Return coordinate-derived thicknesses for fixed coordinate systems function getStaticThickness( CS, SSH, depth ) diff --git a/src/ALE/coord_adapt.F90 b/src/ALE/coord_adapt.F90 index 3b6a068f66..960befe6d2 100644 --- a/src/ALE/coord_adapt.F90 +++ b/src/ALE/coord_adapt.F90 @@ -5,16 +5,67 @@ !> Regrid columns for the adaptive coordinate module coord_adapt +use MOM_coms, only : reproducing_sum use MOM_EOS, only : calculate_density_derivs -use MOM_error_handler, only : MOM_error, FATAL +use MOM_error_handler, only : MOM_error, FATAL, WARNING use MOM_unit_scaling, only : unit_scale_type use MOM_variables, only : ocean_grid_type, thermo_var_ptrs use MOM_verticalGrid, only : verticalGrid_type +use filter_utils, only : filter_CS, filtered_grid_motion +use coord_zlike, only : init_coord_zlike, end_coord_zlike, zlike_CS, set_zlike_params, build_zstar_column implicit none ; private #include +!> AG coordinate diagnostic control structure +type, public :: adapt_diag_CS + !> Along-coordinate i-gradient of density (used for density term) + real, dimension(:,:,:), allocatable :: slope_u + !> Along-coordinate j-gradient of density (used for density term) + real, dimension(:,:,:), allocatable :: slope_v + + !> Denominator used for calculating density displacement, i-direction + real, dimension(:,:,:), allocatable :: denom_u + !> Denominator used for calculating density displacement, j-direction + real, dimension(:,:,:), allocatable :: denom_v + + !> Physical-space slope of interface along i-direction (used for density weighting) + real, dimension(:,:,:), allocatable :: phys_u + !> Physical-space slope of interface along j-direction (used for density weighting) + real, dimension(:,:,:), allocatable :: phys_v + + !> Coordinate-space slope of interface along i-direction (used for density weighting) + real, dimension(:,:,:), allocatable :: coord_u + !> Coordinate-space slope of interface along j-direction (used for density weighting) + real, dimension(:,:,:), allocatable :: coord_v + + !> Amount of limiting applied to density (before weighting) + real, dimension(:,:,:), allocatable :: limiting_density + !> Amount of limiting applied to smoothing (before weighting) + real, dimension(:,:,:), allocatable :: limiting_smoothing + + !> The adjustment provided by the convective adjustment term + real, dimension(:,:,:), allocatable :: w_adjust + + !> Interface displacement due to density term + real, dimension(:,:,:), allocatable :: disp_density + !> Interface displacement due to smoothing term + real, dimension(:,:,:), allocatable :: disp_smoothing + !> Interface displacement due to unlimited smoothing term + real, dimension(:,:,:), allocatable :: disp_unlimited + + !>@{ Diagnostics IDs + integer :: id_slope_u, id_slope_v + integer :: id_denom_u, id_denom_v + integer :: id_phys_u, id_phys_v + integer :: id_coord_u, id_coord_v + integer :: id_limiting_density, id_limiting_smoothing + integer :: id_w_adjust + integer :: id_disp_density, id_disp_smoothing, id_disp_unlimited + !>@} +end type adapt_diag_CS + !> Control structure for adaptive coordinates (coord_adapt). type, public :: adapt_CS ; private @@ -22,60 +73,93 @@ module coord_adapt integer :: nk !> Nominal near-surface resolution [H ~> m or kg m-2] - real, allocatable, dimension(:) :: coordinateResolution + real, allocatable, dimension(:) :: coordinate_resolution + + !> If positive, a manual coefficient for the density adaptivity term. + !! If negative, either density or pressure adaptivity are chosen, + !! depending on the local coordinate slope, with a minimum of min_smooth + !! going toward the pressure term. + real :: alpha_rho - !> Ratio of optimisation and diffusion timescales [nondim] - real :: adaptTimeRatio + !> The complement of alpha_rho: a positive value is a manually-specified + !! coefficient; a negative value is automatically-determined, with a + !! value of at least min_smooth. + real :: alpha_p - !> Nondimensional coefficient determining how much optimisation to apply [nondim] - real :: adaptAlpha + !> Minimum weighting of the pressure adaptivity (smoothing) term, used + !! when alpha_rho and alpha_p are negative. + real :: min_smooth - !> Near-surface zooming depth [H ~> m or kg m-2] - real :: adaptZoom + !> The timescale over which to apply the diffusive adaptivity terms. [T ~> s] + real :: adaptivity_timescale - !> Near-surface zooming coefficient [nondim] - real :: adaptZoomCoeff + !> The timescale over which to restore towards the calculated + !! or pre-defined target coordinate. [T ~> s] + real :: restoring_timescale - !> Stratification-dependent diffusion coefficient [nondim] - real :: adaptBuoyCoeff + !> Interface slope cutoff for defining stratified/unstratified regions. + real :: slope_cutoff - !> Reference density difference for stratification-dependent diffusion [R ~> kg m-3] - real :: adaptDrho0 + !> If true, use the uniform mean of thicknesses where required. + !! Otherwise, use the "upstream" thickness in the direction of + !! interface movement due to adaptivity. + logical :: use_mean_h - !> If true, form a HYCOM1-like mixed layet by preventing interfaces - !! from becoming shallower than the depths set by coordinateResolution - logical :: adaptDoMin = .false. + !> If true, the on-interface density gradient is calculated in the layers + !! above and below. They must agree on sign to prevent a null mode, and the + !! minimum is chosen, to prefer smoothing. + !! Otherwise, the gradient is calculated directly on the interface. + logical :: use_twin_gradient + + !> If true, calculate the slope in physical space (taking into account the + !! vertical distance between adjacent points). Otherwise, the slope is only + !! calculated along the interface. + logical :: use_physical_slope + + !> If true, restore towards the dynamically-determined mean position of + !! a given interface. Otherwise, use the specified coordinate locations. + logical :: do_restore_mean + + !> The non-dimensional scale for the adjustment performed for diagonal + !! convective instabilities. + real :: adjustment_scale + + !> Used if do_restore_mean is .false.: delegate to a zlike coordinate + !! for the restoring term target. + type(zlike_CS), pointer :: zlike_CS => null() + + !> Used for outputting diagnostics from within the regridding routine. + type(adapt_diag_CS), pointer :: diag_CS => null() end type adapt_CS -public init_coord_adapt, set_adapt_params, build_adapt_column, end_coord_adapt +public init_coord_adapt, set_adapt_params, build_adapt_grid, end_coord_adapt +public associate_adapt_diag, get_adapt_diag_CS contains !> Initialise an adapt_CS with parameters -subroutine init_coord_adapt(CS, nk, coordinateResolution, m_to_H, kg_m3_to_R) +subroutine init_coord_adapt(CS, nk, coordinate_resolution) type(adapt_CS), pointer :: CS !< Unassociated pointer to hold the control structure integer, intent(in) :: nk !< Number of layers in the grid - real, dimension(:), intent(in) :: coordinateResolution !< Nominal near-surface resolution [m] or + real, dimension(:), intent(in) :: coordinate_resolution !< Nominal near-surface resolution [m] or !! other units specified with m_to_H - real, intent(in) :: m_to_H !< A conversion factor from m to the units of thicknesses, - !! perhaps in units of [H m-1 ~> 1 or kg m-3] - real, intent(in) :: kg_m3_to_R !< A conversion factor from kg m-3 to the units of density, - !! perhaps in units of [R m3 kg-1 ~> 1] if (associated(CS)) call MOM_error(FATAL, "init_coord_adapt: CS already associated") allocate(CS) - allocate(CS%coordinateResolution(nk)) + allocate(CS%coordinate_resolution(nk)) CS%nk = nk - CS%coordinateResolution(:) = coordinateResolution(:) + CS%coordinate_resolution(:) = coordinate_resolution(:) + + CS%alpha_rho = -1.0 + CS%alpha_p = -1.0 - ! Set real parameter default values - CS%adaptTimeRatio = 1e-1 ! Nondim. - CS%adaptAlpha = 1.0 ! Nondim. - CS%adaptZoom = 200.0 * m_to_H ! [H ~> m or kg m-2] - CS%adaptZoomCoeff = 0.0 ! Nondim. - CS%adaptBuoyCoeff = 0.0 ! Nondim. - CS%adaptDrho0 = 0.5 * kg_m3_to_R ! [R ~> kg m-3] + CS%use_mean_h = .false. + CS%use_twin_gradient = .true. + CS%use_physical_slope = .true. + CS%do_restore_mean = .false. + + call init_coord_zlike(CS%zlike_CS, nk, coordinate_resolution) end subroutine init_coord_adapt @@ -85,222 +169,863 @@ subroutine end_coord_adapt(CS) ! nothing to do if (.not. associated(CS)) return - deallocate(CS%coordinateResolution) + + call end_coord_zlike(CS%zlike_CS) + + if (associated(CS%diag_CS)) deallocate(CS%diag_CS) + + deallocate(CS%coordinate_resolution) deallocate(CS) end subroutine end_coord_adapt !> This subtroutine can be used to set the parameters for coord_adapt module -subroutine set_adapt_params(CS, adaptTimeRatio, adaptAlpha, adaptZoom, adaptZoomCoeff, & - adaptBuoyCoeff, adaptDrho0, adaptDoMin) +subroutine set_adapt_params(CS, alpha_rho, alpha_p, adaptivity_timescale, use_mean_h, & + use_twin_gradient, slope_cutoff, min_smooth, use_physical_slope, & + restoring_timescale, do_restore_mean, & + adjustment_scale) + type(adapt_CS), pointer :: CS !< The control structure for this module - real, optional, intent(in) :: adaptTimeRatio !< Ratio of optimisation and diffusion timescales [nondim] - real, optional, intent(in) :: adaptAlpha !< Nondimensional coefficient determining - !! how much optimisation to apply [nondim] - real, optional, intent(in) :: adaptZoom !< Near-surface zooming depth [H ~> m or kg m-2] - real, optional, intent(in) :: adaptZoomCoeff !< Near-surface zooming coefficient [nondim] - real, optional, intent(in) :: adaptBuoyCoeff !< Stratification-dependent diffusion coefficient [nondim] - real, optional, intent(in) :: adaptDrho0 !< Reference density difference for - !! stratification-dependent diffusion [R ~> kg m-3] - logical, optional, intent(in) :: adaptDoMin !< If true, form a HYCOM1-like mixed layer by - !! preventing interfaces from becoming shallower than - !! the depths set by coordinateResolution + real, optional, intent(in) :: alpha_rho !< Density adaptivity coefficient + real, optional, intent(in) :: alpha_p !< Pressure adaptivity coefficient + real, optional, intent(in) :: adaptivity_timescale !< Adaptivity timescale + logical, optional, intent(in) :: use_mean_h !< Use uniform or "upstream" mean thickness? + logical, optional, intent(in) :: use_twin_gradient !< Calculate interface density gradient layers above and below + real, optional, intent(in) :: slope_cutoff !< Stratified/unstratified cutoff + real, optional, intent(in) :: min_smooth !< Minimum pressure adaptivity contribution + logical, optional, intent(in) :: use_physical_slope !< Use physical or along-interface slope + real, optional, intent(in) :: restoring_timescale !< Timescale for restoring term + logical, optional, intent(in) :: do_restore_mean !< Restore to the mean height? + real, optional, intent(in) :: adjustment_scale !< Hydrostatic adjustment scale if (.not. associated(CS)) call MOM_error(FATAL, "set_adapt_params: CS not associated") - if (present(adaptTimeRatio)) CS%adaptTimeRatio = adaptTimeRatio - if (present(adaptAlpha)) CS%adaptAlpha = adaptAlpha - if (present(adaptZoom)) CS%adaptZoom = adaptZoom - if (present(adaptZoomCoeff)) CS%adaptZoomCoeff = adaptZoomCoeff - if (present(adaptBuoyCoeff)) CS%adaptBuoyCoeff = adaptBuoyCoeff - if (present(adaptDrho0)) CS%adaptDrho0 = adaptDrho0 - if (present(adaptDoMin)) CS%adaptDoMin = adaptDoMin + if (present(alpha_rho)) CS%alpha_rho = alpha_rho + if (present(alpha_p)) CS%alpha_p = alpha_p + if (present(adaptivity_timescale)) CS%adaptivity_timescale = adaptivity_timescale + if (present(use_mean_h)) CS%use_mean_h = use_mean_h + if (present(use_twin_gradient)) CS%use_twin_gradient = use_twin_gradient + if (present(slope_cutoff)) CS%slope_cutoff = slope_cutoff + if (present(min_smooth)) CS%min_smooth = min_smooth + if (present(use_physical_slope)) CS%use_physical_slope = use_physical_slope + if (present(restoring_timescale)) CS%restoring_timescale = restoring_timescale + if (present(do_restore_mean)) CS%do_restore_mean = do_restore_mean + if (present(adjustment_scale)) CS%adjustment_scale = adjustment_scale end subroutine set_adapt_params -subroutine build_adapt_column(CS, G, GV, US, tv, i, j, zInt, tInt, sInt, h, nom_depth_H, zNext) - type(adapt_CS), intent(in) :: CS !< The control structure for this module +!> Associate a diagnostic control structure with an existing +!! AG control structure -- used to get around the circular +!! dependency of diagnostics depending on coordinates. +subroutine associate_adapt_diag(CS, diag_CS) + type(adapt_CS), pointer :: CS + type(adapt_diag_CS), target :: diag_CS + + if (associated(CS%diag_CS)) deallocate(CS%diag_CS) + CS%diag_CS => diag_CS +end subroutine associate_adapt_diag + +!> Return the associated diagnostic control structure for an +!! AG control structure +function get_adapt_diag_CS(CS) + type(adapt_CS), intent(in) :: CS + type(adapt_diag_CS), pointer :: get_adapt_diag_CS + + get_adapt_diag_CS => CS%diag_CS +end function get_adapt_diag_CS + +!> Calculate the along-coordinate density derivatives +!! and the physical analogue thereof. The derivatives can +!! be calculated in the i- or j-direction, depending on the +!! value of di/dj. +subroutine calc_derivs(G, GV, CS, US, h, z_int, tv, i, j, k, & + di, dj, dk_sig_int, alpha, beta, Idx, mask, hd_sig, hd_sig_phys) + type(ocean_grid_type), intent(in) :: G + type(verticalGrid_type), intent(in) :: GV + type(adapt_CS), intent(in) :: CS + type(unit_scale_type), intent(in) :: US + real, dimension(SZI_(G), SZJ_(G), SZK_(GV)), intent(in) :: h + real, dimension(SZI_(G), SZJ_(G), SZK_(GV)+1), intent(in) :: z_int + type(thermo_var_ptrs), intent(in) :: tv + integer, intent(in) :: i, j, k, di, dj + real, dimension(SZI_(G), SZJ_(G)), intent(in) :: dk_sig_int + real, intent(in) :: alpha, beta, Idx, mask + real, intent(out) :: hd_sig, hd_sig_phys + + real :: H_to_L + real :: d_sig_up, d_sig_dn, d_sig, dk_sig, h_interp + + H_to_L = GV%H_to_Z * US%Z_to_L + + if (CS%use_twin_gradient) then + d_sig_up = alpha * (tv%t(i+di,j+dj,k-1) - tv%t(i,j,k-1)) & + + beta * (tv%s(i+di,j+dj,k-1) - tv%s(i,j,k-1)) + d_sig_dn = alpha * (tv%t(i+di,j+dj,k) - tv%t(i,j,k)) & + + beta * (tv%s(i+di,j+dj,k) - tv%s(i,j,k)) + + if (d_sig_up * d_sig_dn <= 0.) then + d_sig = 0. + else + d_sig = sign(min(abs(d_sig_up), abs(d_sig_dn)), d_sig_up) + end if + end if + + dk_sig = 0.5 * (dk_sig_int(i,j) + dk_sig_int(i+di,j+dj)) + + if (d_sig * dk_sig < 0.) then + h_interp = 0.5 * (h(i,j,k-1) + h(i+di,j+dj,k)) + else + h_interp = 0.5 * (h(i,j,k) + h(i+di,j+dj,k-1)) + end if + + if (CS%use_mean_h) & + h_interp = 0.25 * ((h(i,j,k-1) + h(i+di,j,k)) + (h(i,j,k) + h(i+di,j,k-1))) + + hd_sig = h_interp * d_sig * Idx * H_to_L * mask + hd_sig_phys = hd_sig - Idx * dk_sig * (z_int(i+di,j+dj,K) - z_int(i,j,K)) * H_to_L * mask +end subroutine calc_derivs + +subroutine build_adapt_grid(G, GV, US, h, tv, dzInterface, CS, fCS, min_thickness, dt) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure - type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type - type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various - !! thermodynamic variables - integer, intent(in) :: i !< The i-index of the column to work on - integer, intent(in) :: j !< The j-index of the column to work on - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: zInt !< Interface heights [H ~> m or kg m-2]. - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: tInt !< Interface temperatures [C ~> degC] - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: sInt !< Interface salinities [S ~> ppt] - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] - real, dimension(SZI_(G),SZJ_(G)), intent(in) :: nom_depth_H !< The bathymetric depth of this column - !! relative to mean sea level or another locally - !! valid reference height, converted to thickness - !! units [H ~> m or kg m-2] - real, dimension(SZK_(GV)+1), intent(inout) :: zNext !< updated interface positions [H ~> m or kg m-2] - - ! Local variables - integer :: k, nz - real :: h_up ! The upwind source grid thickness based on the direction of the - ! adjustive fluxes [H ~> m or kg m-2] - real :: b1 ! The inverse of the tridiagonal denominator [nondim] - real :: b_denom_1 ! The leading term in the tridiagonal denominator [nondim] - real :: d1 ! A term in the tridiagonal expressions [nondim] - real :: depth ! Depth in thickness units [H ~> m or kg m-2] - real :: nominal_z ! A nominal interface position in thickness units [H ~> m or kg m-2] - real :: stretching ! A stretching factor for the water column [nondim] - real :: drdz ! The vertical density gradient [R H-1 ~> kg m-4 or m-1] - real, dimension(SZK_(GV)+1) :: alpha ! drho/dT [R C-1 ~> kg m-3 degC-1] - real, dimension(SZK_(GV)+1) :: beta ! drho/dS [R S-1 ~> kg m-3 ppt-1] - real, dimension(SZK_(GV)+1) :: del2sigma ! Laplacian of in situ density times grid spacing [R ~> kg m-3] - real, dimension(SZK_(GV)+1) :: dh_d2s ! Thickness change in response to del2sigma [H ~> m or kg m-2] - real, dimension(SZK_(GV)) :: kGrid ! grid diffusivity on layers [nondim] - real, dimension(SZK_(GV)) :: c1 ! A tridiagonal work array [nondim] - - nz = CS%nk - - ! set bottom and surface of zNext - zNext(1) = 0. - zNext(nz+1) = zInt(i,j,nz+1) - - ! local depth for scaling diffusivity - depth = nom_depth_H(i,j) - - ! initialize del2sigma and the thickness change response to it zero - del2sigma(:) = 0.0 ; dh_d2s(:) = 0.0 - - ! calculate del-squared of neutral density by a - ! stencilled finite difference - ! TODO: this needs to be adjusted to account for vanished layers near topography - - ! up (j-1) - if (G%mask2dT(i,j-1) > 0.0) then - call calculate_density_derivs( & - 0.5 * (tInt(i,j,2:nz) + tInt(i,j-1,2:nz)), & - 0.5 * (sInt(i,j,2:nz) + sInt(i,j-1,2:nz)), & - 0.5 * (zInt(i,j,2:nz) + zInt(i,j-1,2:nz)) * (GV%H_to_RZ * GV%g_Earth), & - alpha, beta, tv%eqn_of_state, (/2,nz/) ) - - del2sigma(2:nz) = del2sigma(2:nz) + & - (alpha(2:nz) * (tInt(i,j-1,2:nz) - tInt(i,j,2:nz)) + & - beta(2:nz) * (sInt(i,j-1,2:nz) - sInt(i,j,2:nz))) - endif - ! down (j+1) - if (G%mask2dT(i,j+1) > 0.0) then - call calculate_density_derivs( & - 0.5 * (tInt(i,j,2:nz) + tInt(i,j+1,2:nz)), & - 0.5 * (sInt(i,j,2:nz) + sInt(i,j+1,2:nz)), & - 0.5 * (zInt(i,j,2:nz) + zInt(i,j+1,2:nz)) * (GV%H_to_RZ * GV%g_Earth), & - alpha, beta, tv%eqn_of_state, (/2,nz/) ) - - del2sigma(2:nz) = del2sigma(2:nz) + & - (alpha(2:nz) * (tInt(i,j+1,2:nz) - tInt(i,j,2:nz)) + & - beta(2:nz) * (sInt(i,j+1,2:nz) - sInt(i,j,2:nz))) - endif - ! left (i-1) - if (G%mask2dT(i-1,j) > 0.0) then - call calculate_density_derivs( & - 0.5 * (tInt(i,j,2:nz) + tInt(i-1,j,2:nz)), & - 0.5 * (sInt(i,j,2:nz) + sInt(i-1,j,2:nz)), & - 0.5 * (zInt(i,j,2:nz) + zInt(i-1,j,2:nz)) * (GV%H_to_RZ * GV%g_Earth), & - alpha, beta, tv%eqn_of_state, (/2,nz/) ) - - del2sigma(2:nz) = del2sigma(2:nz) + & - (alpha(2:nz) * (tInt(i-1,j,2:nz) - tInt(i,j,2:nz)) + & - beta(2:nz) * (sInt(i-1,j,2:nz) - sInt(i,j,2:nz))) - endif - ! right (i+1) - if (G%mask2dT(i+1,j) > 0.0) then - call calculate_density_derivs( & - 0.5 * (tInt(i,j,2:nz) + tInt(i+1,j,2:nz)), & - 0.5 * (sInt(i,j,2:nz) + sInt(i+1,j,2:nz)), & - 0.5 * (zInt(i,j,2:nz) + zInt(i+1,j,2:nz)) * (GV%H_to_RZ * GV%g_Earth), & - alpha, beta, tv%eqn_of_state, (/2,nz/) ) - - del2sigma(2:nz) = del2sigma(2:nz) + & - (alpha(2:nz) * (tInt(i+1,j,2:nz) - tInt(i,j,2:nz)) + & - beta(2:nz) * (sInt(i+1,j,2:nz) - sInt(i,j,2:nz))) - endif - - ! at this point, del2sigma contains the local neutral density curvature at - ! h-points, on interfaces - ! we need to divide by drho/dz to give an interfacial displacement - ! - ! a positive curvature means we're too light relative to adjacent columns, - ! so del2sigma needs to be positive too (push the interface deeper) - call calculate_density_derivs(tInt(i,j,:), sInt(i,j,:), zInt(i,j,:) * (GV%H_to_RZ * GV%g_Earth), & - alpha, beta, tv%eqn_of_state, (/1,nz+1/) ) - do K = 2, nz - ! TODO make lower bound here configurable - dh_d2s(K) = del2sigma(K) * (0.5 * (h(i,j,k-1) + h(i,j,k))) / & - max(alpha(K) * (tv%T(i,j,k) - tv%T(i,j,k-1)) + & - beta(K) * (tv%S(i,j,k) - tv%S(i,j,k-1)), 1e-20*US%kg_m3_to_R) - - ! don't move the interface so far that it would tangle with another - ! interface in the direction we're moving (or exceed a Nyquist limit - ! that could cause oscillations of the interface) - h_up = merge(h(i,j,k), h(i,j,k-1), dh_d2s(K) > 0.) - dh_d2s(K) = 0.5 * CS%adaptAlpha * & - sign(min(abs(del2sigma(K)), 0.5 * h_up), dh_d2s(K)) - - ! update interface positions so we can diffuse them - zNext(K) = zInt(i,j,K) + dh_d2s(K) - enddo + type(unit_scale_type), intent(in) :: US !< The dimensional unit scaling type + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses, in H (usually m or kg m-2) + type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to thermodynamic variables + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: dzInterface !< The interface changes + type(adapt_CS), intent(in) :: CS !< Regridding control structure + type(filter_CS), intent(in) :: fCS !< Filtering control structure + real, intent(in) :: min_thickness !< ALE layer minimum thickness + real, optional, intent(in) :: dt !< The intended timestep over which this + !! regridding operation applies - ! solve diffusivity equation to smooth grid - ! upper diagonal coefficients: -kGrid(2:nz) - ! lower diagonal coefficients: -kGrid(1:nz-1) - ! diagonal coefficients: 1 + (kGrid(1:nz-1) + kGrid(2:nz)) - ! - ! first, calculate the diffusivities within layers - do k = 1, nz - ! calculate the dr bit of drdz - drdz = 0.5 * (alpha(K) + alpha(K+1)) * (tInt(i,j,K+1) - tInt(i,j,K)) + & - 0.5 * (beta(K) + beta(K+1)) * (sInt(i,j,K+1) - sInt(i,j,K)) - ! divide by dz from the new interface positions - drdz = drdz / (zNext(K) - zNext(K+1) + GV%H_subroundoff) - ! don't do weird stuff in unstably-stratified regions - drdz = max(drdz, 0.) - - ! set vertical grid diffusivity - kGrid(k) = (CS%adaptTimeRatio * nz**2 * depth) * & - ( CS%adaptZoomCoeff / (CS%adaptZoom + 0.5*(zNext(K) + zNext(K+1))) + & - (CS%adaptBuoyCoeff * drdz / CS%adaptDrho0) + & - max(1.0 - CS%adaptZoomCoeff - CS%adaptBuoyCoeff, 0.0) / depth) - enddo + ! local variables + integer :: i, j, k, k2, kt, nz ! indices and dimension lengths - ! initial denominator (first diagonal element) - b1 = 1.0 - ! initial Q_1 = 1 - q_1 = 1 - 0/1 - d1 = 1.0 - ! work on all interior interfaces - do K = 2, nz - ! calculate numerator of Q_k - b_denom_1 = 1. + d1 * kGrid(k-1) - ! update denominator for k - b1 = 1.0 / (b_denom_1 + kGrid(k)) - - ! calculate q_k - c1(K) = kGrid(k) * b1 - ! update Q_k = 1 - q_k - d1 = b_denom_1 * b1 - - ! update RHS - zNext(K) = b1 * (zNext(K) + kGrid(k-1)*zNext(K-1)) - enddo - ! final substitution - do K = nz, 2, -1 - zNext(K) = zNext(K) + c1(K)*zNext(K+1) + ! interface heights + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: z_int, z_new, h_int + ! drho/dt and drho/ds on interfaces + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: alpha_int + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: beta_int + ! vertical gradient in sigma + real, dimension(SZI_(G),SZJ_(G)) :: dk_sig_int + ! final change in interface height + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: dz_a, dz_p, dz_r + + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h_upd + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: w + + ! interface position after adaptivity, mean interface position across basin + real, dimension(SZK_(GV)+1) :: z_mean, h_col, z_col, z_upd, dz_col + + ! numerator of density term and upstreamed h + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: hdi_sig, hdi_sig_phys + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: hdj_sig, hdj_sig_phys + ! temporary alpha/beta interpolated to velocity points + real :: alpha, beta + ! some temporary quantities + real :: eps, weight, weight2, h_interp, i_denom, j_denom + ! numerator (and intermediates) of density term before multiplication by h + real :: di_sig, di_sig_up, di_sig_dn + real :: dj_sig, dj_sig_up, dj_sig_dn + ! difference quantities interpolated to other locations + real :: hdi_sig_u, hdj_sig_u, hdi_sig_v, hdj_sig_v, dk_sig_u, dk_sig_v + real :: ts_ratio, slope, phys_slope + real :: global_z_sum, global_h_sum + real :: dz_p_unlim + real :: tmp, dir, CFL + real :: dsig_horiz, dsig_vert_up, dsig_vert_down + real :: H_to_L, L_to_H + + logical :: do_diag + + character(len=11) :: fname + + ! we could probably assume some limit without a specified timestep + if (.not. present(dt)) then + dzInterface(:,:,:) = 0.0 + return + end if + + eps = 1. ; eps = epsilon(eps) + nz = GV%ke + + L_to_H = US%L_to_Z * GV%Z_to_H + + call set_zlike_params(CS%zlike_CS, min_thickness=min_thickness) + + ! zero out diagnostic arrays + do_diag = .true. + if (.not. associated(CS%diag_CS)) then + call MOM_error(WARNING, 'build_adapt_grid expected diag_CS associated') + do_diag = .false. + end if + + if (do_diag) then + if (allocated(CS%diag_CS%phys_u)) CS%diag_CS%phys_u(:,:,:) = 0.0 + if (allocated(CS%diag_CS%phys_v)) CS%diag_CS%phys_v(:,:,:) = 0.0 + if (allocated(CS%diag_CS%slope_u)) CS%diag_CS%slope_u(:,:,:) = 0.0 + if (allocated(CS%diag_CS%slope_v)) CS%diag_CS%slope_v(:,:,:) = 0.0 + if (allocated(CS%diag_CS%denom_u)) CS%diag_CS%denom_u(:,:,:) = 0.0 + if (allocated(CS%diag_CS%denom_v)) CS%diag_CS%denom_v(:,:,:) = 0.0 + if (allocated(CS%diag_CS%coord_u)) CS%diag_CS%coord_u(:,:,:) = 0.0 + if (allocated(CS%diag_CS%coord_v)) CS%diag_CS%coord_v(:,:,:) = 0.0 + if (allocated(CS%diag_CS%limiting_smoothing)) CS%diag_CS%limiting_smoothing(:,:,:) = 0.0 + if (allocated(CS%diag_CS%limiting_density)) CS%diag_CS%limiting_density(:,:,:) = 0.0 + if (allocated(CS%diag_CS%w_adjust)) CS%diag_CS%w_adjust(:,:,:) = 0.0 + if (allocated(CS%diag_CS%disp_density)) CS%diag_CS%disp_density(:,:,:) = 0.0 + if (allocated(CS%diag_CS%disp_smoothing)) CS%diag_CS%disp_smoothing(:,:,:) = 0.0 + if (allocated(CS%diag_CS%disp_unlimited)) CS%diag_CS%disp_unlimited(:,:,:) = 0.0 + end if + + ! sum from free surface downward + z_int(:,:,1) = sum(h, 3) - G%bathyT(:,:) * GV%Z_to_H ! free-surface + do K = 1,nz + z_int(:,:,K+1) = z_int(:,:,K) - h(:,:,k) enddo - if (CS%adaptDoMin) then - nominal_z = 0. - stretching = zInt(i,j,nz+1) / depth + if (CS%do_restore_mean) then + ! calculate geometric mean of thicknesses on interfaces + ! we only need to do this in our own domain because this + ! is a global sum + z_new(:,:,:) = 0. ; h_int(:,:,:) = 0. + do j = G%jsc,G%jec + do i = G%isc,G%iec + h_int(i,j,2:nz) = (h(i,j,2:nz) * h(i,j,1:nz-1)) / & + (h(i,j,2:nz) + h(i,j,1:nz-1) + GV%H_subroundoff) + ! we don't really want to volume-weight this, we just want to discount vanished layers + ! this way, we won't bias towards thick layers + h_int(i,j,2:nz) = max(GV%H_to_m * h_int(i,j,2:nz), 1.0) + h_int(i,j,2:nz) = h_int(i,j,2:nz) * (G%areaT(i,j) * G%mask2dT(i,j)) + ! weight height by thickness + z_new(i,j,2:nz) = z_int(i,j,2:nz) * h_int(i,j,2:nz) + enddo + enddo + global_z_sum = reproducing_sum(z_new, G%isc, G%iec, G%jsc, G%jec, sums=z_mean) + global_h_sum = reproducing_sum(h_int, G%isc, G%iec, G%jsc, G%jec, sums=h_col) + z_mean(2:nz) = z_mean(2:nz) / h_col(2:nz) - do k = 2, nz+1 - nominal_z = nominal_z + CS%coordinateResolution(k-1) * stretching - ! take the deeper of the calculated and nominal positions - zNext(K) = max(zNext(K), nominal_z) - ! interface can't go below topography - zNext(K) = min(zNext(K), zInt(i,j,nz+1)) + do K = 2,nz-1 + if (z_mean(K) < z_mean(K+1)) then + print *, z_mean + call MOM_error(FATAL, 'tangled z_mean') + endif enddo - endif -end subroutine build_adapt_column + else + ! we'll restore to the predefined coordinate resolution + z_mean(1) = 0. + do K = 2,nz + z_mean(K) = z_mean(K-1) - CS%coordinate_resolution(k-1) * GV%Z_to_H + end do + end if + + ! the top and bottom interfaces don't move + dz_a(:,:,1) = 0. ; dz_a(:,:,nz+1) = 0. + dz_p(:,:,1) = 0. ; dz_p(:,:,nz+1) = 0. + dz_r(:,:,1) = 0. ; dz_r(:,:,nz+1) = 0. + w(:,:,1) = 0. ; w(:,:,nz+1) = 0. + + h_upd(:,:,:) = 0. + + ! nondimensionalise the adaptivity timescale wrt. the ALE timescale + ! to get a scaling for diffusive adaptivity + ts_ratio = dt / CS%adaptivity_timescale + ts_ratio = min(ts_ratio, 1.0) + + !$omp parallel default(none) & + !$omp shared(tv, GV, G, CS, US, z_int, h, alpha_int, beta_int) & + !$omp shared(hdi_sig, hdj_sig, hdi_sig_phys, hdj_sig_phys) & + !$omp shared(L_to_H, ts_ratio, dz_a, dz_p, do_diag, eps, nz) & + !$omp private(i, j, k, dk_sig_int, alpha, beta) & + !$omp private(hdi_sig_u, hdj_sig_u, dk_sig_u) & + !$omp private(hdi_sig_v, hdj_sig_v, dk_sig_v) & + !$omp private(i_denom, j_denom, dz_p_unlim, slope, phys_slope) & + !$omp private(weight, weight2) + block + use MOM_domains, only : pass_var, EAST_FACE, NORTH_FACE + + ! for some reason we get a segfault if these are brought in as private to the + ! parallel block, so instead we allocate them locally (they'll be deallocated at the + ! end of the block anyway, but annoying to have to use heap space) + real, allocatable, dimension(:,:) :: t_int, s_int + real, allocatable, dimension(:,:) :: dz_s_i, dz_s_j, dz_p_i, dz_p_j, dz_i, dz_j + real, allocatable, dimension(:,:) :: weight_adapt_i, weight_smooth_i, weight_adapt_j, weight_smooth_j + ! number of active points in stencil, and stencil position + integer :: np, ni, nj + integer, parameter :: filter_width = 3 + + allocate(t_int(SZI_(G),SZJ_(G)), s_int(SZI_(G),SZJ_(G))) + allocate(dz_s_i(SZIB_(G),SZJ_(G)), dz_s_j(SZI_(G),SZJB_(G))) + allocate(dz_p_i(SZIB_(G),SZJ_(G)), dz_p_j(SZI_(G),SZJB_(G))) + allocate(dz_i(SZIB_(G),SZJ_(G)), dz_j(SZI_(G),SZJB_(G))) + allocate(weight_adapt_i(SZIB_(G),SZJ_(G)), weight_smooth_i(SZIB_(G),SZJ_(G))) + allocate(weight_adapt_j(SZI_(G),SZJB_(G)), weight_smooth_j(SZI_(G),SZJB_(G))) + + !$omp do + do K = 2,nz + dz_s_i(:,:) = 0. ; dz_s_j(:,:) = 0. + dz_p_i(:,:) = 0. ; dz_p_j(:,:) = 0. + dz_i(:,:) = 0. ; dz_j(:,:) = 0. + weight_adapt_i(:,:) = 0. ; weight_smooth_i(:,:) = 0. + weight_adapt_j(:,:) = 0. ; weight_smooth_j(:,:) = 0. + + do j = G%jsc-3,G%jec+3 + do i = G%isc-3,G%iec+3 + t_int(i,j) = ( & + tv%t(i,j,k-1) * (h(i,j,k) + GV%H_subroundoff) + & + tv%t(i,j,k) * (h(i,j,k-1) + GV%H_subroundoff)) / & + (h(i,j,k-1) + h(i,j,k) + 2*GV%H_subroundoff) + s_int(i,j) = ( & + tv%s(i,j,k-1) * (h(i,j,k) + GV%H_subroundoff) + & + tv%s(i,j,k) * (h(i,j,k-1) + GV%H_subroundoff)) / & + (h(i,j,k-1) + h(i,j,k) + 2*GV%H_subroundoff) + enddo + + call calculate_density_derivs(t_int(:,j), s_int(:,j), -z_int(:,j,K) * GV%H_to_Pa, & + alpha_int(:,j,K), beta_int(:,j,K), & + G%isc-3, G%iec+3 - (G%isc-3) + 1, tv%eqn_of_state) + + do i = G%isc-3,G%iec+3 + dk_sig_int(i,j) = alpha_int(i,j,K) * (tv%t(i,j,k) - tv%t(i,j,k-1)) + & + beta_int(i,j,K) * (tv%s(i,j,k) - tv%s(i,j,k-1)) + enddo + enddo + + ! calculate horizontal derivatives on i-points + ! reduce I-halo 2 -> 1 + do j = G%jsc-2,G%jec+2 + do I = G%IscB-2,G%IecB+2 + alpha = 0.5 * (alpha_int(i,j,K) + alpha_int(i+1,j,K)) + beta = 0.5 * (beta_int(i,j,K) + beta_int(i+1,j,K)) + + call calc_derivs(G, GV, CS, US, h, z_int, tv, I, j, k, 1, 0, dk_sig_int, alpha, beta, G%IdxCu(I,j), & + G%mask2dCu(I,j), hdi_sig(I,j,K), hdi_sig_phys(I,j,K)) + enddo + enddo + + ! calculate horizontal derivatives on j-points + ! reduce J-halo 2 -> 1 + do J = G%JscB-2,G%JecB+2 + do i = G%isc-2,G%iec+2 + alpha = 0.5 * (alpha_int(i,j,K) + alpha_int(i,j+1,K)) + beta = 0.5 * (beta_int(i,j,K) + beta_int(i,j+1,K)) + + call calc_derivs(G, GV, CS, US, h, z_int, tv, i, J, k, 0, 1, dk_sig_int, alpha, beta, G%IdyCv(i,J), & + G%mask2dCv(i,J), hdj_sig(i,J,K), hdj_sig_phys(i,J,K)) + enddo + enddo + + ! u-points + do j = G%jsc-1,G%jec+1 + do I = G%IscB-1,G%IecB+1 + if (G%mask2dCu(I,j) == 0) then + dz_i(I,j) = 0. + dz_s_i(I,j) = 0. + dz_p_i(I,j) = 0. + cycle + endif + + ! interpolate terms in the denominator onto the u-point + hdi_sig_u = hdi_sig(I,j,K)**2 + hdj_sig_u = 0.25 * ((hdj_sig(i,J,K)**2 + hdj_sig(i+1,J-1,K)**2) + & + (hdj_sig(i+1,J,K)**2 + hdj_sig(i,J-1,K)**2)) + dk_sig_u = 0.5 * (dk_sig_int(i,j)**2 + dk_sig_int(i+1,j)**2) + + i_denom = hdi_sig_u + hdj_sig_u + dk_sig_u + if (abs(i_denom) < eps .or. dk_sig_int(i,j) < 0.0 .or. dk_sig_int(i+1,j) < 0.0) then + ! if gradients in all directions are exactly zero, we don't want any flux + dz_s_i(I,j) = 0. + else + dz_s_i(I,j) = hdi_sig(I,j,K) / sign(sqrt(i_denom), dk_sig_u) + end if + + if (do_diag) then + ! DIAG: slope_u + if (allocated(CS%diag_CS%slope_u)) CS%diag_CS%slope_u(I,j,K) = dz_s_i(I,j) + ! DIAG: denom_u + if (allocated(CS%diag_CS%denom_u)) CS%diag_CS%denom_u(I,j,K) = sqrt(i_denom) + end if + + ! to convert from the density gradient to the flux, flip the sign and multiply by + ! kappa*dt + dz_s_i(I,j) = -dz_s_i(I,j) * G%dxCu(I,j)**2 * ts_ratio * L_to_H**2 + + dz_p_unlim = dz_s_i(I,j) + + ! limit slope based on adjacent layers + ! dz_s_i has opposite sign to hdi_sig + if (dz_s_i(I,j) < 0.) then + ! hdi_sig positive -- left down, right up + dz_s_i(I,j) = max(dz_s_i(I,j), -0.125 * min( & + h(i,j,k) * G%areaT(i,j), & + h(i+1,j,k-1) * G%areaT(i+1,j)) * G%IdyCu(I,j) * L_to_H) + else + ! hdi_sig negative -- left up, right down + dz_s_i(I,j) = min(dz_s_i(I,j), 0.125 * min( & + h(i,j,k-1) * G%areaT(i,j), & + h(i+1,j,k) * G%areaT(i+1,j)) * G%IdyCu(I,j) * L_to_H) + end if + + if (do_diag) then + ! DIAG: limiting_density + ! difference between the unlimited slope flux and the limited, across the participating adjacent cells + if (allocated(CS%diag_CS%limiting_density)) then + CS%diag_CS%limiting_density(i,j,K) = CS%diag_CS%limiting_density(i,j,K) + & + (dz_s_i(I,j) - dz_p_unlim) + CS%diag_CS%limiting_density(i+1,j,K) = CS%diag_CS%limiting_density(i+1,j,K) + & + (dz_s_i(I,j) - dz_p_unlim) + end if + end if + + ! we also calculate the difference in pressure (interface position) + dz_p_i(I,j) = (z_int(i+1,j,K) - z_int(i,j,K)) * G%dxCu(I,j) * ts_ratio * L_to_H + dz_p_unlim = dz_p_i(I,j) + ! dz_p_i positive => left is further down than right + ! => move left up, right down + + if (dz_p_i(I,j) < 0.) then + ! dz_p_i negative -- right up, left down + dz_p_i(I,j) = max(dz_p_i(I,j), -0.125 * min( & + h(i,j,k) * G%areaT(i,j), & + h(i+1,j,k-1) * G%areaT(i+1,j)) * G%IdyCu(I,j) * L_to_H) + else + ! dz_p_i positive -- left up, right down + dz_p_i(I,j) = min(dz_p_i(I,j), 0.125 * min( & + h(i,j,k-1) * G%areaT(i,j), & + h(i+1,j,k) * G%areaT(i+1,j)) * G%IdyCu(I,j) * L_to_H) + end if + + if (do_diag) then + ! DIAG: limiting_smoothing + ! similar to limiting_density, but applied on the pressure (smoothing) term + if (allocated(CS%diag_CS%limiting_smoothing)) then + CS%diag_CS%limiting_smoothing(i,j,K) = CS%diag_CS%limiting_smoothing(i,j,K) + & + (dz_p_i(I,j) - dz_p_unlim) + CS%diag_CS%limiting_smoothing(i+1,j,K) = CS%diag_CS%limiting_smoothing(i+1,j,K) + & + (dz_p_i(I,j) - dz_p_unlim) + end if + end if + + ! calculate and diagnose along-coordinate slope + if (abs(i_denom) < eps .or. dk_sig_int(i,j) < 0.0 .or. dk_sig_int(i+1,j) < 0.0) then + slope = 1.0 + else + slope = (hdi_sig_u + hdj_sig_u) / i_denom + endif + + ! calculate physical slope + hdi_sig_u = hdi_sig_phys(I,j,K)**2 + hdj_sig_u = 0.25 * ((hdj_sig_phys(i,J,K)**2 + hdj_sig_phys(i+1,J-1,K)**2) + & + (hdj_sig_phys(i+1,J,K)**2 + hdj_sig_phys(i,J-1,K)**2)) + i_denom = hdi_sig_u + hdj_sig_u + dk_sig_u + + if (abs(i_denom) < eps .or. dk_sig_int(i,j) < 0.0 .or. dk_sig_int(i+1,j) < 0.0) then + ! unstratified limit + phys_slope = 1.0 + else + phys_slope = (hdi_sig_u + hdj_sig_u) / i_denom + endif + + if (do_diag) then + ! DIAG: coord_u + if (allocated(CS%diag_CS%coord_u)) CS%diag_CS%coord_u(I,j,K) = slope + ! DIAG: phys_u + if (allocated(CS%diag_CS%phys_u)) CS%diag_CS%phys_u(I,j,K) = phys_slope + end if + + ! use physical slope or not? + if (CS%use_physical_slope) slope = phys_slope + + ! calculate weighting between density and pressure terms + ! by a cutoff value on the local normalised stratification + if (slope <= CS%slope_cutoff**2 .and. k > 2) then + weight = 1.0 - CS%min_smooth ; weight2 = 0. + else + weight = 0.0 ; weight2 = 1.0 - CS%min_smooth + endif + + ! override weights if required + if (CS%alpha_rho >= 0.) then + weight = CS%alpha_rho + + if (CS%alpha_p < 0.) then + weight2 = 1.0 - CS%alpha_rho + else + weight2 = CS%alpha_p + endif + else if (CS%alpha_p >= 0.) then + weight2 = CS%alpha_p + weight = 1.0 - CS%alpha_p + endif + + weight_adapt_i(I,j) = weight + weight_smooth_i(I,j) = weight2 + end do + end do + + ! v-points + do J = G%JscB-1,G%JecB+1 + do i = G%isc-1,G%iec+1 + if (G%mask2dCv(i,J) == 0) then + dz_j(i,J) = 0. + dz_s_j(i,J) = 0. + dz_p_j(i,J) = 0. + cycle + endif + + hdj_sig_v = hdj_sig(i,J,K)**2 + hdi_sig_v = 0.25 * ((hdi_sig(I,j,K)**2 + hdi_sig(I-1,j+1,K)**2) + & + (hdi_sig(I,j+1,K)**2 + hdi_sig(I-1,j,K)**2)) + dk_sig_v = 0.5 * (dk_sig_int(i,j)**2 + dk_sig_int(i,j+1)**2) + + j_denom = hdj_sig_v + hdi_sig_v + dk_sig_v + if (abs(j_denom) < eps .or. dk_sig_int(i,j) < 0.0 .or. dk_sig_int(i,j+1) < 0.0) then + dz_s_j(i,J) = 0. + else + dz_s_j(i,J) = hdj_sig(i,J,K) / sign(sqrt(j_denom), dk_sig_v) + end if + + if (do_diag) then + ! DIAG: slope_v + if (allocated(CS%diag_CS%slope_v)) CS%diag_CS%slope_v(i,J,K) = dz_s_j(i,J) + ! DIAG: denom_v + if (allocated(CS%diag_CS%denom_v)) CS%diag_CS%denom_v(i,J,K) = sqrt(j_denom) + end if + + ! dz_s_j beforehand is unitless (ratio of densities) + dz_s_j(i,J) = -dz_s_j(i,J) * G%dyCv(i,J)**2 * ts_ratio * L_to_H**2 + ! dz_s_j is now [m2] + + dz_p_unlim = dz_s_j(i,J) + + ! density limiter + ! dz_s_j [m2] + if (dz_s_j(i,J) < 0.) then + ! hdj_sig positive -- left down, right up + dz_s_j(i,J) = max(dz_s_j(i,J), -0.125 * min( & + h(i,j,k) * G%areaT(i,j), & + h(i,j+1,k-1) * G%areaT(i,j+1)) * G%IdxCv(i,J) * L_to_H) + else + ! hdj_sig negative -- left up, right down + dz_s_j(i,J) = min(dz_s_j(i,J), 0.125 * min( & + h(i,j,k-1) * G%areaT(i,j), & + h(i,j+1,k) * G%areaT(i,j+1)) * G%IdxCv(i,J) * L_to_H) + end if + + if (do_diag) then + ! DIAG: limiting_density + ! see u-point loop for explanation + if (allocated(CS%diag_CS%limiting_density)) then + CS%diag_CS%limiting_density(i,j,K) = CS%diag_CS%limiting_density(i,j,K) + & + (dz_s_j(i,J) - dz_p_unlim) + CS%diag_CS%limiting_density(i,j+1,K) = CS%diag_CS%limiting_density(i,j+1,K) + & + (dz_s_j(i,J) - dz_p_unlim) + end if + end if + + dz_p_j(i,J) = (z_int(i,j+1,K) - z_int(i,j,K)) * G%dyCv(i,J) * ts_ratio * L_to_H + dz_p_unlim = dz_p_j(i,J) + + if (dz_p_j(i,J) < 0.) then + dz_p_j(i,J) = max(dz_p_j(i,J), -0.125 * min( & + h(i,j,k) * G%areaT(i,j), & + h(i,j+1,k-1) * G%areaT(i,j+1)) * G%IdxCv(i,J) * L_to_H) + else + dz_p_j(i,J) = min(dz_p_j(i,J), 0.125 * min( & + h(i,j,k-1) * G%areaT(i,j), & + h(i,j+1,k) * G%areaT(i,j+1)) * G%IdxCv(i,J) * L_to_H) + end if + + if (do_diag) then + ! DIAG: limiting_smoothing + if (allocated(CS%diag_CS%limiting_smoothing)) then + CS%diag_CS%limiting_smoothing(i,j,K) = CS%diag_CS%limiting_smoothing(i,j,K) + & + (dz_p_j(i,J) - dz_p_unlim) + CS%diag_CS%limiting_smoothing(i,j+1,K) = CS%diag_CS%limiting_smoothing(i,j+1,K) + & + (dz_p_j(i,J) - dz_p_unlim) + end if + end if + + ! diagnose along-coordinate slope + if (abs(j_denom) < eps .or. dk_sig_int(i,j) < 0.0 .or. dk_sig_int(i,j+1) < 0.0) then + slope = 1.0 + else + slope = (hdi_sig_v + hdj_sig_v) / j_denom + endif + + hdj_sig_v = hdj_sig_phys(i,J,K)**2 + hdi_sig_v = 0.25 * ((hdi_sig_phys(I,j,K)**2 + hdi_sig_phys(I-1,j+1,K)**2) + & + (hdi_sig_phys(I,j+1,K)**2 + hdi_sig_phys(I-1,j,K)**2)) + j_denom = hdi_sig_v + hdj_sig_v + dk_sig_v + + if (abs(j_denom) < eps .or. dk_sig_int(i,j) < 0.0 .or. dk_sig_int(i,j+1) < 0.0) then + phys_slope = 1.0 + else + phys_slope = (hdi_sig_v + hdj_sig_v) / j_denom + endif + + if (do_diag) then + ! DIAG: coord_v + if (allocated(CS%diag_CS%coord_v)) CS%diag_CS%coord_v(i,J,K) = slope + ! DIAG: phys_v + if (allocated(CS%diag_CS%phys_v)) CS%diag_CS%phys_v(i,J,K) = phys_slope + end if + + if (CS%use_physical_slope) slope = phys_slope + + if (slope <= CS%slope_cutoff**2 .and. k > 2) then + weight = 1.0 - CS%min_smooth ; weight2 = 0. + else + weight = 0.0 ; weight2 = 1.0 - CS%min_smooth + endif + + ! override weights if required + if (CS%alpha_rho >= 0.) then + weight = CS%alpha_rho + + if (CS%alpha_p < 0.) then + weight2 = 1.0 - CS%alpha_rho + else + weight2 = CS%alpha_p + endif + else if (CS%alpha_p >= 0.) then + weight2 = CS%alpha_p + weight = 1.0 - CS%alpha_p + endif + + weight_adapt_j(i,J) = weight + weight_smooth_j(i,J) = weight2 + end do + end do + + call pass_var(weight_adapt_i, G%Domain, position=EAST_FACE) + call pass_var(weight_smooth_i, G%Domain, position=EAST_FACE) + call pass_var(weight_adapt_j, G%Domain, position=NORTH_FACE) + call pass_var(weight_smooth_j, G%Domain, position=NORTH_FACE) + + do j = G%jsc-1,G%jec+1 + do I = G%IscB-1,G%IecB+1 + if (G%mask2dCu(I,j) == 0) cycle + + weight = 0 ; weight2 = 0 ; np = 0 + + do nj = -filter_width,filter_width ; do ni = -filter_width,filter_width + ! filter point is oob or masked: don't add it to our stencil average + if (i+ni < G%IsdB .or. i+ni > G%IedB .or. & + j+nj < G%jsd .or. j+nj > G%jed .or. & + G%mask2dCu(I+ni,j+nj) == 0) cycle + weight = weight + weight_adapt_i(I+ni,j+nj) + weight2 = weight2 + weight_smooth_i(I+ni,j+nj) + np = np + 1 + end do; end do + + dz_s_i(I,j) = dz_s_i(I,j) * weight / np + dz_p_i(I,j) = dz_p_i(I,j) * weight2 / np + + ! combining density and pressure fluxes + ! and re-apply limiter -- with a full cut-off this isn't necessary + dz_i(I,j) = dz_s_i(I,j) + dz_p_i(I,j) + if (dz_i(I,j) < 0.) then + ! hdi_sig positive -- left down, right up + dz_i(I,j) = max(dz_i(I,j), -0.125 * min( & + h(i,j,k) * G%areaT(i,j), & + h(i+1,j,k-1) * G%areaT(i+1,j)) * G%IdyCu(I,j) * L_to_H) + else + ! hdi_sig negative -- left up, right down + dz_i(I,j) = min(dz_i(I,j), 0.125 * min( & + h(i,j,k-1) * G%areaT(i,j), & + h(i+1,j,k) * G%areaT(i+1,j)) * G%IdyCu(I,j) * L_to_H) + end if + end do + end do + + do J = G%JscB-1,G%JecB+1 + do i = G%isc-1,G%iec+1 + if (G%mask2dCv(i,J) == 0) cycle + + weight = 0 ; weight2 = 0 ; np = 0 + + do nj = -filter_width,filter_width ; do ni = -filter_width,filter_width + if (i+ni < G%isd .or. i+ni > G%ied .or. & + j+nj < G%JsdB .or. j+nj > G%JedB .or. & + G%mask2dCv(i+ni,J+nj) == 0) cycle + weight = weight + weight_adapt_j(i+ni,J+nj) + weight2 = weight2 + weight_smooth_j(i+ni,J+nj) + np = np + 1 + end do; end do + + dz_s_j(i,J) = dz_s_j(i,J) * weight / np + dz_p_j(i,J) = dz_p_j(i,J) * weight2 / np + + dz_j(i,J) = dz_s_j(i,J) + dz_p_j(i,J) + if (dz_j(i,J) < 0.) then + ! hdj_sig positive -- left down, right up + dz_j(i,J) = max(dz_j(i,J), -0.125 * min( & + h(i,j,k) * G%areaT(i,j), & + h(i,j+1,k-1) * G%areaT(i,j+1)) * G%IdxCv(i,J) * L_to_H) + else + ! hdj_sig negative -- left up, right down + dz_j(i,J) = min(dz_j(i,J), 0.125 * min( & + h(i,j,k-1) * G%areaT(i,j), & + h(i,j+1,k) * G%areaT(i,j+1)) * G%IdxCv(i,J) * L_to_H) + end if + end do + end do + + do j = G%jsc-1,G%jec+1 + do i = G%isc-1,G%iec+1 + ! prior to this point, dz_a and dz_p should be limited such that they + ! can't cause any tangling. however, they may still lead to some grid-scale + ! checkerboarding, so we reduce by another factor of 2 + dz_a(i,j,K) = 0.25 * G%IareaT(i,j) / L_to_H & + * ((G%dyCu(I,j) * dz_i(I,j) - G%dyCu(I-1,j) * dz_i(I-1,j)) & + + (G%dxCv(i,J) * dz_j(i,J) - G%dxCv(i,J-1) * dz_j(i,J-1))) + + ! apply the change in interface position due to this flux immediately + z_int(i,j,K) = z_int(i,j,K) + dz_a(i,j,K) + end do + end do + + if (do_diag) then + ! DIAG: disp_density + if (allocated(CS%diag_CS%disp_density)) then + do j = G%jsc-1,G%jec+1 + do i = G%isc-1,G%iec+1 + CS%diag_CS%disp_density(i,j,K) = 0.25 * G%IareaT(i,j) / L_to_H & + * ((G%dyCu(I,j) * dz_s_i(I,j) - G%dyCu(I-1,j) * dz_s_i(I-1,j)) & + + (G%dxCv(i,J) * dz_s_j(i,J) - G%dxCv(i,J-1) * dz_s_j(i,J-1))) + end do + end do + end if + ! DIAG: disp_smoothing + if (allocated(CS%diag_CS%disp_smoothing)) then + do j = G%jsc-1,G%jec+1 + do i = G%isc-1,G%iec+1 + CS%diag_CS%disp_smoothing(i,j,K) = 0.25 * G%IareaT(i,j) / L_to_H & + * ((G%dyCu(I,j) * dz_p_i(I,j) - G%dyCu(I-1,j) * dz_p_i(I-1,j)) & + + (G%dxCv(i,J) * dz_p_j(i,J) - G%dxCv(i,J-1) * dz_p_j(i,J-1))) + end do + end do + end if + end if + + ! calculate the z-smoothing fluxes and apply in a second step + ! this lets us use a "barotropic" limiter, which should be much less + ! restrictive than the layer-based one + do j = G%jsc-1,G%jec+1 + do I = G%IscB-1,G%IecB+1 + if (G%mask2dCu(I,j) == 0) then + dz_p_i(I,j) = 0. + cycle + endif + + dz_p_i(I,j) = (z_int(i+1,j,K) - z_int(i,j,K)) * G%dxCu(I,j) * ts_ratio * L_to_H + ! dz_p_i positive => left is further down than right + ! => move left up, right down + + ! XXX this becomes a barotropic limiter + if (dz_p_i(I,j) < 0.) then + ! dz_p_i negative -- right up, left down + dz_p_i(I,j) = max(dz_p_i(I,j), -min( & + (z_int(i,j,K) - z_int(i,j,nz+1)) * G%areaT(i,j), & + (z_int(i+1,j,1) - z_int(i+1,j,K)) * G%areaT(i+1,j)) * G%IdyCu(I,j) * L_to_H) + else + ! dz_p_i positive -- left up, right down + dz_p_i(I,j) = min(dz_p_i(I,j), min( & + (z_int(i,j,1) - z_int(i,j,K)) * G%areaT(i,j), & + (z_int(i+1,j,K) - z_int(i+1,j,nz+1)) * G%areaT(i+1,j)) * G%IdyCu(I,j) * L_to_H) + end if + dz_p_i(I,j) = dz_p_i(I,j) * CS%min_smooth + end do + end do + + do J = G%JscB-1,G%JecB+1 + do i = G%isc-1,G%iec+1 + if (G%mask2dCv(i,J) == 0) then + dz_p_j(i,J) = 0. + cycle + endif + + dz_p_j(i,J) = (z_int(i,j+1,K) - z_int(i,j,K)) * G%dyCv(i,J) * ts_ratio * L_to_H + + if (dz_p_j(i,J) < 0.) then + dz_p_j(i,J) = max(dz_p_j(i,J), -min( & + (z_int(i,j,K) - z_int(i,j,nz+1)) * G%areaT(i,j), & + (z_int(i,j+1,1) - z_int(i,j+1,K)) * G%areaT(i,j+1)) * G%IdxCv(i,J) * L_to_H) + else + dz_p_j(i,J) = min(dz_p_j(i,J), min( & + (z_int(i,j,1) - z_int(i,j,K)) * G%areaT(i,j), & + (z_int(i,j+1,K) - z_int(i,j+1,nz+1)) * G%areaT(i,j+1)) * G%IdxCv(i,J) * L_to_H) + end if + dz_p_j(i,J) = dz_p_j(i,J) * CS%min_smooth + end do + end do + + ! calculate flux due to barotropically-limited smoothing term + do j = G%jsc-1,G%jec+1 + do i = G%isc-1,G%iec+1 + dz_p(i,j,K) = 0.5 * 0.25 * G%IareaT(i,j) / L_to_H & + * ((G%dyCu(I,j) * dz_p_i(I,j) - G%dyCu(I-1,j) * dz_p_i(I-1,j)) & + + (G%dxCv(i,J) * dz_p_j(i,J) - G%dxCv(i,J-1) * dz_p_j(i,J-1))) + end do + end do + end do + !$omp end do + end block + !$omp end parallel + + if (do_diag) then + ! DIAG: disp_unlimited + if (allocated(CS%diag_CS%disp_unlimited)) & + CS%diag_CS%disp_unlimited(:,:,:) = dz_p(:,:,:) + end if + + ts_ratio = dt / CS%restoring_timescale + !$omp parallel do private(z_upd, z_col, i, j, k) + do j = G%jsc-1,G%jec+1 + do i = G%isc-1,G%iec+1 + dzInterface(i,j,:) = 0. + ! for land points, leave interfaecs undisturbed (possibly doesn't matter) + if (G%mask2dT(i,j) == 0) cycle + + ! calculate change in interface position due to restoring term + ! z_int has already been updated by layer-limited fluxes + ! add the barotropically limited flux too + z_upd(:) = z_int(i,j,:) + dz_p(i,j,:) + + if (fCS%depth_of_time_filter_shallow > 0.) then + ! build a z-star column + call build_zstar_column(CS%zlike_CS, G%bathyT(i,j) * GV%Z_to_H, sum(h(i,j,:)), z_mean, zScale=GV%Z_to_H) + + ! filtered_grid_motion will fail if z_upd and z_mean are tangled with each other + ! this basically means that every pair (z_upd(K),z_mean(K)) should be adjacent in a sorted list + ! we can't (shouldn't?) change z_upd, so we can only tweak z_mean to ensure this condition is met + ! restore with depth-dependent profile + z_col(:) = z_mean(:) + + call filtered_grid_motion(fCS, CS%nk, nz, z_upd, z_col, dz_col) + ! dz_col is the additional displacement on top of the interface displacement we already had + dzInterface(i,j,2:nz) = dz_a(i,j,2:nz) + dz_p(i,j,2:nz) + dz_col(2:nz) + else + do K = 2,nz + dz_r(i,j,K) = ts_ratio * (max(min(z_mean(K), z_upd(1)), z_upd(nz+1)) - z_upd(K)) & + / (1.0 + ts_ratio) + + ! using filtered_grid_motion to obtain our dzInterface leads to a loss of precision: + ! we effectively add the depth of the ocean and immediately subtract it out, losing + ! about 4-5 orders of magnitude! + ! instead, we just apply the calculated value directly + ! combine both the layer-limited and barotropically-limited fluxes + dzInterface(i,j,K) = dz_a(i,j,K) + dz_p(i,j,K) + + if (CS%restoring_timescale > 0.) & + dzInterface(i,j,K) = dzInterface(i,j,K) + dz_r(i,j,K) + enddo + endif + + ! update h from previous steps in preparation for adjustment + do k = 1,nz + h_upd(i,j,k) = h(i,j,k) + (dzInterface(i,j,K) - dzInterface(i,j,K+1)) + enddo + enddo + enddo +end subroutine build_adapt_grid end module coord_adapt diff --git a/src/ALE/filter_utils.F90 b/src/ALE/filter_utils.F90 new file mode 100644 index 0000000000..e4d371b3f0 --- /dev/null +++ b/src/ALE/filter_utils.F90 @@ -0,0 +1,182 @@ +!> Contains utility functions for filtering ALE grids +module filter_utils + +use MOM_error_handler, only : MOM_error, FATAL + +implicit none + +!> Coordinate filtering control structure +type, public :: filter_CS + !> Weight given to old coordinate when blending between new and old grids [nondim] + !! Used only below depth_of_time_filter_shallow, with a cubic variation + !! from zero to full effect between depth_of_time_filter_shallow and + !! depth_of_time_filter_deep. + real :: old_grid_weight = 0. + + !> Depth above which no time-filtering of grid is applied [H ~> m or kg m-2] + real :: depth_of_time_filter_shallow = 0. + + !> Depth below which time-filtering of grid is applied at full effect [H ~> m or kg m-2] + real :: depth_of_time_filter_deep = 0. +end type filter_CS + +public filtered_grid_motion + +contains + +!> Returns the change in interface position motion after filtering and +!! assuming the top and bottom interfaces do not move. The filtering is +!! a function of depth, and is applied as the integrated average filtering +!! over the trajectory of the interface. By design, this code can not give +!! tangled interfaces provided that z_old and z_new are not already tangled. +subroutine filtered_grid_motion( CS, nkt, nk, z_old, z_new, dz_g ) + type(filter_CS), intent(in) :: CS !< Regridding control structure + integer, intent(in) :: nkt !< Number of cells in target grid + integer, intent(in) :: nk !< Number of cells in source grid + real, dimension(nk+1), intent(in) :: z_old !< Old grid position [H ~> m or kg m-2] + real, dimension(nkt+1), intent(in) :: z_new !< New grid position before filtering [H ~> m or kg m-2] + real, dimension(nkt+1), intent(inout) :: dz_g !< Change in interface positions + !! including the effects of filtering [H ~> m or kg m-2] + ! Local variables + real :: sgn ! The sign convention for downward [nondim]. + real :: dz_tgt ! The target grid movement of the unfiltered grid [H ~> m or kg m-2] + real :: zr1 ! The old grid position of an interface relative to the surface [H ~> m or kg m-2] + real :: z_old_k ! The corrected position of the old grid [H ~> m or kg m-2] + real :: Aq ! A temporary variable related to the grid weights [nondim] + real :: Bq ! A temporary variable used in the linear term in the quadratic expression for the + ! filtered grid movement [H ~> m or kg m-2] + real :: z0, dz0 ! Together these give the position of an interface relative to a reference hieght + ! that may be adjusted for numerical accuracy in a solver [H ~> m or kg m-2] + real :: F0 ! An estimated grid movement [H ~> m or kg m-2] + real :: zs ! The depth at which the shallow filtering timescale applies [H ~> m or kg m-2] + real :: zd ! The depth at which the deep filtering timescale applies [H ~> m or kg m-2] + real :: dzwt ! The depth range over which the transition in the filtering timescale occurs [H ~> m or kg m-2] + real :: Idzwt ! The Adcroft reciprocal of dzwt [H-1 ~> m-1 or m2 kg-1] + real :: wtd ! The weight given to the new grid when time filtering [nondim] + real :: Iwtd ! The inverse of wtd [nondim] + real :: Int_zs ! A depth integral of the weights in [H ~> m or kg m-2] + real :: Int_zd ! A depth integral of the weights in [H ~> m or kg m-2] + real :: dInt_zs_zd ! The depth integral of the weights between the deep and shallow depths in [H ~> m or kg m-2] +! For debugging: + real, dimension(nk+1) :: z_act ! The final grid positions after the filtered movement [H ~> m or kg m-2] +! real, dimension(nk+1) :: ddz_g_s, ddz_g_d + logical :: debug = .false. + integer :: k + + if ((z_old(nk+1) - z_old(1)) * (z_new(nkt+1) - z_new(1)) < 0.0) then + call MOM_error(FATAL, "filtered_grid_motion: z_old and z_new use different sign conventions.") + elseif ((z_old(nk+1) - z_old(1)) * (z_new(nkt+1) - z_new(1)) == 0.0) then + ! This is a massless column, so do nothing and return. + do k=1,nkt+1 ; dz_g(k) = 0.0 ; enddo ; return + elseif ((z_old(nk+1) - z_old(1)) + (z_new(nkt+1) - z_new(1)) > 0.0) then + sgn = 1.0 + else + sgn = -1.0 + endif + + if (debug) then + do k=2,nkt+1 + if (sgn*(z_new(k)-z_new(k-1)) < -5e-16*(abs(z_new(k))+abs(z_new(k-1))) ) & + call MOM_error(FATAL, "filtered_grid_motion: z_new is tangled.") + enddo + do k=2,nk+1 + if (sgn*(z_old(k)-z_old(k-1)) < -5e-16*(abs(z_old(k))+abs(z_old(k-1))) ) & + call MOM_error(FATAL, "filtered_grid_motion: z_old is tangled.") + enddo + ! ddz_g_s(:) = 0.0 ; ddz_g_d(:) = 0.0 + endif + + zs = CS%depth_of_time_filter_shallow + zd = CS%depth_of_time_filter_deep + wtd = 1.0 - CS%old_grid_weight + Iwtd = 1.0 / wtd + + dzwt = (zd - zs) + Idzwt = 0.0 ; if (abs(zd - zs) > 0.0) Idzwt = 1.0 / (zd - zs) + dInt_zs_zd = 0.5*(1.0 + Iwtd) * (zd - zs) + Aq = 0.5*(Iwtd - 1.0) + + dz_g(1) = 0.0 + z_old_k = z_old(1) + do k = 2,nkt+1 + if (k<=nk+1) z_old_k = z_old(k) ! This allows for virtual z_old interface at bottom of the model + ! zr1 is positive and increases with depth, and dz_tgt is positive downward. + dz_tgt = sgn*(z_new(k) - z_old_k) + zr1 = sgn*(z_old_k - z_old(1)) + + ! First, handle the two simple and common cases that do not pass through + ! the adjustment rate transition zone. + if ((zr1 > zd) .and. (zr1 + wtd * dz_tgt > zd)) then + dz_g(k) = sgn * wtd * dz_tgt + elseif ((zr1 < zs) .and. (zr1 + dz_tgt < zs)) then + dz_g(k) = sgn * dz_tgt + else + ! Find the new value by inverting the equation + ! integral(0 to dz_new) Iwt(z) dz = dz_tgt + ! This is trivial where Iwt is a constant, and agrees with the two limits above. + + ! Take test values at the transition points to figure out which segment + ! the new value will be found in. + if (zr1 >= zd) then + Int_zd = Iwtd*(zd - zr1) + Int_zs = Int_zd - dInt_zs_zd + elseif (zr1 <= zs) then + Int_zs = (zs - zr1) + Int_zd = dInt_zs_zd + (zs - zr1) + else +! Int_zd = (zd - zr1) * (Iwtd + 0.5*(1.0 - Iwtd) * (zd - zr1) / (zd - zs)) + Int_zd = (zd - zr1) * (Iwtd*(0.5*(zd+zr1) - zs) + 0.5*(zd - zr1)) * Idzwt + Int_zs = (zs - zr1) * (0.5*Iwtd * ((zr1 - zs)) + (zd - 0.5*(zr1+zs))) * Idzwt + ! It has been verified that Int_zs = Int_zd - dInt_zs_zd to within roundoff. + endif + + if (dz_tgt >= Int_zd) then ! The new location is in the deep, slow region. + dz_g(k) = sgn * ((zd-zr1) + wtd*(dz_tgt - Int_zd)) + elseif (dz_tgt <= Int_zs) then ! The new location is in the shallow region. + dz_g(k) = sgn * ((zs-zr1) + (dz_tgt - Int_zs)) + else ! We need to solve a quadratic equation for z_new. + ! For accuracy, do the integral from the starting depth or the nearest + ! edge of the transition region. The results with each choice are + ! mathematically equivalent, but differ in roundoff, and this choice + ! should minimize the likelihood of inadvertently overlapping interfaces. + if (zr1 <= zs) then ; dz0 = zs-zr1 ; z0 = zs ; F0 = dz_tgt - Int_zs + elseif (zr1 >= zd) then ; dz0 = zd-zr1 ; z0 = zd ; F0 = dz_tgt - Int_zd + else ; dz0 = 0.0 ; z0 = zr1 ; F0 = dz_tgt ; endif + + Bq = (dzwt + 2.0*Aq*(z0-zs)) + ! Solve the quadratic: Aq*(zn-z0)**2 + Bq*(zn-z0) - F0*dzwt = 0 + ! Note that b>=0, and the two terms in the standard form cancel for the right root. + dz_g(k) = sgn * (dz0 + 2.0*F0*dzwt / (Bq + sqrt(Bq**2 + 4.0*Aq*F0*dzwt) )) + +! if (debug) then +! dz0 = zs-zr1 ; z0 = zs ; F0 = dz_tgt - Int_zs ; Bq = (dzwt + 2.0*Aq*(z0-zs)) +! ddz_g_s(k) = sgn * (dz0 + 2.0*F0*dzwt / (Bq + sqrt(Bq**2 + 4.0*Aq*F0*dzwt) )) - dz_g(k) +! dz0 = zd-zr1 ; z0 = zd ; F0 = dz_tgt - Int_zd ; Bq = (dzwt + 2.0*Aq*(z0-zs)) +! ddz_g_d(k) = sgn * (dz0 + 2.0*F0*dzwt / (Bq + sqrt(Bq**2 + 4.0*Aq*F0*dzwt) )) - dz_g(k) +! +! if (abs(ddz_g_s(k)) > 1e-12*(abs(dz_g(k)) + abs(dz_g(k)+ddz_g_s(k)))) & +! call MOM_error(WARNING, "filtered_grid_motion: Expect z_output to be tangled (sc).") +! if (abs(ddz_g_d(k) - ddz_g_s(k)) > 1e-12*(abs(dz_g(k)+ddz_g_d(k)) + abs(dz_g(k)+ddz_g_s(k)))) & +! call MOM_error(WARNING, "filtered_grid_motion: Expect z_output to be tangled.") +! endif + endif + + endif + enddo + !dz_g(nkt+1) = 0.0 + + if (debug) then + z_old_k = z_old(1) + do k=1,nkt+1 + if (k<=nk+1) z_old_k = z_old(k) ! This allows for virtual z_old interface at bottom of the model + z_act(k) = z_old_k + dz_g(k) + enddo + do k=2,nkt+1 + if (sgn*((z_act(k))-z_act(k-1)) < -1e-15*(abs(z_act(k))+abs(z_act(k-1))) ) & + call MOM_error(FATAL, "filtered_grid_motion: z_output is tangled.") + enddo + endif + +end subroutine filtered_grid_motion + +end module filter_utils diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index e167054cb6..50a7389ba1 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -24,6 +24,7 @@ target_sources(mom6shared PRIVATE ALE/coord_rho.F90 ALE/coord_sigma.F90 ALE/coord_zlike.F90 + ALE/filter_utils.F90 ALE/MOM_ALE.F90 ALE/MOM_hybgen_regrid.F90 ALE/MOM_hybgen_remap.F90 diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index c79b47f5db..04c1541efd 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -59,6 +59,7 @@ module MOM use MOM_ALE, only : ALE_remap_set_h_vel, ALE_remap_set_h_vel_via_dz use MOM_ALE, only : ALE_update_regrid_weights, pre_ALE_diagnostics, ALE_register_diags use MOM_ALE, only : ALE_set_extrap_boundaries +use MOM_ALE, only : ALE_register_coord_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 @@ -1787,10 +1788,10 @@ subroutine ALE_regridding_and_remapping(CS, G, GV, US, u, v, h, tv, dtdia, Time_ call cpu_clock_begin(id_clock_pass) if (associated(tv%T)) & - call create_group_pass(pass_T_S_h, tv%T, G%Domain, To_All+Omit_Corners, halo=1) + call create_group_pass(pass_T_S_h, tv%T, G%Domain, To_All+Omit_Corners, halo=2) if (associated(tv%S)) & - call create_group_pass(pass_T_S_h, tv%S, G%Domain, To_All+Omit_Corners, halo=1) - call create_group_pass(pass_T_S_h, h, G%Domain, To_All+Omit_Corners, halo=1) + call create_group_pass(pass_T_S_h, tv%S, G%Domain, To_All+Omit_Corners, halo=2) + call create_group_pass(pass_T_S_h, h, G%Domain, To_All+Omit_Corners, halo=2) call do_group_pass(pass_T_S_h, G%Domain) call cpu_clock_end(id_clock_pass) @@ -1817,9 +1818,9 @@ subroutine ALE_regridding_and_remapping(CS, G, GV, US, u, v, h, tv, dtdia, Time_ call diag_update_remap_grids(CS%diag) if (use_ice_shelf) then - call ALE_regrid(G, GV, US, h, h_new, dzRegrid, tv, CS%ALE_CSp, CS%frac_shelf_h, PCM_cell) + call ALE_regrid(G, GV, US, h, h_new, dzRegrid, tv, CS%ALE_CSp, CS%frac_shelf_h, PCM_cell, dt=dtdia) else - call ALE_regrid(G, GV, US, h, h_new, dzRegrid, tv, CS%ALE_CSp, PCM_cell=PCM_cell) + call ALE_regrid(G, GV, US, h, h_new, dzRegrid, tv, CS%ALE_CSp, PCM_cell=PCM_cell, dt=dtdia) endif if (showCallTree) call callTree_waypoint("new grid generated") @@ -3627,6 +3628,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, & CS%use_ALE_algorithm, use_KPP) if (CS%use_ALE_algorithm) then call ALE_register_diags(Time, G, GV, US, diag, CS%ALE_CSp) + call ALE_register_coord_diags(Time, G, GV, US, diag, CS%ALE_CSp) endif ! Do any necessary halo updates on any auxiliary variables that have been initialized. diff --git a/src/user/DOME2d_initialization.F90 b/src/user/DOME2d_initialization.F90 index 96cb779eb5..111bc767cc 100644 --- a/src/user/DOME2d_initialization.F90 +++ b/src/user/DOME2d_initialization.F90 @@ -18,6 +18,7 @@ module DOME2d_initialization use regrid_consts, only : coordinateMode, DEFAULT_COORDINATE_MODE use regrid_consts, only : REGRIDDING_LAYER, REGRIDDING_ZSTAR use regrid_consts, only : REGRIDDING_RHO, REGRIDDING_SIGMA +use regrid_consts, only : REGRIDDING_ADAPTIVE implicit none ; private @@ -196,7 +197,7 @@ subroutine DOME2d_initialize_thickness ( h, depth_tot, G, GV, US, param_file, ju ! ! enddo ; enddo - case ( REGRIDDING_ZSTAR ) + case ( REGRIDDING_ZSTAR, REGRIDDING_ADAPTIVE ) do j=js,je ; do i=is,ie eta1D(nz+1) = -depth_tot(i,j) @@ -284,7 +285,7 @@ subroutine DOME2d_initialize_temperature_salinity ( T, S, h, G, GV, US, param_fi select case ( coordinateMode(verticalCoordinate) ) - case ( REGRIDDING_ZSTAR, REGRIDDING_SIGMA ) + case ( REGRIDDING_ZSTAR, REGRIDDING_SIGMA, REGRIDDING_ADAPTIVE ) do j=js,je ; do i=is,ie xi0 = 0.0