MOM_ALE.F90

1! This file is part of MOM6, the Modular Ocean Model version 6.

2! See the LICENSE file for licensing information.

3! SPDX-License-Identifier: Apache-2.0

4

5!> This module contains the main regridding routines.

6!!

7!! Regridding comprises two steps:

8!! 1. Interpolation and creation of a new grid based on target interface

9!! densities (or any other criterion).

10!! 2. Remapping of quantities between old grid and new grid.

11!!

12!! Original module written by Laurent White, 2008.06.09

13module mom_ale

14

15use mom_debugging, only : check_column_integrals

16use mom_diag_mediator, only : register_diag_field, post_data, diag_ctrl

17use mom_diag_mediator, only : time_type, diag_update_remap_grids, query_averaging_enabled

18use mom_domains, only : create_group_pass, do_group_pass, group_pass_type

19use mom_error_handler, only : mom_error, fatal, warning

20use mom_error_handler, only : calltree_showquery

21use mom_error_handler, only : calltree_enter, calltree_leave, calltree_waypoint

22use mom_hybgen_unmix, only : hybgen_unmix, init_hybgen_unmix, end_hybgen_unmix, hybgen_unmix_cs

23use mom_hybgen_regrid, only : hybgen_regrid_cs

24use mom_file_parser, only : get_param, param_file_type, log_param

25use mom_interface_heights,only : find_eta, calc_derived_thermo

26use mom_open_boundary, only : ocean_obc_type, obc_direction_e, obc_direction_w

27use mom_open_boundary, only : obc_direction_n, obc_direction_s

28use mom_regridding, only : initialize_regridding, regridding_main, end_regridding

29use mom_regridding, only : uniformresolution

30use mom_regridding, only : inflate_vanished_layers_old

31use mom_regridding, only : regridding_preadjust_reqs, convective_adjustment

32use mom_regridding, only : set_target_densities_from_gv, set_target_densities

33use mom_regridding, only : regriddingcoordinatemodedoc, default_coordinate_mode

34use mom_regridding, only : regriddinginterpschemedoc, regriddingdefaultinterpscheme

35use mom_regridding, only : regriddingdefaultboundaryextrapolation

36use mom_regridding, only : regriddingdefaultminthickness

37use mom_regridding, only : regridding_cs, set_regrid_params, write_regrid_file

38use mom_regridding, only : getcoordinateinterfaces

39use mom_regridding, only : getcoordinateunits, getcoordinateshortname

40use mom_regridding, only : getstaticthickness

41use mom_remapping, only : initialize_remapping, end_remapping

42use mom_remapping, only : remapping_core_h, remapping_core_w

43use mom_remapping, only : remappingschemesdoc, remappingdefaultscheme

44use mom_remapping, only : interpolate_column, reintegrate_column

45use mom_remapping, only : remapping_cs, dzfromh1h2, remapping_set_param

46use mom_string_functions, only : uppercase, extractword, extract_integer

47use mom_tracer_registry, only : tracer_registry_type, tracer_type, mom_tracer_chkinv

48use mom_unit_scaling, only : unit_scale_type

49use mom_variables, only : ocean_grid_type, thermo_var_ptrs

50use mom_verticalgrid, only : get_thickness_units, verticalgrid_type

51

52!use regrid_consts, only : coordinateMode, DEFAULT_COORDINATE_MODE

53use regrid_consts, only : coordinateunits, coordinatemode, state_dependent

54use regrid_edge_values, only : edge_values_implicit_h4

55use plm_functions, only : plm_reconstruction, plm_boundary_extrapolation

56use plm_functions, only : plm_extrapolate_slope, plm_monotonized_slope, plm_slope_wa

57use ppm_functions, only : ppm_reconstruction, ppm_boundary_extrapolation

58use recon1d_plm_wls, only : plm_wls

59

60implicit none ; private

61#include <MOM_memory.h>

62

63

64!> ALE control structure

65type, public :: ale_cs ; private

66 logical :: remap_uv_using_old_alg !< If true, uses the old "remapping via a delta z"

67 !! method. If False, uses the new method that

68 !! remaps between grids described by h.

69 logical :: partial_cell_vel_remap !< If true, use partial cell thicknesses at velocity points

70 !! that are masked out where they extend below the shallower

71 !! of the neighboring bathymetry for remapping velocity.

72

73 real :: regrid_time_scale !< The time-scale used in blending between the current (old) grid

74 !! and the target (new) grid [T ~> s]

75

76 type(regridding_cs) :: regridcs !< Regridding parameters and work arrays

77 type(remapping_cs) :: remapcs !< Remapping parameters and work arrays

78 type(remapping_cs) :: vel_remapcs !< Remapping parameters for velocities and work arrays

79

80 type(hybgen_unmix_cs), pointer :: hybgen_unmixcs => null() !< Parameters for hybgen remapping

81

82 logical :: use_hybgen_unmix !< If true, use the hybgen unmixing code before regridding

83 logical :: do_conv_adj !< If true, do convective adjustment before regridding

84

85 integer :: nk !< Used only for queries, not directly by this module

86 real :: bbl_h_vel_mask !< The thickness of a bottom boundary layer within which velocities in

87 !! thin layers are zeroed out after remapping, following practice with

88 !! Hybgen remapping, or a negative value to avoid such filtering

89 !! altogether, in [H ~> m or kg m-2].

90 real :: h_vel_mask !< A thickness at velocity points below which near-bottom layers are

91 !! zeroed out after remapping, following the practice with Hybgen

92 !! remapping, or a negative value to avoid such filtering altogether,

93 !! in [H ~> m or kg m-2].

94

95 logical :: remap_after_initialization !< Indicates whether to regrid/remap after initializing the state.

96

97 integer :: answer_date !< The vintage of the expressions and order of arithmetic to use for

98 !! remapping. Values below 20190101 result in the use of older, less

99 !! accurate expressions that were in use at the end of 2018. Higher

100 !! values result in the use of more robust and accurate forms of

101 !! mathematically equivalent expressions.

102

103 logical :: conserve_ke !< Apply a correction to the baroclinic velocity after remapping to

104 !! conserve KE.

105

106 logical :: debug !< If true, write verbose checksums for debugging purposes.

107 logical :: show_call_tree !< For debugging

108

109 ! for diagnostics

110 type(diag_ctrl), pointer :: diag !< structure to regulate output

111 integer, dimension(:), allocatable :: id_tracer_remap_tendency !< diagnostic id

112 integer, dimension(:), allocatable :: id_htracer_remap_tendency !< diagnostic id

113 integer, dimension(:), allocatable :: id_htracer_remap_tendency_2d !< diagnostic id

114 logical, dimension(:), allocatable :: do_tendency_diag !< flag for doing diagnostics

115 integer :: id_dzregrid = -1 !< diagnostic id

116

117 ! diagnostic for fields prior to applying ALE remapping

118 integer :: id_u_preale = -1 !< diagnostic id for zonal velocity before ALE.

119 integer :: id_v_preale = -1 !< diagnostic id for meridional velocity before ALE.

120 integer :: id_h_preale = -1 !< diagnostic id for layer thicknesses before ALE.

121 integer :: id_t_preale = -1 !< diagnostic id for temperatures before ALE.

122 integer :: id_s_preale = -1 !< diagnostic id for salinities before ALE.

123 integer :: id_e_preale = -1 !< diagnostic id for interface heights before ALE.

124 integer :: id_vert_remap_h = -1 !< diagnostic id for layer thicknesses used for remapping

125 integer :: id_vert_remap_h_tendency = -1 !< diagnostic id for layer thickness tendency due to ALE

126 integer :: id_remap_delta_integ_u2 = -1 !< Change in depth-integrated rho0*u**2/2

127 integer :: id_remap_delta_integ_v2 = -1 !< Change in depth-integrated rho0*v**2/2

128

129end type

130

131! Publicly available functions

132public ale_init

133public ale_end

134public ale_regrid

135public ale_offline_inputs

136public ale_regrid_accelerated

137public ale_remap_scalar

138public ale_remap_tracers

139public ale_remap_velocities

140public ale_remap_set_h_vel, ale_remap_set_h_vel_via_dz

141public ale_remap_interface_vals

142public ale_remap_vertex_vals

143public ale_plm_edge_values

144public ts_plm_edge_values

145public ts_ppm_edge_values

146public ts_plm_wls_edge_values

147public adjustgridforintegrity

148public ale_initregridding

149public ale_getcoordinate

150public ale_getcoordinateunits

151public ale_writecoordinatefile

152public ale_updateverticalgridtype

153public ale_initthicknesstocoord

154public ale_update_regrid_weights

155public pre_ale_diagnostics

156public pre_ale_adjustments

157public ale_remap_init_conds

158public ale_register_diags

159public ale_set_extrap_boundaries

160

161! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional

162! consistency testing. These are noted in comments with units like Z, H, L, and T, along with

163! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units

164! vary with the Boussinesq approximation, the Boussinesq variant is given first.

165

166contains

167

168!> This routine is typically called (from initialize_MOM in file MOM.F90)

169!! before the main time integration loop to initialize the regridding stuff.

170!! We read the MOM_input file to register the values of different

171!! regridding/remapping parameters.

172subroutine ale_init( param_file, G, GV, US, max_depth, CS)

173 type(param_file_type), intent(in) :: param_file !< Parameter file

174 type(ocean_grid_type), intent(in) :: g !< Grid structure

175 type(verticalgrid_type), intent(in) :: gv !< Ocean vertical grid structure

176 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type

177 real, intent(in) :: max_depth !< The maximum depth of the ocean [Z ~> m].

178 type(ale_cs), pointer :: cs !< Module control structure

179

180 ! Local variables

181 character(len=40) :: mdl = "MOM_ALE" ! This module's name.

182 character(len=80) :: string, vel_string ! Temporary strings

183 real :: filter_shallow_depth, filter_deep_depth ! Depth ranges of filtering [H ~> m or kg m-2]

184 integer :: default_answer_date ! The default setting for the various ANSWER_DATE flags.

185 logical :: check_reconstruction

186 logical :: check_remapping

187 logical :: force_bounds_in_subcell

188 logical :: local_logical

189 logical :: remap_boundary_extrap

190 logical :: init_boundary_extrap

191 logical :: om4_remap_via_sub_cells

192 type(hybgen_regrid_cs), pointer :: hybgen_regridcs => null() ! Control structure for hybgen regridding

193 ! for sharing parameters.

194 real :: h_neglect, h_neglect_edge ! small thicknesses [H ~> m or kg m-2]

195

196 if (associated(cs)) then

197 call mom_error(warning, "ALE_init called with an associated "// &

198 "control structure.")

199 return

200 endif

201 allocate(cs)

202

203 cs%show_call_tree = calltree_showquery()

204 if (cs%show_call_tree) call calltree_enter("ALE_init(), MOM_ALE.F90")

205

206 call get_param(param_file, mdl, "REMAP_UV_USING_OLD_ALG", cs%remap_uv_using_old_alg, &

207 "If true, uses the old remapping-via-a-delta-z method for "//&

208 "remapping u and v. If false, uses the new method that remaps "//&

209 "between grids described by an old and new thickness.", &

210 default=.false.)

211

212 ! Initialize and configure regridding

213 call ale_initregridding(g, gv, us, max_depth, param_file, mdl, cs%regridCS)

214 call regridding_preadjust_reqs(cs%regridCS, cs%do_conv_adj, cs%use_hybgen_unmix, &

215 hybgen_cs=hybgen_regridcs)

216

217 ! Initialize and configure remapping that is orchestrated by ALE.

218 call get_param(param_file, mdl, "REMAPPING_SCHEME", string, &

219 "This sets the reconstruction scheme used "//&

220 "for vertical remapping for all variables. "//&

221 "It can be one of the following schemes: \n"//&

222 trim(remappingschemesdoc), default=remappingdefaultscheme)

223 call get_param(param_file, mdl, "VELOCITY_REMAPPING_SCHEME", vel_string, &

224 "This sets the reconstruction scheme used for vertical remapping "//&

225 "of velocities. By default it is the same as REMAPPING_SCHEME. "//&

226 "It can be one of the following schemes: \n"//&

227 trim(remappingschemesdoc), default=trim(string))

228 call get_param(param_file, mdl, "FATAL_CHECK_RECONSTRUCTIONS", check_reconstruction, &

229 "If true, cell-by-cell reconstructions are checked for "//&

230 "consistency and if non-monotonicity or an inconsistency is "//&

231 "detected then a FATAL error is issued.", default=.false.)

232 call get_param(param_file, mdl, "FATAL_CHECK_REMAPPING", check_remapping, &

233 "If true, the results of remapping are checked for "//&

234 "conservation and new extrema and if an inconsistency is "//&

235 "detected then a FATAL error is issued.", default=.false.)

236 call get_param(param_file, mdl, "REMAP_BOUND_INTERMEDIATE_VALUES", force_bounds_in_subcell, &

237 "If true, the values on the intermediate grid used for remapping "//&

238 "are forced to be bounded, which might not be the case due to "//&

239 "round off.", default=.false.)

240 call get_param(param_file, mdl, "REMAP_BOUNDARY_EXTRAP", remap_boundary_extrap, &

241 "If true, values at the interfaces of boundary cells are "//&

242 "extrapolated instead of piecewise constant", default=.false.)

243 call get_param(param_file, mdl, "INIT_BOUNDARY_EXTRAP", init_boundary_extrap, &

244 "If true, values at the interfaces of boundary cells are "//&

245 "extrapolated instead of piecewise constant during initialization. "//&

246 "Defaults to REMAP_BOUNDARY_EXTRAP.", default=remap_boundary_extrap)

247 call get_param(param_file, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, &

248 "This sets the default value for the various _ANSWER_DATE parameters.", &

249 default=99991231)

250 call get_param(param_file, mdl, "REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, &

251 "This selects the remapping algorithm used in OM4 that does not use "//&

252 "the full reconstruction for the top- and lower-most sub-layers, but instead "//&

253 "assumes they are always vanished (untrue) and so just uses their edge values. "//&

254 "We recommend setting this option to false.", default=.true.)

255 call get_param(param_file, mdl, "REMAPPING_ANSWER_DATE", cs%answer_date, &

256 "The vintage of the expressions and order of arithmetic to use for remapping. "//&

257 "Values below 20190101 result in the use of older, less accurate expressions "//&

258 "that were in use at the end of 2018. Higher values result in the use of more "//&

259 "robust and accurate forms of mathematically equivalent expressions.", &

260 default=default_answer_date, do_not_log=.not.gv%Boussinesq)

261 if (.not.gv%Boussinesq) cs%answer_date = max(cs%answer_date, 20230701)

262

263 if (cs%answer_date >= 20190101) then

264 h_neglect = gv%H_subroundoff ; h_neglect_edge = gv%H_subroundoff

265 elseif (gv%Boussinesq) then

266 h_neglect = gv%m_to_H * 1.0e-30 ; h_neglect_edge = gv%m_to_H * 1.0e-10

267 else

268 h_neglect = gv%kg_m2_to_H * 1.0e-30 ; h_neglect_edge = gv%kg_m2_to_H * 1.0e-10

269 endif

270

271 call initialize_remapping( cs%remapCS, string, nk=gv%ke, &

272 boundary_extrapolation=init_boundary_extrap, &

273 check_reconstruction=check_reconstruction, &

274 check_remapping=check_remapping, &

275 force_bounds_in_subcell=force_bounds_in_subcell, &

276 om4_remap_via_sub_cells=om4_remap_via_sub_cells, &

277 answer_date=cs%answer_date, &

278 h_neglect=h_neglect, h_neglect_edge=h_neglect_edge)

279 call initialize_remapping( cs%vel_remapCS, vel_string, nk=gv%ke, &

280 boundary_extrapolation=init_boundary_extrap, &

281 check_reconstruction=check_reconstruction, &

282 check_remapping=check_remapping, &

283 force_bounds_in_subcell=force_bounds_in_subcell, &

284 om4_remap_via_sub_cells=om4_remap_via_sub_cells, &

285 answer_date=cs%answer_date, &

286 h_neglect=h_neglect, h_neglect_edge=h_neglect_edge)

287

288 call get_param(param_file, mdl, "PARTIAL_CELL_VELOCITY_REMAP", cs%partial_cell_vel_remap, &

289 "If true, use partial cell thicknesses at velocity points that are masked out "//&

290 "where they extend below the shallower of the neighboring bathymetry for "//&

291 "remapping velocity.", default=.false.)

292

293 call get_param(param_file, mdl, "REMAP_AFTER_INITIALIZATION", cs%remap_after_initialization, &

294 "If true, applies regridding and remapping immediately after "//&

295 "initialization so that the state is ALE consistent. This is a "//&

296 "legacy step and should not be needed if the initialization is "//&

297 "consistent with the coordinate mode.", default=.true.)

298

299 call get_param(param_file, mdl, "REGRID_USE_DEPTH_BASED_TIME_FILTER", local_logical, &

300 "If true, always uses depth-based time filtering code that updates the "//&

301 "generated grid using REGRID_TIME_SCALE, REGRID_FILTER_SHALLOW_DEPTH, "//&

302 "REGRID_FILTER_DEEP_DEPTH parameters. Setting to True always uses "//&

303 "filtering but setting to False bypasses calculations when filter times = 0.", &

304 default=.true.)

305 call set_regrid_params(cs%regridCS, use_depth_based_time_filter=local_logical)

306 call get_param(param_file, mdl, "REGRID_TIME_SCALE", cs%regrid_time_scale, &

307 "The time-scale used in blending between the current (old) grid "//&

308 "and the target (new) grid. A short time-scale favors the target "//&

309 "grid (0. or anything less than DT_THERM) has no memory of the old "//&

310 "grid. A very long time-scale makes the model more Lagrangian.", &

311 units="s", default=0., scale=us%s_to_T)

312 call get_param(param_file, mdl, "REGRID_FILTER_SHALLOW_DEPTH", filter_shallow_depth, &

313 "The depth above which no time-filtering is applied. Above this depth "//&

314 "final grid exactly matches the target (new) grid.", &

315 units="m", default=0., scale=gv%m_to_H)

316 call get_param(param_file, mdl, "REGRID_FILTER_DEEP_DEPTH", filter_deep_depth, &

317 "The depth below which full time-filtering is applied with time-scale "//&

318 "REGRID_TIME_SCALE. Between depths REGRID_FILTER_SHALLOW_DEPTH and "//&

319 "REGRID_FILTER_DEEP_DEPTH the filter weights adopt a cubic profile.", &

320 units="m", default=0., scale=gv%m_to_H)

321 call set_regrid_params(cs%regridCS, depth_of_time_filter_shallow=filter_shallow_depth, &

322 depth_of_time_filter_deep=filter_deep_depth)

323 call get_param(param_file, mdl, "REGRID_USE_OLD_DIRECTION", local_logical, &

324 "If true, the regridding integrates upwards from the bottom for "//&

325 "interface positions, much as the main model does. If false "//&

326 "regridding integrates downward, consistent with the remapping code.", &

327 default=.true., do_not_log=.true.)

328 call set_regrid_params(cs%regridCS, integrate_downward_for_e=.not.local_logical)

329

330 call get_param(param_file, mdl, "REMAP_VEL_MASK_BBL_THICK", cs%BBL_h_vel_mask, &

331 "A thickness of a bottom boundary layer below which velocities in thin layers "//&

332 "are zeroed out after remapping, following practice with Hybgen remapping, "//&

333 "or a negative value to avoid such filtering altogether.", &

334 default=-0.001, units="m", scale=gv%m_to_H)

335 call get_param(param_file, mdl, "REMAP_VEL_MASK_H_THIN", cs%h_vel_mask, &

336 "A thickness at velocity points below which near-bottom layers are zeroed out "//&

337 "after remapping, following practice with Hybgen remapping, "//&

338 "or a negative value to avoid such filtering altogether.", &

339 default=1.0e-6, units="m", scale=gv%m_to_H, do_not_log=(cs%BBL_h_vel_mask<=0.0))

340

341 if (cs%use_hybgen_unmix) &

342 call init_hybgen_unmix(cs%hybgen_unmixCS, gv, us, param_file, hybgen_regridcs)

343

344 call get_param(param_file, mdl, "REMAP_VEL_CONSERVE_KE", cs%conserve_ke, &

345 "If true, a correction is applied to the baroclinic component of velocity "//&

346 "after remapping so that total KE is conserved. KE may not be conserved "//&

347 .and."when (CS%BBL_h_vel_mask > 0.0) (CS%h_vel_mask > 0.0)", &

348 default=.false.)

349 call get_param(param_file, "MOM", "DEBUG", cs%debug, &

350 "If true, write out verbose debugging data.", &

351 default=.false., debuggingparam=.true.)

352

353 ! Keep a record of values for subsequent queries

354 cs%nk = gv%ke

355

356 if (cs%show_call_tree) call calltree_leave("ALE_init()")

357end subroutine ale_init

358

359!> Sets the boundary extrapolation set for the remapping type.

360subroutine ale_set_extrap_boundaries( param_file, CS)

361 type(param_file_type), intent(in) :: param_file !< Parameter file

362 type(ale_cs), pointer :: cs !< Module control structure

363

364 logical :: remap_boundary_extrap

365 call get_param(param_file, "MOM_ALE", "REMAP_BOUNDARY_EXTRAP", remap_boundary_extrap, &

366 "If true, values at the interfaces of boundary cells are "//&

367 "extrapolated instead of piecewise constant", default=.false.)

368 call remapping_set_param(cs%remapCS, boundary_extrapolation=remap_boundary_extrap)

369end subroutine ale_set_extrap_boundaries

370

371!> Sets the remapping algorithm to that of OM4

372!!

373!! The remapping aglorithm used in OM4 made poor assumptions about the reconstructions

374!! in the top/bottom layers, namely that they were always vanished and could be

375!! represented solely by their upper/lower edge value respectively.

376!! Passing .false. here uses the full reconstruction of those top and bottom layers

377!! and properly sample those layers.

378subroutine ale_set_om4_remap_algorithm( CS, om4_remap_via_sub_cells )

379 type(ale_cs), pointer :: CS !< Module control structure

380 logical, intent(in) :: om4_remap_via_sub_cells !< If true, use OM4 remapping algorithm

381

382 call remapping_set_param(cs%remapCS, om4_remap_via_sub_cells=om4_remap_via_sub_cells )

383

384end subroutine ale_set_om4_remap_algorithm

385

386!> Initialize diagnostics for the ALE module.

387subroutine ale_register_diags(Time, G, GV, US, diag, CS)

388 type(time_type),target, intent(in) :: time !< Time structure

389 type(ocean_grid_type), intent(in) :: g !< Grid structure

390 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type

391 type(verticalgrid_type), intent(in) :: gv !< Ocean vertical grid structure

392 type(diag_ctrl), target, intent(in) :: diag !< Diagnostics control structure

393 type(ale_cs), pointer :: cs !< Module control structure

394

395 ! Local variables

396 character(len=48) :: thickness_units

397

398 cs%diag => diag

399 thickness_units = get_thickness_units(gv)

400

401 ! These diagnostics of the state variables before ALE are useful for

402 ! debugging the ALE code.

403 cs%id_u_preale = register_diag_field('ocean_model', 'u_preale', diag%axesCuL, time, &

404 'Zonal velocity before remapping', 'm s-1', conversion=us%L_T_to_m_s)

405 cs%id_v_preale = register_diag_field('ocean_model', 'v_preale', diag%axesCvL, time, &

406 'Meridional velocity before remapping', 'm s-1', conversion=us%L_T_to_m_s)

407 cs%id_h_preale = register_diag_field('ocean_model', 'h_preale', diag%axesTL, time, &

408 'Layer Thickness before remapping', thickness_units, conversion=gv%H_to_MKS, &

409 v_extensive=.true.)

410 cs%id_T_preale = register_diag_field('ocean_model', 'T_preale', diag%axesTL, time, &

411 'Temperature before remapping', 'degC', conversion=us%C_to_degC)

412 cs%id_S_preale = register_diag_field('ocean_model', 'S_preale', diag%axesTL, time, &

413 'Salinity before remapping', 'PSU', conversion=us%S_to_ppt)

414 cs%id_e_preale = register_diag_field('ocean_model', 'e_preale', diag%axesTi, time, &

415 'Interface Heights before remapping', 'm', conversion=us%Z_to_m)

416

417 cs%id_dzRegrid = register_diag_field('ocean_model', 'dzRegrid', diag%axesTi, time, &

418 'Change in interface height due to ALE regridding', 'm', conversion=gv%H_to_m)

419 cs%id_vert_remap_h = register_diag_field('ocean_model', 'vert_remap_h', diag%axestl, time, &

420 'layer thicknesses after ALE regridding and remapping', &

421 thickness_units, conversion=gv%H_to_MKS, v_extensive=.true.)

422 cs%id_vert_remap_h_tendency = register_diag_field('ocean_model', &

423 'vert_remap_h_tendency', diag%axestl, time, &

424 'Layer thicknesses tendency due to ALE regridding and remapping', &

425 trim(thickness_units)//" s-1", conversion=gv%H_to_MKS*us%s_to_T, v_extensive=.true.)

426 cs%id_remap_delta_integ_u2 = register_diag_field('ocean_model', 'ale_u2', diag%axesCu1, time, &

427 'Rate of change in half rho0 times depth integral of squared zonal '//&

428 'velocity by remapping. If REMAP_VEL_CONSERVE_KE is .true. then '//&

429 'this measures the change before the KE-conserving correction is applied.', &

430 'W m-2', conversion=us%RZ3_T3_to_W_m2*gv%H_to_RZ*us%L_to_Z**2)

431 cs%id_remap_delta_integ_v2 = register_diag_field('ocean_model', 'ale_v2', diag%axesCv1, time, &

432 'Rate of change in half rho0 times depth integral of squared meridional '//&

433 'velocity by remapping. If REMAP_VEL_CONSERVE_KE is .true. then '//&

434 'this measures the change before the KE-conserving correction is applied.', &

435 'W m-2', conversion=us%RZ3_T3_to_W_m2*gv%H_to_RZ*us%L_to_Z**2)

436

437end subroutine ale_register_diags

438

439!> Crudely adjust (initial) grid for integrity.

440!! This routine is typically called (from initialize_MOM in file MOM.F90)

441!! before the main time integration loop to initialize the regridding stuff.

442!! We read the MOM_input file to register the values of different

443!! regridding/remapping parameters.

444subroutine adjustgridforintegrity( CS, G, GV, h )

445 type(ale_cs), intent(in) :: cs !< Regridding parameters and options

446 type(ocean_grid_type), intent(in) :: g !< Ocean grid informations

447 type(verticalgrid_type), intent(in) :: gv !< Ocean vertical grid structure

448 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h !< Current 3D grid thickness that

449 !! are to be adjusted [H ~> m or kg m-2]

450 call inflate_vanished_layers_old( cs%regridCS, g, gv, h(:,:,:) )

451

452end subroutine adjustgridforintegrity

453

454

455!> End of regridding (memory deallocation).

456!! This routine is typically called (from MOM_end in file MOM.F90)

457!! after the main time integration loop to deallocate the regridding stuff.

458subroutine ale_end(CS)

459 type(ale_cs), pointer :: cs !< module control structure

460

461 ! Deallocate memory used for the regridding

462 call end_remapping( cs%remapCS )

463

464 if (cs%use_hybgen_unmix) call end_hybgen_unmix( cs%hybgen_unmixCS )

465 call end_regridding( cs%regridCS )

466

467 deallocate(cs)

468

469end subroutine ale_end

470

471!> Save any diagnostics of the state before ALE remapping. These diagnostics are

472!! mostly used for debugging.

473subroutine pre_ale_diagnostics(G, GV, US, h, u, v, tv, CS)

474 type(ocean_grid_type), intent(in) :: g !< Ocean grid informations

475 type(verticalgrid_type), intent(in) :: gv !< Ocean vertical grid structure

476 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type

477 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h !< Current 3D grid obtained after the

478 !! last time step [H ~> m or kg m-2]

479 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(inout) :: u !< Zonal velocity field [L T-1 ~> m s-1]

480 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(inout) :: v !< Meridional velocity field [L T-1 ~> m s-1]

481 type(thermo_var_ptrs), intent(inout) :: tv !< Thermodynamic variable structure

482 type(ale_cs), pointer :: cs !< Regridding parameters and options

483

484 ! Local variables

485 real :: eta_preale(szi_(g),szj_(g),szk_(gv)+1) ! Interface heights before remapping [Z ~> m]

486

487 if (cs%id_u_preale > 0) call post_data(cs%id_u_preale, u, cs%diag)

488 if (cs%id_v_preale > 0) call post_data(cs%id_v_preale, v, cs%diag)

489 if (cs%id_h_preale > 0) call post_data(cs%id_h_preale, h, cs%diag)

490 if (cs%id_T_preale > 0) call post_data(cs%id_T_preale, tv%T, cs%diag)

491 if (cs%id_S_preale > 0) call post_data(cs%id_S_preale, tv%S, cs%diag)

492 if (cs%id_e_preale > 0) then

493 call find_eta(h, tv, g, gv, us, eta_preale, dzref=g%Z_ref)

494 call post_data(cs%id_e_preale, eta_preale, cs%diag)

495 endif

496

497end subroutine pre_ale_diagnostics

498

499

500!> Potentially do some preparatory work, such as convective adjustment, to clean up the model

501!! state before regridding.

502subroutine pre_ale_adjustments(G, GV, US, h, tv, Reg, CS, u, v)

503 type(ocean_grid_type), intent(in) :: g !< Ocean grid informations

504 type(verticalgrid_type), intent(in) :: gv !< Ocean vertical grid structure

505 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type

506 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h !< Current 3D grid obtained after the

507 !! last time step [H ~> m or kg m-2]

508 type(thermo_var_ptrs), intent(inout) :: tv !< Thermodynamic variable structure

509 type(tracer_registry_type), pointer :: reg !< Tracer registry structure

510 type(ale_cs), pointer :: cs !< Regridding parameters and options

511 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &

512 optional, intent(inout) :: u !< Zonal velocity field [L T-1 ~> m s-1]

513 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &

514 optional, intent(inout) :: v !< Meridional velocity field [L T-1 ~> m s-1]

515

516 integer :: ntr

517

518 ! Do column-wise convective adjustment.

519 ! Tracers and velocities should probably also undergo consistent adjustments.

520 if (cs%do_conv_adj) call convective_adjustment(g, gv, h, tv)

521

522 if (cs%use_hybgen_unmix) then

523 ntr = 0 ; if (associated(reg)) ntr = reg%ntr

524 call hybgen_unmix(g, gv, us, cs%hybgen_unmixCS, tv, reg, ntr, h)

525 endif

526

527end subroutine pre_ale_adjustments

528

529!> Takes care of building a new grid. The creation of the new grid can be based on z coordinates,

530!! target interface densities, sigma coordinates or any arbitrary coordinate system.

531subroutine ale_regrid( G, GV, US, h, h_new, dzRegrid, tv, CS, frac_shelf_h, PCM_cell)

532 type(ocean_grid_type), intent(in) :: g !< Ocean grid informations

533 type(verticalgrid_type), intent(in) :: gv !< Ocean vertical grid structure

534 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type

535 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses in 3D grid before

536 !! regridding [H ~> m or kg m-2]

537 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: h_new !< Layer thicknesses in 3D grid after

538 !! regridding [H ~> m or kg m-2]

539 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(out) :: dzregrid !< The change in grid interface positions

540 !! due to regridding, in the same units as

541 !! thicknesses [H ~> m or kg m-2]

542 type(thermo_var_ptrs), intent(inout) :: tv !< Thermodynamic variable structure

543 type(ale_cs), pointer :: cs !< Regridding parameters and options

544 real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: frac_shelf_h !< Fractional ice shelf coverage [nondim]

545 logical, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &

546 optional, intent(out) :: pcm_cell !< If true, use PCM remapping in a cell.

547

548 ! Local variables

549 logical :: showcalltree

550

551 showcalltree = calltree_showquery()

552

553 if (showcalltree) call calltree_enter("ALE_regrid(), MOM_ALE.F90")

554

555 ! Build the new grid and store it in h_new. The old grid is retained as h.

556 ! Both are needed for the subsequent remapping of variables.

557 dzregrid(:,:,:) = 0.0

558 call regridding_main( cs%remapCS, cs%regridCS, g, gv, us, h, tv, h_new, dzregrid, &

559 frac_shelf_h=frac_shelf_h, pcm_cell=pcm_cell)

560

561 if (cs%id_dzRegrid>0) then ; if (query_averaging_enabled(cs%diag)) then

562 call post_data(cs%id_dzRegrid, dzregrid, cs%diag, alt_h=h_new)

563 endif ; endif

564

565 if (showcalltree) call calltree_leave("ALE_regrid()")

566

567end subroutine ale_regrid

568

569!> Regrid/remap stored fields used for offline tracer integrations. These input fields are assumed to have

570!! the same layer thicknesses at the end of the last offline interval (which should be a Zstar grid). This

571!! routine builds a grid on the runtime specified vertical coordinate

572subroutine ale_offline_inputs(CS, G, GV, US, h, tv, Reg, uhtr, vhtr, Kd, debug, OBC)

573 type(ale_cs), pointer :: cs !< Regridding parameters and options

574 type(ocean_grid_type), intent(in ) :: g !< Ocean grid informations

575 type(verticalgrid_type), intent(in ) :: gv !< Ocean vertical grid structure

576 type(unit_scale_type), intent(in ) :: us !< A dimensional unit scaling type

577 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h !< Layer thicknesses [H ~> m or kg m-2]

578 type(thermo_var_ptrs), intent(inout) :: tv !< Thermodynamic variable structure

579 type(tracer_registry_type), pointer :: reg !< Tracer registry structure

580 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(inout) :: uhtr !< Zonal mass fluxes [H L2 ~> m3 or kg]

581 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(inout) :: vhtr !< Meridional mass fluxes [H L2 ~> m3 or kg]

582 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: kd !< Input diffusivities

583 !! [H Z T-1 ~> m2 s-1 or kg m-1 s-1]

584 logical, intent(in ) :: debug !< If true, then turn checksums

585 type(ocean_obc_type), pointer :: obc !< Open boundary structure

586 ! Local variables

587 integer :: nk, i, j, k, isc, iec, jsc, jec

588 real, dimension(SZI_(G), SZJ_(G), SZK_(GV)) :: h_new ! Layer thicknesses after regridding [H ~> m or kg m-2]

589 real, dimension(SZI_(G), SZJ_(G), SZK_(GV)+1) :: dzregrid ! The change in grid interface positions [H ~> m or kg m-2]

590 real, dimension(SZK_(GV)) :: h_src ! Source grid thicknesses at velocity points [H ~> m or kg m-2]

591 real, dimension(SZK_(GV)) :: h_dest ! Destination grid thicknesses at velocity points [H ~> m or kg m-2]

592 real, dimension(SZK_(GV)) :: temp_vec ! Transports on the destination grid [H L2 ~> m3 or kg]

593

594 isc = g%isc ; iec = g%iec ; jsc = g%jsc ; jec = g%jec ; nk = gv%ke

595 dzregrid(:,:,:) = 0.0

596 h_new(:,:,:) = 0.0

597

598 if (debug) call mom_tracer_chkinv("Before ALE_offline_inputs", g, gv, h, reg%Tr, reg%ntr)

599

600 ! Build new grid from the Zstar state onto the requested vertical coordinate. The new grid is stored

601 ! in h_new. The old grid is h. Both are needed for the subsequent remapping of variables. Convective

602 ! adjustment right now is not used because it is unclear what to do with vanished layers

603 call regridding_main( cs%remapCS, cs%regridCS, g, gv, us, h, tv, h_new, dzregrid)

604 if (cs%show_call_tree) call calltree_waypoint("new grid generated (ALE_offline_inputs)")

605

606 ! Remap all variables from old grid h onto new grid h_new

607 call ale_remap_tracers(cs, g, gv, h, h_new, reg, debug=cs%show_call_tree)

608 if (allocated(tv%SpV_avg)) tv%valid_SpV_halo = -1 ! Record that SpV_avg is no longer valid.

609 if (cs%show_call_tree) call calltree_waypoint("state remapped (ALE_inputs)")

610

611 ! Reintegrate mass transports from Zstar to the offline vertical coordinate

612 do j=jsc,jec ; do i=g%iscB,g%iecB

613 if (g%mask2dCu(i,j)>0.) then

614 h_src(:) = 0.5 * (h(i,j,:) + h(i+1,j,:))

615 h_dest(:) = 0.5 * (h_new(i,j,:) + h_new(i+1,j,:))

616 call reintegrate_column(nk, h_src, uhtr(i,j,:), nk, h_dest, temp_vec)

617 uhtr(i,j,:) = temp_vec

618 endif

619 enddo ; enddo

620 do j=g%jscB,g%jecB ; do i=isc,iec

621 if (g%mask2dCv(i,j)>0.) then

622 h_src(:) = 0.5 * (h(i,j,:) + h(i,j+1,:))

623 h_dest(:) = 0.5 * (h_new(i,j,:) + h_new(i,j+1,:))

624 call reintegrate_column(nk, h_src, vhtr(i,j,:), nk, h_dest, temp_vec)

625 vhtr(i,j,:) = temp_vec

626 endif

627 enddo ; enddo

628

629 do j=jsc,jec ; do i=isc,iec

630 if (g%mask2dT(i,j)>0.) then

631 if (check_column_integrals(nk, h_src, nk, h_dest)) then

632 call mom_error(fatal, "ALE_offline_inputs: Kd interpolation columns do not match")

633 endif

634 call interpolate_column(nk, h(i,j,:), kd(i,j,:), nk, h_new(i,j,:), kd(i,j,:), .true.)

635 endif

636 enddo ; enddo

637

638 call ale_remap_scalar(cs%remapCS, g, gv, nk, h, tv%T, h_new, tv%T)

639 call ale_remap_scalar(cs%remapCS, g, gv, nk, h, tv%S, h_new, tv%S)

640

641 if (debug) call mom_tracer_chkinv("After ALE_offline_inputs", g, gv, h_new, reg%Tr, reg%ntr)

642

643 ! Copy over the new layer thicknesses

644 do k = 1,nk ; do j = jsc-1,jec+1 ; do i = isc-1,iec+1

645 h(i,j,k) = h_new(i,j,k)

646 enddo ; enddo ; enddo

647

648 if (allocated(tv%SpV_avg)) tv%valid_SpV_halo = -1 ! Record that SpV_avg is no longer valid.

649

650 if (cs%show_call_tree) call calltree_leave("ALE_offline_inputs()")

651end subroutine ale_offline_inputs

652

653

654!> For a state-based coordinate, accelerate the process of regridding by

655!! repeatedly applying the grid calculation algorithm

656subroutine ale_regrid_accelerated(CS, G, GV, US, h, tv, n_itt, u, v, OBC, Reg, dt, &

657 dzRegrid, initial)

658 type(ale_cs), pointer :: cs !< ALE control structure

659 type(ocean_grid_type), intent(inout) :: g !< Ocean grid

660 type(verticalgrid_type), intent(in) :: gv !< Vertical grid

661 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type

662 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &

663 intent(inout) :: h !< Original thicknesses [H ~> m or kg m-2]

664 type(thermo_var_ptrs), intent(inout) :: tv !< Thermo vars (T/S/EOS)

665 integer, intent(in) :: n_itt !< Number of times to regrid

666 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &

667 intent(inout) :: u !< Zonal velocity [L T-1 ~> m s-1]

668 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &

669 intent(inout) :: v !< Meridional velocity [L T-1 ~> m s-1]

670 type(ocean_obc_type), pointer :: obc !< Open boundary structure

671 type(tracer_registry_type), &

672 optional, pointer :: reg !< Tracer registry to remap onto new grid

673 real, optional, intent(in) :: dt !< Model timestep to provide a timescale for regridding [T ~> s]

674 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &

675 optional, intent(inout) :: dzregrid !< Final change in interface positions [H ~> m or kg m-2]

676 logical, optional, intent(in) :: initial !< Whether we're being called from an initialization

677 !! routine (and expect diagnostics to work)

678

679 ! Local variables

680 integer :: i, j, itt, nz

681 type(thermo_var_ptrs) :: tv_local ! local/intermediate temp/salt

682 type(group_pass_type) :: pass_t_s_h ! group pass if the coordinate has a stencil

683 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h_loc ! A working copy of layer thicknesses [H ~> m or kg m-2]

684 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h_orig ! The original layer thicknesses [H ~> m or kg m-2]

685 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), target :: t ! local temporary temperatures [C ~> degC]

686 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), target :: s ! local temporary salinities [S ~> ppt]

687 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: h_old_u ! Source grid thickness at zonal

688 ! velocity points [H ~> m or kg m-2]

689 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: h_old_v ! Source grid thickness at meridional

690 ! velocity points [H ~> m or kg m-2]

691 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: h_new_u ! Destination grid thickness at zonal

692 ! velocity points [H ~> m or kg m-2]

693 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: h_new_v ! Destination grid thickness at meridional

694 ! velocity points [H ~> m or kg m-2]

695

696 ! we have to keep track of the total dzInterface if for some reason

697 ! we're using the old remapping algorithm for u/v

698 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: dzinterface ! Interface height changes within

699 ! an iteration [H ~> m or kg m-2]

700 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: dzinttotal ! Cumulative interface position changes [H ~> m or kg m-2]

701

702 nz = gv%ke

703

704 ! initial total interface displacement due to successive regridding

705 if (cs%remap_uv_using_old_alg) &

706 dzinttotal(:,:,:) = 0.

707

708 call create_group_pass(pass_t_s_h, t, g%domain)

709 call create_group_pass(pass_t_s_h, s, g%domain)

710 call create_group_pass(pass_t_s_h, h_loc, g%domain)

711

712 ! copy original temp/salt and set our local tv_pointers to them

713 tv_local = tv

714 t(:,:,:) = tv%T(:,:,:)

715 s(:,:,:) = tv%S(:,:,:)

716 tv_local%T => t

717 tv_local%S => s

718

719 ! get local copy of thickness and save original state for remapping

720 h_loc(:,:,:) = h(:,:,:)

721 h_orig(:,:,:) = h(:,:,:)

722

723 ! Apply timescale to regridding (for e.g. filtered_grid_motion)

724 if (present(dt)) &

725 call ale_update_regrid_weights(dt, cs)

726

727 do itt = 1, n_itt

728

729 call do_group_pass(pass_t_s_h, g%domain)

730

731 ! generate new grid

732 if (cs%do_conv_adj) call convective_adjustment(g, gv, h_loc, tv_local)

733

734 ! Update the layer specific volumes if necessary

735 if (allocated(tv_local%SpV_avg)) call calc_derived_thermo(tv_local, h, g, gv, us, halo=1)

736

737 call regridding_main(cs%remapCS, cs%regridCS, g, gv, us, h_loc, tv_local, h, dzinterface)

738 if (cs%remap_uv_using_old_alg) &

739 dzinttotal(:,:,:) = dzinttotal(:,:,:) + dzinterface(:,:,:)

740

741 ! remap from original grid onto new grid

742 do j = g%jsc-1,g%jec+1 ; do i = g%isc-1,g%iec+1

743 call remapping_core_h(cs%remapCS, nz, h_orig(i,j,:), tv%S(i,j,:), nz, h(i,j,:), &

744 tv_local%S(i,j,:))

745 call remapping_core_h(cs%remapCS, nz, h_orig(i,j,:), tv%T(i,j,:), nz, h(i,j,:), &

746 tv_local%T(i,j,:))

747 enddo ; enddo

748

749 ! starting grid for next iteration

750 h_loc(:,:,:) = h(:,:,:)

751 enddo

752

753 ! remap all state variables (including those that weren't needed for regridding)

754 call ale_remap_tracers(cs, g, gv, h_orig, h, reg)

755

756 call ale_remap_set_h_vel(cs, g, gv, h_orig, h_old_u, h_old_v, obc)

757 if (cs%remap_uv_using_old_alg) then

758 call ale_remap_set_h_vel_via_dz(cs, g, gv, h, h_new_u, h_new_v, obc, h_orig, dzinttotal)

759 else

760 call ale_remap_set_h_vel(cs, g, gv, h, h_new_u, h_new_v, obc)

761 endif

762

763 call ale_remap_velocities(cs, g, gv, h_old_u, h_old_v, h_new_u, h_new_v, u, v)

764

765 ! save total dzregrid for diags if needed?

766 if (present(dzregrid)) dzregrid(:,:,:) = dzinttotal(:,:,:)

767

768 if (allocated(tv%SpV_avg)) tv%valid_SpV_halo = -1 ! Record that SpV_avg is no longer valid.

769

770end subroutine ale_regrid_accelerated

771

772!> This routine takes care of remapping all tracer variables between the old and the

773!! new grids. This routine is called during initialization of the model at time=0, to

774!! remap initial conditions to the model grid. It is also called during a

775!! time step to update the state.

776subroutine ale_remap_tracers(CS, G, GV, h_old, h_new, Reg, debug, dt, PCM_cell)

777 type(ale_cs), intent(in) :: cs !< ALE control structure

778 type(ocean_grid_type), intent(in) :: g !< Ocean grid structure

779 type(verticalgrid_type), intent(in) :: gv !< Ocean vertical grid structure

780 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_old !< Thickness of source grid

781 !! [H ~> m or kg m-2]

782 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_new !< Thickness of destination grid

783 !! [H ~> m or kg m-2]

784 type(tracer_registry_type), pointer :: reg !< Tracer registry structure

785 logical, optional, intent(in) :: debug !< If true, show the call tree

786 real, optional, intent(in) :: dt !< time step for diagnostics [T ~> s]

787 logical, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &

788 optional, intent(in) :: pcm_cell !< Use PCM remapping in cells where true

789

790 ! Local variables

791 real :: tr_column(gv%ke) ! A column of updated tracer concentrations [CU ~> Conc]

792 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: work_conc ! The rate of change of concentrations [Conc T-1 ~> Conc s-1]

793 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: work_cont ! The rate of change of cell-integrated tracer

794 ! content [Conc H T-1 ~> Conc m s-1 or Conc kg m-2 s-1] or

795 ! cell thickness [H T-1 ~> m s-1 or kg m-2 s-1]

796 real, dimension(SZI_(G),SZJ_(G)) :: work_2d ! The rate of change of column-integrated tracer

797 ! content [Conc H T-1 ~> Conc m s-1 or Conc kg m-2 s-1]

798 logical :: pcm(gv%ke) ! If true, do PCM remapping from a cell.

799 real :: idt ! The inverse of the timestep [T-1 ~> s-1]

800 real :: h1(gv%ke) ! A column of source grid layer thicknesses [H ~> m or kg m-2]

801 real :: h2(gv%ke) ! A column of target grid layer thicknesses [H ~> m or kg m-2]

802 logical :: show_call_tree

803 type(tracer_type), pointer :: tr => null()

804 integer :: i, j, k, m, nz, ntr

805

806 show_call_tree = .false.

807 if (present(debug)) show_call_tree = debug

808

809 if (show_call_tree) call calltree_enter("ALE_remap_tracers(), MOM_ALE.F90")

810

811 nz = gv%ke

812

813 ntr = 0 ; if (associated(reg)) ntr = reg%ntr

814

815 if (present(dt)) then

816 idt = 1.0/dt

817 work_conc(:,:,:) = 0.0

818 work_cont(:,:,:) = 0.0

819 endif

820

821 ! Remap all registered tracers, including temperature and salinity.

822 if (ntr>0) then

823 if (show_call_tree) call calltree_waypoint("remapping tracers (ALE_remap_tracers)")

824 !$OMP parallel do default(shared) private(h1,h2,tr_column,Tr,PCM,work_conc,work_cont,work_2d)

825 do m=1,ntr ! For each tracer

826 tr => reg%Tr(m)

827 do j = g%jsc,g%jec ; do i = g%isc,g%iec ; if (g%mask2dT(i,j)>0.) then

828 ! Build the start and final grids

829 h1(:) = h_old(i,j,:)

830 h2(:) = h_new(i,j,:)

831 if (present(pcm_cell)) then

832 pcm(:) = pcm_cell(i,j,:)

833 call remapping_core_h(cs%remapCS, nz, h1, tr%t(i,j,:), nz, h2, tr_column, pcm_cell=pcm)

834 else

835 call remapping_core_h(cs%remapCS, nz, h1, tr%t(i,j,:), nz, h2, tr_column)

836 endif

837

838 ! Possibly underflow any very tiny tracer concentrations to 0. Note that this is not conservative!

839 if (tr%conc_underflow > 0.0) then ; do k=1,gv%ke

840 if (abs(tr_column(k)) < tr%conc_underflow) tr_column(k) = 0.0

841 enddo ; endif

842

843 ! Intermediate steps for tendency of tracer concentration and tracer content.

844 if (present(dt)) then

845 if (tr%id_remap_conc > 0) then

846 do k=1,gv%ke

847 work_conc(i,j,k) = (tr_column(k) - tr%t(i,j,k)) * idt

848 enddo

849 endif

850 if (tr%id_remap_cont > 0 .or. tr%id_remap_cont_2d > 0) then

851 do k=1,gv%ke

852 work_cont(i,j,k) = (tr_column(k)*h2(k) - tr%t(i,j,k)*h1(k)) * idt

853 enddo

854 endif

855 endif

856

857 ! update tracer concentration

858 tr%t(i,j,:) = tr_column(:)

859 endif ; enddo ; enddo

860

861 ! tendency diagnostics.

862 if (present(dt)) then

863 if (tr%id_remap_conc > 0) then

864 call post_data(tr%id_remap_conc, work_conc, cs%diag)

865 endif

866 if (tr%id_remap_cont > 0) then

867 call post_data(tr%id_remap_cont, work_cont, cs%diag)

868 endif

869

870 if (tr%id_remap_cont_2d > 0) then

871 do j = g%jsc,g%jec ; do i = g%isc,g%iec

872 work_2d(i,j) = 0.0

873 do k = 1,gv%ke

874 work_2d(i,j) = work_2d(i,j) + work_cont(i,j,k)

875 enddo

876 enddo ; enddo

877 call post_data(tr%id_remap_cont_2d, work_2d, cs%diag)

878 endif

879 endif

880 enddo ! m=1,ntr

881

882 endif ! endif for ntr > 0

883

884

885 if (cs%id_vert_remap_h > 0) call post_data(cs%id_vert_remap_h, h_old, cs%diag)

886 if ((cs%id_vert_remap_h_tendency > 0) .and. present(dt)) then

887 do k = 1, nz ; do j = g%jsc,g%jec ; do i = g%isc,g%iec

888 work_cont(i,j,k) = (h_new(i,j,k) - h_old(i,j,k))*idt

889 enddo ; enddo ; enddo

890 call post_data(cs%id_vert_remap_h_tendency, work_cont, cs%diag)

891 endif

892

893 if (show_call_tree) call calltree_leave("ALE_remap_tracers(), MOM_ALE.F90")

894

895end subroutine ale_remap_tracers

896

897!> This routine sets the thicknesses at velocity points used for vertical remapping.

898subroutine ale_remap_set_h_vel(CS, G, GV, h_new, h_u, h_v, OBC, debug)

899 type(ale_cs), intent(in) :: cs !< ALE control structure

900 type(ocean_grid_type), intent(in) :: g !< Ocean grid structure

901 type(verticalgrid_type), intent(in) :: gv !< Ocean vertical grid structure

902 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_new !< Thickness at tracer points of the

903 !! grid being interpolated to velocity

904 !! points [H ~> m or kg m-2]

905 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &

906 intent(inout) :: h_u !< Grid thickness at zonal velocity

907 !! points [H ~> m or kg m-2]

908 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &

909 intent(inout) :: h_v !< Grid thickness at meridional velocity

910 !! points [H ~> m or kg m-2]

911 type(ocean_obc_type), pointer :: obc !< Open boundary structure

912 logical, optional, intent(in) :: debug !< If true, show the call tree

913

914 ! Local variables

915 logical :: show_call_tree

916 integer :: i, j, k

917

918 show_call_tree = .false.

919 if (present(debug)) show_call_tree = debug

920 if (show_call_tree) call calltree_enter("ALE_remap_set_h_vel()")

921

922 ! Build the u- and v-velocity grid thicknesses for remapping.

923

924 !$OMP parallel do default(shared)

925 do k=1,gv%ke ; do j=g%jsc,g%jec ; do i=g%IscB,g%IecB ; if (g%mask2dCu(i,j)>0.) then

926 h_u(i,j,k) = 0.5*(h_new(i,j,k) + h_new(i+1,j,k))

927 endif ; enddo ; enddo ; enddo

928 !$OMP parallel do default(shared)

929 do k=1,gv%ke ; do j=g%JscB,g%JecB ; do i=g%isc,g%iec ; if (g%mask2dCv(i,j)>0.) then

930 h_v(i,j,k) = 0.5*(h_new(i,j,k) + h_new(i,j+1,k))

931 endif ; enddo ; enddo ; enddo

932

933 ! Mask out blocked portions of velocity cells.

934 if (cs%partial_cell_vel_remap) call ale_remap_set_h_vel_partial(cs, g, gv, h_new, h_u, h_v)

935

936 ! Take open boundary conditions into account.

937 if (associated(obc)) call ale_remap_set_h_vel_obc(g, gv, h_new, h_u, h_v, obc)

938

939 if (show_call_tree) call calltree_leave("ALE_remap_set_h_vel()")

940

941end subroutine ale_remap_set_h_vel

942

943!> This routine sets the thicknesses at velocity points used for vertical remapping using a

944!! combination of the old grid and interface movements.

945subroutine ale_remap_set_h_vel_via_dz(CS, G, GV, h_new, h_u, h_v, OBC, h_old, dzInterface, debug)

946 type(ale_cs), intent(in) :: cs !< ALE control structure

947 type(ocean_grid_type), intent(in) :: g !< Ocean grid structure

948 type(verticalgrid_type), intent(in) :: gv !< Ocean vertical grid structure

949 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_new !< Thickness at tracer points of the

950 !! grid being interpolated to velocity

951 !! points [H ~> m or kg m-2]

952 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &

953 intent(inout) :: h_u !< Grid thickness at zonal velocity

954 !! points [H ~> m or kg m-2]

955 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &

956 intent(inout) :: h_v !< Grid thickness at meridional velocity

957 !! points [H ~> m or kg m-2]

958 type(ocean_obc_type), pointer :: obc !< Open boundary structure

959 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &

960 intent(in) :: h_old !< Thickness of source grid when generating

961 !! the destination grid via the old

962 !! algorithm [H ~> m or kg m-2]

963 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &

964 intent(in) :: dzinterface !< Change in interface position

965 !! [H ~> m or kg m-2]

966 logical, optional, intent(in) :: debug !< If true, show the call tree

967

968 ! Local variables

969 logical :: show_call_tree

970 integer :: i, j, k

971

972 show_call_tree = .false.

973 if (present(debug)) show_call_tree = debug

974 if (show_call_tree) call calltree_enter("ALE_remap_set_h_vel()")

975

976 ! Build the u- and v-velocity grid thicknesses for remapping using the old grid and interface movement.

977

978 !$OMP parallel do default(shared)

979 do k=1,gv%ke ; do j=g%jsc,g%jec ; do i=g%IscB,g%IecB ; if (g%mask2dCu(i,j)>0.) then

980 h_u(i,j,k) = max( 0., 0.5*(h_old(i,j,k) + h_old(i+1,j,k)) + &

981 0.5 * (( dzinterface(i,j,k) + dzinterface(i+1,j,k) ) - &

982 ( dzinterface(i,j,k+1) + dzinterface(i+1,j,k+1) )) )

983 endif ; enddo ; enddo ; enddo

984

985 !$OMP parallel do default(shared)

986 do k=1,gv%ke ; do j=g%JscB,g%JecB ; do i=g%isc,g%iec ; if (g%mask2dCv(i,j)>0.) then

987 h_v(i,j,k) = max( 0., 0.5*(h_old(i,j,k) + h_old(i,j+1,k)) + &

988 0.5 * (( dzinterface(i,j,k) + dzinterface(i,j+1,k) ) - &

989 ( dzinterface(i,j,k+1) + dzinterface(i,j+1,k+1) )) )

990 endif ; enddo ; enddo ; enddo

991

992 ! Mask out blocked portions of velocity cells.

993 if (cs%partial_cell_vel_remap) call ale_remap_set_h_vel_partial(cs, g, gv, h_old, h_u, h_v)

994

995 ! Take open boundary conditions into account.

996 if (associated(obc)) call ale_remap_set_h_vel_obc(g, gv, h_new, h_u, h_v, obc)

997

998 if (show_call_tree) call calltree_leave("ALE_remap_set_h_vel()")

999

1000end subroutine ale_remap_set_h_vel_via_dz

1001

1002!> Mask out the thicknesses at velocity points where they are below the minimum depth

1003!! at adjacent tracer points

1004subroutine ale_remap_set_h_vel_partial(CS, G, GV, h_mask, h_u, h_v)

1005 type(ale_cs), intent(in) :: CS !< ALE control structure

1006 type(ocean_grid_type), intent(in) :: G !< Ocean grid structure

1007 type(verticalgrid_type), intent(in) :: GV !< Ocean vertical grid structure

1008 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_mask !< Thickness at tracer points

1009 !! used to apply the partial

1010 !! cell masking [H ~> m or kg m-2]

1011 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &

1012 intent(inout) :: h_u !< Grid thickness at zonal velocity

1013 !! points [H ~> m or kg m-2]

1014 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &

1015 intent(inout) :: h_v !< Grid thickness at meridional velocity

1016 !! points [H ~> m or kg m-2]

1017 ! Local variables

1018 real, dimension(SZI_(G),SZJ_(G)) :: h_tot ! The vertically summed thicknesses [H ~> m or kg m-2]

1019 real :: h_mask_vel ! A depth below which the thicknesses at a velocity point are masked out [H ~> m or kg m-2]

1020 integer :: i, j, k

1021

1022 h_tot(:,:) = 0.0

1023 do k=1,gv%ke ; do j=g%jsc-1,g%jec+1 ; do i=g%isc-1,g%iec+1

1024 h_tot(i,j) = h_tot(i,j) + h_mask(i,j,k)

1025 enddo ; enddo ; enddo

1026

1027 !$OMP parallel do default(shared) private(h_mask_vel)

1028 do j=g%jsc,g%jec ; do i=g%IscB,g%IecB ; if (g%mask2dCu(i,j)>0.) then

1029 h_mask_vel = min(h_tot(i,j), h_tot(i+1,j))

1030 call apply_partial_cell_mask(h_u(i,j,:), h_mask_vel)

1031 endif ; enddo ; enddo

1032

1033 !$OMP parallel do default(shared) private(h_mask_vel)

1034 do j=g%JscB,g%JecB ; do i=g%isc,g%iec ; if (g%mask2dCv(i,j)>0.) then

1035 h_mask_vel = min(h_tot(i,j), h_tot(i,j+1))

1036 call apply_partial_cell_mask(h_v(i,j,:), h_mask_vel)

1037 endif ; enddo ; enddo

1038

1039end subroutine ale_remap_set_h_vel_partial

1040

1041! Reset thicknesses at velocity points on open boundary condition segments

1042subroutine ale_remap_set_h_vel_obc(G, GV, h_new, h_u, h_v, OBC)

1043 type(ocean_grid_type), intent(in) :: G !< Ocean grid structure

1044 type(verticalgrid_type), intent(in) :: GV !< Ocean vertical grid structure

1045 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_new !< Thickness at tracer points of the

1046 !! grid being interpolated to velocity

1047 !! points [H ~> m or kg m-2]

1048 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &

1049 intent(inout) :: h_u !< Grid thickness at zonal velocity

1050 !! points [H ~> m or kg m-2]

1051 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &

1052 intent(inout) :: h_v !< Grid thickness at meridional velocity

1053 !! points [H ~> m or kg m-2]

1054 type(ocean_obc_type), pointer :: OBC !< Open boundary structure

1055

1056 ! Local variables

1057 integer :: i, j, k, nz, is_OBC, ie_OBC, js_OBC, je_OBC

1058

1059 if (.not.associated(obc)) return

1060

1061 nz = gv%ke

1062

1063 ! Take open boundary conditions into account.

1064 if (obc%u_E_OBCs_on_PE) then

1065 js_obc = max(g%jsc, obc%js_u_E_obc) ; je_obc = min(g%jec, obc%je_u_E_obc)

1066 is_obc = max(g%IscB, obc%Is_u_E_obc) ; ie_obc = min(g%IecB, obc%Ie_u_E_obc)

1067 !$OMP parallel do default(shared)

1068 do j=js_obc,je_obc ; do i=is_obc,ie_obc ; if (obc%segnum_u(i,j) > 0) then ! OBC_DIRECTION_E

1069 do k=1,nz ; h_u(i,j,k) = h_new(i,j,k) ; enddo

1070 endif ; enddo ; enddo

1071 endif

1072 if (obc%u_W_OBCs_on_PE) then

1073 js_obc = max(g%jsc, obc%js_u_W_obc) ; je_obc = min(g%jec, obc%je_u_W_obc)

1074 is_obc = max(g%IscB, obc%Is_u_W_obc) ; ie_obc = min(g%IecB, obc%Ie_u_W_obc)

1075 !$OMP parallel do default(shared)

1076 do j=js_obc,je_obc ; do i=is_obc,ie_obc ; if (obc%segnum_u(i,j) < 0) then ! OBC_DIRECTION_W

1077 do k=1,nz ; h_u(i,j,k) = h_new(i+1,j,k) ; enddo

1078 endif ; enddo ; enddo

1079 endif

1080

1081 if (obc%v_N_OBCs_on_PE) then

1082 js_obc = max(g%JscB, obc%Js_v_N_obc) ; je_obc = min(g%JecB, obc%Je_v_N_obc)

1083 is_obc = max(g%isc, obc%is_v_N_obc) ; ie_obc = min(g%iec, obc%ie_v_N_obc)

1084 !$OMP parallel do default(shared)

1085 do j=js_obc,je_obc ; do i=is_obc,ie_obc ; if (obc%segnum_v(i,j) > 0) then ! OBC_DIRECTION_N

1086 do k=1,nz ; h_v(i,j,k) = h_new(i,j,k) ; enddo

1087 endif ; enddo ; enddo

1088 endif

1089 if (obc%v_S_OBCs_on_PE) then

1090 js_obc = max(g%JscB, obc%Js_v_S_obc) ; je_obc = min(g%JecB, obc%Je_v_S_obc)

1091 is_obc = max(g%isc, obc%is_v_S_obc) ; ie_obc = min(g%iec, obc%ie_v_S_obc)

1092 !$OMP parallel do default(shared)

1093 do j=js_obc,je_obc ; do i=is_obc,ie_obc ; if (obc%segnum_v(i,j) < 0) then ! OBC_DIRECTION_S

1094 do k=1,nz ; h_v(i,j,k) = h_new(i,j+1,k) ; enddo

1095 endif ; enddo ; enddo

1096 endif

1097

1098end subroutine ale_remap_set_h_vel_obc

1099

1100!> This routine remaps velocity components between the old and the new grids,

1101!! with thicknesses at velocity points taken to be arithmetic averages of tracer thicknesses.

1102!! This routine may be called during initialization of the model at time=0, to

1103!! remap initial conditions to the model grid. It is also called during a

1104!! time step to update the state.

1105subroutine ale_remap_velocities(CS, G, GV, h_old_u, h_old_v, h_new_u, h_new_v, u, v, debug, &

1106 dt, allow_preserve_variance)

1107 type(ale_cs), intent(in) :: cs !< ALE control structure

1108 type(ocean_grid_type), intent(in) :: g !< Ocean grid structure

1109 type(verticalgrid_type), intent(in) :: gv !< Ocean vertical grid structure

1110 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &

1111 intent(in) :: h_old_u !< Source grid thickness at zonal

1112 !! velocity points [H ~> m or kg m-2]

1113 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &

1114 intent(in) :: h_old_v !< Source grid thickness at meridional

1115 !! velocity points [H ~> m or kg m-2]

1116 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &

1117 intent(in) :: h_new_u !< Destination grid thickness at zonal

1118 !! velocity points [H ~> m or kg m-2]

1119 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &

1120 intent(in) :: h_new_v !< Destination grid thickness at meridional

1121 !! velocity points [H ~> m or kg m-2]

1122 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &

1123 intent(inout) :: u !< Zonal velocity [L T-1 ~> m s-1]

1124 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &

1125 intent(inout) :: v !< Meridional velocity [L T-1 ~> m s-1]

1126 logical, optional, intent(in) :: debug !< If true, show the call tree

1127 real, optional, intent(in) :: dt !< time step for diagnostics [T ~> s]

1128 logical, optional, intent(in) :: allow_preserve_variance !< If true, enables ke-conserving

1129 !! correction

1130

1131 ! Local variables

1132 real :: h_mask_vel ! A depth below which the thicknesses at a velocity point are masked out [H ~> m or kg m-2]

1133 real :: u_src(gv%ke) ! A column of u-velocities on the source grid [L T-1 ~> m s-1]

1134 real :: u_tgt(gv%ke) ! A column of u-velocities on the target grid [L T-1 ~> m s-1]

1135 real :: v_src(gv%ke) ! A column of v-velocities on the source grid [L T-1 ~> m s-1]

1136 real :: v_tgt(gv%ke) ! A column of v-velocities on the target grid [L T-1 ~> m s-1]

1137 real :: h1(gv%ke) ! A column of source grid layer thicknesses [H ~> m or kg m-2]

1138 real :: h2(gv%ke) ! A column of target grid layer thicknesses [H ~> m or kg m-2]

1139 real :: rescale_coef ! Factor that scales the baroclinic velocity to conserve ke [nondim]

1140 real :: u_bt, v_bt ! Depth-averaged velocity components [L T-1 ~> m s-1]

1141 real :: ke_c_src, ke_c_tgt ! \int [u_c or v_c]^2 dz on src and tgt grids [H L2 T-2 ~> m3 s-2]

1142 real, dimension(SZIB_(G),SZJ_(G)) :: du2h_tot ! The rate of change of vertically integrated

1143 ! 0.5 * rho0 * u**2 [R Z L2 T-3 ~> W m-2]

1144 real, dimension(SZI_(G),SZJB_(G)) :: dv2h_tot ! The rate of change of vertically integrated

1145 ! 0.5 * rho0 * v**2 [R Z L2 T-3 ~> W m-2]

1146 real :: u2h_tot, v2h_tot ! The vertically integrated u**2 and v**2 [H L2 T-2 ~> m3 s-2 or kg s-2]

1147 real :: i_dt ! 1 / dt [T-1 ~> s-1]

1148 logical :: variance_option ! Contains the value of allow_preserve_variance when present, else false

1149 logical :: show_call_tree

1150 integer :: i, j, k, nz

1151

1152 show_call_tree = .false.

1153 if (present(debug)) show_call_tree = debug

1154 if (show_call_tree) call calltree_enter("ALE_remap_velocities()")

1155

1156 ! Setup related to KE conservation

1157 variance_option = .false.

1158 if (present(allow_preserve_variance)) variance_option=allow_preserve_variance

1159 if (present(dt)) i_dt = 1.0 / dt

1160

1161 if (cs%id_remap_delta_integ_u2>0) du2h_tot(:,:) = 0.

1162 if (cs%id_remap_delta_integ_v2>0) dv2h_tot(:,:) = 0.

1163

1164 if (((cs%id_remap_delta_integ_u2>0) .or. (cs%id_remap_delta_integ_v2>0)) .and. .not.present(dt))&

1165 call mom_error(fatal, "ALE KE diagnostics requires passing dt into ALE_remap_velocities")

1166

1167 nz = gv%ke

1168

1169 ! --- Remap u profiles from the source vertical grid onto the new target grid.

1170

1171 !$OMP parallel do default(shared) private(h1,h2,u_src,h_mask_vel,u_tgt, &

1172 !$OMP u_bt,ke_c_src,ke_c_tgt,rescale_coef, &

1173 !$OMP u2h_tot,v2h_tot)

1174 do j=g%jsc,g%jec ; do i=g%IscB,g%IecB ; if (g%mask2dCu(i,j)>0.) then

1175 ! Make a 1-d copy of the start and final grids and the source velocity

1176 do k=1,nz

1177 h1(k) = h_old_u(i,j,k)

1178 h2(k) = h_new_u(i,j,k)

1179 u_src(k) = u(i,j,k)

1180 enddo

1181

1182 if (cs%id_remap_delta_integ_u2>0) then

1183 u2h_tot = 0.

1184 do k=1,nz

1185 u2h_tot = u2h_tot - h1(k) * (u_src(k)**2)

1186 enddo

1187 endif

1188

1189 call remapping_core_h(cs%vel_remapCS, nz, h1, u_src, nz, h2, u_tgt)

1190

1191 if (variance_option .and. cs%conserve_ke) then

1192 ! Conserve ke_u by correcting baroclinic component.

1193 ! Assumes total depth doesn't change during remap, and

1194 ! that \int u(z) dz doesn't change during remap.

1195 ! First get barotropic component

1196 u_bt = 0.0

1197 do k=1,nz

1198 u_bt = u_bt + h2(k) * u_tgt(k) ! Dimensions [H L T-1 ~> m2 s-1 or kg m-1 s-1]

1199 enddo

1200 u_bt = u_bt / (sum(h2(1:nz)) + gv%H_subroundoff) ! Dimensions return to [L T-1 ~> m s-1]

1201 ! Next get baroclinic ke = \int (u-u_bt)^2 from source and target

1202 ke_c_src = 0.0

1203 ke_c_tgt = 0.0

1204 do k=1,nz

1205 ke_c_src = ke_c_src + h1(k) * (u_src(k) - u_bt)**2

1206 ke_c_tgt = ke_c_tgt + h2(k) * (u_tgt(k) - u_bt)**2

1207 enddo

1208 ! Next rescale baroclinic component on target grid to conserve ke

1209 ! The values 1.5625 = 1.25**2 and 1.25 below mean that the KE-conserving

1210 ! correction cannot amplify the baroclinic part of velocity by more

1211 ! than 25%. This threshold is somewhat arbitrary. It was added to

1212 ! prevent unstable behavior when the amplification factor is large.

1213 if (ke_c_src < 1.5625 * ke_c_tgt) then

1214 rescale_coef = sqrt(ke_c_src / ke_c_tgt)

1215 else

1216 rescale_coef = 1.25

1217 endif

1218 do k=1,nz

1219 u_tgt(k) = u_bt + rescale_coef * (u_tgt(k) - u_bt)

1220 enddo

1221 endif

1222

1223 if (cs%id_remap_delta_integ_u2>0) then

1224 do k=1,nz

1225 u2h_tot = u2h_tot + h2(k) * (u_tgt(k)**2)

1226 enddo

1227 du2h_tot(i,j) = u2h_tot * i_dt

1228 endif

1229

1230 if ((cs%BBL_h_vel_mask > 0.0) .and. (cs%h_vel_mask > 0.0)) &

1231 call mask_near_bottom_vel(u_tgt, h2, cs%BBL_h_vel_mask, cs%h_vel_mask, nz)

1232

1233 ! Copy the column of new velocities back to the 3-d array

1234 do k=1,nz

1235 u(i,j,k) = u_tgt(k)

1236 enddo !k

1237 endif ; enddo ; enddo

1238

1239 if (cs%id_remap_delta_integ_u2>0) call post_data(cs%id_remap_delta_integ_u2, du2h_tot, cs%diag)

1240

1241 if (show_call_tree) call calltree_waypoint("u remapped (ALE_remap_velocities)")

1242

1243

1244 ! --- Remap v profiles from the source vertical grid onto the new target grid.

1245

1246 !$OMP parallel do default(shared) private(h1,h2,v_src,h_mask_vel,v_tgt, &

1247 !$OMP v_bt,ke_c_src,ke_c_tgt,rescale_coef, &

1248 !$OMP u2h_tot,v2h_tot)

1249 do j=g%JscB,g%JecB ; do i=g%isc,g%iec ; if (g%mask2dCv(i,j)>0.) then

1250

1251 do k=1,nz

1252 h1(k) = h_old_v(i,j,k)

1253 h2(k) = h_new_v(i,j,k)

1254 v_src(k) = v(i,j,k)

1255 enddo

1256

1257 if (cs%id_remap_delta_integ_v2>0) then

1258 v2h_tot = 0.

1259 do k=1,nz

1260 v2h_tot = v2h_tot - h1(k) * (v_src(k)**2)

1261 enddo

1262 endif

1263

1264 call remapping_core_h(cs%vel_remapCS, nz, h1, v_src, nz, h2, v_tgt)

1265

1266 if (variance_option .and. cs%conserve_ke) then

1267 ! Conserve ke_v by correcting baroclinic component.

1268 ! Assumes total depth doesn't change during remap, and

1269 ! that \int v(z) dz doesn't change during remap.

1270 ! First get barotropic component

1271 v_bt = 0.0

1272 do k=1,nz

1273 v_bt = v_bt + h2(k) * v_tgt(k) ! Dimensions [H L T-1 ~> m2 s-1 or kg m-1 s-1]

1274 enddo

1275 v_bt = v_bt / (sum(h2(1:nz)) + gv%H_subroundoff) ! Dimensions return to [L T-1 ~> m s-1]

1276 ! Next get baroclinic ke = \int (u-u_bt)^2 from source and target

1277 ke_c_src = 0.0

1278 ke_c_tgt = 0.0

1279 do k=1,nz

1280 ke_c_src = ke_c_src + h1(k) * (v_src(k) - v_bt)**2

1281 ke_c_tgt = ke_c_tgt + h2(k) * (v_tgt(k) - v_bt)**2

1282 enddo

1283 ! Next rescale baroclinic component on target grid to conserve ke

1284 if (ke_c_src < 1.5625 * ke_c_tgt) then

1285 rescale_coef = sqrt(ke_c_src / ke_c_tgt)

1286 else

1287 rescale_coef = 1.25

1288 endif

1289 do k=1,nz

1290 v_tgt(k) = v_bt + rescale_coef * (v_tgt(k) - v_bt)

1291 enddo

1292 endif

1293

1294 if (cs%id_remap_delta_integ_v2>0) then

1295 do k=1,nz

1296 v2h_tot = v2h_tot + h2(k) * (v_tgt(k)**2)

1297 enddo

1298 dv2h_tot(i,j) = v2h_tot * i_dt

1299 endif

1300

1301 if ((cs%BBL_h_vel_mask > 0.0) .and. (cs%h_vel_mask > 0.0)) then

1302 call mask_near_bottom_vel(v_tgt, h2, cs%BBL_h_vel_mask, cs%h_vel_mask, nz)

1303 endif

1304

1305 ! Copy the column of new velocities back to the 3-d array

1306 do k=1,nz

1307 v(i,j,k) = v_tgt(k)

1308 enddo !k

1309 endif ; enddo ; enddo

1310

1311 if (cs%id_remap_delta_integ_v2>0) call post_data(cs%id_remap_delta_integ_v2, dv2h_tot, cs%diag)

1312

1313 if (show_call_tree) call calltree_waypoint("v remapped (ALE_remap_velocities)")

1314 if (show_call_tree) call calltree_leave("ALE_remap_velocities()")

1315

1316end subroutine ale_remap_velocities

1317

1318!> Interpolate to find an updated array of values at interfaces after remapping.

1319subroutine ale_remap_interface_vals(CS, G, GV, h_old, h_new, int_val)

1320 type(ale_cs), intent(in) :: cs !< ALE control structure

1321 type(ocean_grid_type), intent(in) :: g !< Ocean grid structure

1322 type(verticalgrid_type), intent(in) :: gv !< Ocean vertical grid structure

1323 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_old !< Thickness of source grid

1324 !! [H ~> m or kg m-2]

1325 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_new !< Thickness of destination grid

1326 !! [H ~> m or kg m-2]

1327 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &

1328 intent(inout) :: int_val !< The interface values to interpolate [A]

1329

1330 real :: val_src(gv%ke+1) ! A column of interface values on the source grid [A]

1331 real :: val_tgt(gv%ke+1) ! A column of interface values on the target grid [A]

1332 real :: h_src(gv%ke) ! A column of source grid layer thicknesses [H ~> m or kg m-2]

1333 real :: h_tgt(gv%ke) ! A column of target grid layer thicknesses [H ~> m or kg m-2]

1334 integer :: i, j, k, nz

1335

1336 nz = gv%ke

1337

1338 do j=g%jsc,g%jec ; do i=g%isc,g%iec ; if (g%mask2dT(i,j)>0.) then

1339 do k=1,nz

1340 h_src(k) = h_old(i,j,k)

1341 h_tgt(k) = h_new(i,j,k)

1342 enddo

1343

1344 do k=1,nz+1

1345 val_src(k) = int_val(i,j,k)

1346 enddo

1347

1348 call interpolate_column(nz, h_src, val_src, nz, h_tgt, val_tgt, .false.)

1349

1350 do k=1,nz+1

1351 int_val(i,j,k) = val_tgt(k)

1352 enddo

1353 endif ; enddo ; enddo

1354

1355end subroutine ale_remap_interface_vals

1356

1357!> Interpolate to find an updated array of values at vertices of tracer cells after remapping.

1358subroutine ale_remap_vertex_vals(CS, G, GV, h_old, h_new, vert_val)

1359 type(ale_cs), intent(in) :: cs !< ALE control structure

1360 type(ocean_grid_type), intent(in) :: g !< Ocean grid structure

1361 type(verticalgrid_type), intent(in) :: gv !< Ocean vertical grid structure

1362 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_old !< Thickness of source grid

1363 !! [H ~> m or kg m-2]

1364 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_new !< Thickness of destination grid

1365 !! [H ~> m or kg m-2]

1366 real, dimension(SZIB_(G),SZJB_(G),SZK_(GV)+1), &

1367 intent(inout) :: vert_val !< The interface values to interpolate [A]

1368

1369 real :: val_src(gv%ke+1) ! A column of interface values on the source grid [A]

1370 real :: val_tgt(gv%ke+1) ! A column of interface values on the target grid [A]

1371 real :: h_src(gv%ke) ! A column of source grid layer thicknesses [H ~> m or kg m-2]

1372 real :: h_tgt(gv%ke) ! A column of target grid layer thicknesses [H ~> m or kg m-2]

1373 real :: i_mask_sum ! The inverse of the tracer point masks surrounding a corner [nondim]

1374 integer :: i, j, k, nz

1375

1376 nz = gv%ke

1377

1378 do j=g%JscB,g%JecB ; do i=g%IscB,g%IecB

1379 if ((g%mask2dT(i,j) + g%mask2dT(i+1,j+1)) + (g%mask2dT(i+1,j) + g%mask2dT(i,j+1)) > 0.0 ) then

1380 i_mask_sum = 1.0 / ((g%mask2dT(i,j) + g%mask2dT(i+1,j+1)) + &

1381 (g%mask2dT(i+1,j) + g%mask2dT(i,j+1)))

1382

1383 do k=1,nz

1384 h_src(k) = ((g%mask2dT(i,j) * h_old(i,j,k) + g%mask2dT(i+1,j+1) * h_old(i+1,j+1,k)) + &

1385 (g%mask2dT(i+1,j) * h_old(i+1,j,k) + g%mask2dT(i,j+1) * h_old(i,j+1,k)) ) * i_mask_sum

1386 h_tgt(k) = ((g%mask2dT(i,j) * h_new(i,j,k) + g%mask2dT(i+1,j+1) * h_new(i+1,j+1,k)) + &

1387 (g%mask2dT(i+1,j) * h_new(i+1,j,k) + g%mask2dT(i,j+1) * h_new(i,j+1,k)) ) * i_mask_sum

1388 enddo

1389

1390 do k=1,nz+1

1391 val_src(k) = vert_val(i,j,k)

1392 enddo

1393

1394 call interpolate_column(nz, h_src, val_src, nz, h_tgt, val_tgt, .false.)

1395

1396 do k=1,nz+1

1397 vert_val(i,j,k) = val_tgt(k)

1398 enddo

1399 endif ; enddo ; enddo

1400

1401end subroutine ale_remap_vertex_vals

1402

1403!> Mask out thicknesses to 0 when their running sum exceeds a specified value.

1404subroutine apply_partial_cell_mask(h1, h_mask)

1405 real, dimension(:), intent(inout) :: h1 !< A column of thicknesses to be masked out after their

1406 !! running vertical sum exceeds h_mask [H ~> m or kg m-2]

1407 real, intent(in) :: h_mask !< The depth after which the thicknesses in h1 are

1408 !! masked out [H ~> m or kg m-2]

1409 ! Local variables

1410 real :: h1_rsum ! The running sum of h1 [H ~> m or kg m-2]

1411 integer :: k

1412

1413 h1_rsum = 0.0

1414 do k=1,size(h1)

1415 if (h1(k) > h_mask - h1_rsum) then

1416 ! This thickness is reduced because it extends below the shallower neighboring bathymetry.

1417 h1(k) = max(h_mask - h1_rsum, 0.0)

1418 h1_rsum = h_mask

1419 else

1420 h1_rsum = h1_rsum + h1(k)

1421 endif

1422 enddo

1423end subroutine apply_partial_cell_mask

1424

1425

1426!> Zero out velocities in a column in very thin layers near the seafloor

1427subroutine mask_near_bottom_vel(vel, h, h_BBL, h_thin, nk)

1428 integer, intent(in) :: nk !< The number of layers in this column

1429 real, intent(inout) :: vel(nk) !< The velocity component being zeroed out [L T-1 ~> m s-1]

1430 real, intent(in) :: h(nk) !< The layer thicknesses at velocity points [H ~> m or kg m-2]

1431 real, intent(in) :: h_BBL !< The thickness of the near-bottom region over which to apply

1432 !! the filtering [H ~> m or kg m-2]

1433 real, intent(in) :: h_thin !< A layer thickness below which the filtering is applied [H ~> m or kg m-2]

1434

1435 ! Local variables

1436 real :: h_from_bot ! The distance between the top of a layer and the seafloor [H ~> m or kg m-2]

1437 integer :: k

1438

1439 if ((h_bbl < 0.0) .or. (h_thin < 0.0)) return

1440

1441 h_from_bot = 0.0

1442 do k=nk,1,-1

1443 h_from_bot = h_from_bot + h(k)

1444 if (h_from_bot > h_bbl) return

1445 ! Set the velocity to zero in thin, near-bottom layers.

1446 if (h(k) <= h_thin) vel(k) = 0.0

1447 enddo !k

1448

1449end subroutine mask_near_bottom_vel

1450

1451

1452!> Remaps a single scalar between grids described by thicknesses h_src and h_dst.

1453!! h_dst must be dimensioned as a model array with GV%ke layers while h_src can

1454!! have an arbitrary number of layers specified by nk_src.

1455subroutine ale_remap_scalar(CS, G, GV, nk_src, h_src, s_src, h_dst, s_dst, all_cells, old_remap)

1456 type(remapping_cs), intent(in) :: cs !< Remapping control structure

1457 type(ocean_grid_type), intent(in) :: g !< Ocean grid structure

1458 type(verticalgrid_type), intent(in) :: gv !< Ocean vertical grid structure

1459 integer, intent(in) :: nk_src !< Number of levels on source grid

1460 real, dimension(SZI_(G),SZJ_(G),nk_src), intent(in) :: h_src !< Level thickness of source grid

1461 !! [H ~> m or kg m-2] or other units

1462 !! if H_neglect is provided

1463 real, dimension(SZI_(G),SZJ_(G),nk_src), intent(in) :: s_src !< Scalar on source grid, in arbitrary units [A]

1464 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)),intent(in) :: h_dst !< Level thickness of destination grid in the

1465 !! same units as h_src, often [H ~> m or kg m-2]

1466 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)),intent(inout) :: s_dst !< Scalar on destination grid, in the same

1467 !! arbitrary units as s_src [A]

1468 logical, optional, intent(in) :: all_cells !< If false, only reconstruct for

1469 !! non-vanished cells. Use all vanished

1470 !! layers otherwise (default).

1471 logical, optional, intent(in) :: old_remap !< If true, use the old "remapping_core_w"

1472 !! method, otherwise use "remapping_core_h".

1473 ! Local variables

1474 integer :: i, j, k, n_points

1475 real :: dx(gv%ke+1) ! Change in interface position [H ~> m or kg m-2]

1476 logical :: ignore_vanished_layers, use_remapping_core_w

1477

1478 ignore_vanished_layers = .false.

1479 if (present(all_cells)) ignore_vanished_layers = .not. all_cells

1480 use_remapping_core_w = .false.

1481 if (present(old_remap)) use_remapping_core_w = old_remap

1482 n_points = nk_src

1483

1484 !$OMP parallel do default(shared) firstprivate(n_points,dx)

1485 do j = g%jsc,g%jec ; do i = g%isc,g%iec

1486 if (g%mask2dT(i,j) > 0.) then

1487 if (ignore_vanished_layers) then

1488 n_points = 0

1489 do k = 1, nk_src

1490 if (h_src(i,j,k)>0.) n_points = n_points + 1

1491 enddo

1492 s_dst(i,j,:) = 0.

1493 endif

1494 if (use_remapping_core_w) then

1495 call dzfromh1h2( n_points, h_src(i,j,1:n_points), gv%ke, h_dst(i,j,:), dx )

1496 call remapping_core_w(cs, n_points, h_src(i,j,1:n_points), s_src(i,j,1:n_points), &

1497 gv%ke, dx, s_dst(i,j,:))

1498 else

1499 call remapping_core_h(cs, n_points, h_src(i,j,1:n_points), s_src(i,j,1:n_points), &

1500 gv%ke, h_dst(i,j,:), s_dst(i,j,:))

1501 endif

1502 else

1503 s_dst(i,j,:) = 0.

1504 endif

1505 enddo ; enddo

1506

1507end subroutine ale_remap_scalar

1508

1509

1510!> Calculate edge values (top and bottom of layer) for T and S consistent with a PLM reconstruction

1511!! in the vertical direction. Boundary reconstructions are PCM unless bdry_extrap is true.

1512subroutine ts_plm_edge_values( CS, S_t, S_b, T_t, T_b, G, GV, tv, h, bdry_extrap )

1513 type(ocean_grid_type), intent(in) :: g !< ocean grid structure

1514 type(verticalgrid_type), intent(in) :: gv !< Ocean vertical grid structure

1515 type(ale_cs), intent(inout) :: cs !< module control structure

1516 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &

1517 intent(inout) :: s_t !< Salinity at the top edge of each layer [S ~> ppt]

1518 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &

1519 intent(inout) :: s_b !< Salinity at the bottom edge of each layer [S ~> ppt]

1520 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &

1521 intent(inout) :: t_t !< Temperature at the top edge of each layer [C ~> degC]

1522 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &

1523 intent(inout) :: t_b !< Temperature at the bottom edge of each layer [C ~> degC]

1524 type(thermo_var_ptrs), intent(in) :: tv !< thermodynamics structure

1525 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &

1526 intent(in) :: h !< layer thickness [H ~> m or kg m-2]

1527 logical, intent(in) :: bdry_extrap !< If true, use high-order boundary

1528 !! extrapolation within boundary cells

1529

1530 call ale_plm_edge_values( cs, g, gv, h, tv%S, bdry_extrap, s_t, s_b )

1531 call ale_plm_edge_values( cs, g, gv, h, tv%T, bdry_extrap, t_t, t_b )

1532

1533end subroutine ts_plm_edge_values

1534

1535!> Calculate edge values (top and bottom of layer) 3d scalar array.

1536!! Boundary reconstructions are PCM unless bdry_extrap is true.

1537subroutine ale_plm_edge_values( CS, G, GV, h, Q, bdry_extrap, Q_t, Q_b )

1538 type(ale_cs), intent(in) :: cs !< module control structure

1539 type(ocean_grid_type), intent(in) :: g !< ocean grid structure

1540 type(verticalgrid_type), intent(in) :: gv !< Ocean vertical grid structure

1541 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &

1542 intent(in) :: h !< layer thickness [H ~> m or kg m-2]

1543 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &

1544 intent(in) :: q !< 3d scalar array, in arbitrary units [A]

1545 logical, intent(in) :: bdry_extrap !< If true, use high-order boundary

1546 !! extrapolation within boundary cells

1547 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &

1548 intent(inout) :: q_t !< Scalar at the top edge of each layer [A]

1549 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &

1550 intent(inout) :: q_b !< Scalar at the bottom edge of each layer [A]

1551 ! Local variables

1552 integer :: i, j, k

1553 real :: slp(gv%ke) ! Tracer slope times the cell width [A]

1554 real :: mslp ! Monotonized tracer slope times the cell width [A]

1555 real :: h_neglect ! Tiny thicknesses used in remapping [H ~> m or kg m-2]

1556

1557 if (cs%answer_date >= 20190101) then

1558 h_neglect = gv%H_subroundoff

1559 elseif (gv%Boussinesq) then

1560 h_neglect = gv%m_to_H*1.0e-30

1561 else

1562 h_neglect = gv%kg_m2_to_H*1.0e-30

1563 endif

1564

1565 !$OMP parallel do default(shared) private(slp,mslp)

1566 do j = g%jsc-1,g%jec+1 ; do i = g%isc-1,g%iec+1

1567 slp(1) = 0.

1568 do k = 2, gv%ke-1

1569 slp(k) = plm_slope_wa(h(i,j,k-1), h(i,j,k), h(i,j,k+1), h_neglect, &

1570 q(i,j,k-1), q(i,j,k), q(i,j,k+1))

1571 enddo

1572 slp(gv%ke) = 0.

1573

1574 do k = 2, gv%ke-1

1575 mslp = plm_monotonized_slope(q(i,j,k-1), q(i,j,k), q(i,j,k+1), slp(k-1), slp(k), slp(k+1))

1576 q_t(i,j,k) = q(i,j,k) - 0.5 * mslp

1577 q_b(i,j,k) = q(i,j,k) + 0.5 * mslp

1578 enddo

1579 if (bdry_extrap) then

1580 mslp = - plm_extrapolate_slope(h(i,j,2), h(i,j,1), h_neglect, q(i,j,2), q(i,j,1))

1581 q_t(i,j,1) = q(i,j,1) - 0.5 * mslp

1582 q_b(i,j,1) = q(i,j,1) + 0.5 * mslp

1583 mslp = plm_extrapolate_slope(h(i,j,gv%ke-1), h(i,j,gv%ke), h_neglect, &

1584 q(i,j,gv%ke-1), q(i,j,gv%ke))

1585 q_t(i,j,gv%ke) = q(i,j,gv%ke) - 0.5 * mslp

1586 q_b(i,j,gv%ke) = q(i,j,gv%ke) + 0.5 * mslp

1587 else

1588 q_t(i,j,1) = q(i,j,1)

1589 q_b(i,j,1) = q(i,j,1)

1590 q_t(i,j,gv%ke) = q(i,j,gv%ke)

1591 q_b(i,j,gv%ke) = q(i,j,gv%ke)

1592 endif

1593

1594 enddo ; enddo

1595

1596end subroutine ale_plm_edge_values

1597

1598!> Calculate edge values (top and bottom of layer) for T and S consistent with a PPM reconstruction

1599!! in the vertical direction. Boundary reconstructions are PCM unless bdry_extrap is true.

1600subroutine ts_ppm_edge_values( CS, S_t, S_b, T_t, T_b, G, GV, tv, h, bdry_extrap )

1601 type(ocean_grid_type), intent(in) :: g !< ocean grid structure

1602 type(verticalgrid_type), intent(in) :: gv !< Ocean vertical grid structure

1603 type(ale_cs), intent(inout) :: cs !< module control structure

1604 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &

1605 intent(inout) :: s_t !< Salinity at the top edge of each layer [S ~> ppt]

1606 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &

1607 intent(inout) :: s_b !< Salinity at the bottom edge of each layer [S ~> ppt]

1608 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &

1609 intent(inout) :: t_t !< Temperature at the top edge of each layer [C ~> degC]

1610 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &

1611 intent(inout) :: t_b !< Temperature at the bottom edge of each layer [C ~> degC]

1612 type(thermo_var_ptrs), intent(in) :: tv !< thermodynamics structure

1613 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &

1614 intent(in) :: h !< layer thicknesses [H ~> m or kg m-2]

1615 logical, intent(in) :: bdry_extrap !< If true, use high-order boundary

1616 !! extrapolation within boundary cells

1617

1618 ! Local variables

1619 integer :: i, j, k

1620 real :: htmp(gv%ke) ! A 1-d copy of h [H ~> m or kg m-2]

1621 real :: tmp(gv%ke) ! A 1-d copy of a column of temperature [C ~> degC] or salinity [S ~> ppt]

1622 real, dimension(CS%nk,2) :: &

1623 ppol_e ! Edge value of polynomial in [C ~> degC] or [S ~> ppt]

1624 real, dimension(CS%nk,3) :: &

1625 ppol_coefs ! Coefficients of polynomial, all in [C ~> degC] or [S ~> ppt]

1626 real :: h_neglect, h_neglect_edge ! Tiny thicknesses [H ~> m or kg m-2]

1627

1628 if (cs%answer_date >= 20190101) then

1629 h_neglect = gv%H_subroundoff ; h_neglect_edge = gv%H_subroundoff

1630 elseif (gv%Boussinesq) then

1631 h_neglect = gv%m_to_H*1.0e-30 ; h_neglect_edge = gv%m_to_H*1.0e-10

1632 else

1633 h_neglect = gv%kg_m2_to_H*1.0e-30 ; h_neglect_edge = gv%kg_m2_to_H*1.0e-10

1634 endif

1635

1636 ! Determine reconstruction within each column

1637 !$OMP parallel do default(shared) private(hTmp,tmp,ppol_E,ppol_coefs)

1638 do j = g%jsc-1,g%jec+1 ; do i = g%isc-1,g%iec+1

1639

1640 ! Build current grid

1641 htmp(:) = h(i,j,:)

1642 tmp(:) = tv%S(i,j,:)

1643

1644 ! Reconstruct salinity profile

1645 ppol_e(:,:) = 0.0

1646 ppol_coefs(:,:) = 0.0

1647 call edge_values_implicit_h4( gv%ke, htmp, tmp, ppol_e, h_neglect=h_neglect_edge, &

1648 answer_date=cs%answer_date )

1649 call ppm_reconstruction( gv%ke, htmp, tmp, ppol_e, ppol_coefs, h_neglect, &

1650 answer_date=cs%answer_date )

1651 if (bdry_extrap) &

1652 call ppm_boundary_extrapolation( gv%ke, htmp, tmp, ppol_e, ppol_coefs, h_neglect )

1653

1654 do k = 1,gv%ke

1655 s_t(i,j,k) = ppol_e(k,1)

1656 s_b(i,j,k) = ppol_e(k,2)

1657 enddo

1658

1659 ! Reconstruct temperature profile

1660 ppol_e(:,:) = 0.0

1661 ppol_coefs(:,:) = 0.0

1662 tmp(:) = tv%T(i,j,:)

1663 if (cs%answer_date < 20190101) then

1664 call edge_values_implicit_h4( gv%ke, htmp, tmp, ppol_e, h_neglect=1.0e-10*gv%m_to_H, &

1665 answer_date=cs%answer_date )

1666 else

1667 call edge_values_implicit_h4( gv%ke, htmp, tmp, ppol_e, h_neglect=gv%H_subroundoff, &

1668 answer_date=cs%answer_date )

1669 endif

1670 call ppm_reconstruction( gv%ke, htmp, tmp, ppol_e, ppol_coefs, h_neglect, &

1671 answer_date=cs%answer_date )

1672 if (bdry_extrap) &

1673 call ppm_boundary_extrapolation(gv%ke, htmp, tmp, ppol_e, ppol_coefs, h_neglect )

1674

1675 do k = 1,gv%ke

1676 t_t(i,j,k) = ppol_e(k,1)

1677 t_b(i,j,k) = ppol_e(k,2)

1678 enddo

1679

1680 enddo ; enddo

1681

1682end subroutine ts_ppm_edge_values

1683

1684!> Calculate edge values (top and bottom of layer) for T and S consistent with a PLM reconstruction

1685!! in the vertical direction that uses weighted least squares for the slope.

1686subroutine ts_plm_wls_edge_values(CS, S_t, S_b, T_t, T_b, G, GV, tv, h)

1687 type(ocean_grid_type), intent(in) :: g !< ocean grid structure

1688 type(verticalgrid_type), intent(in) :: gv !< Ocean vertical grid structure

1689 type(ale_cs), intent(inout) :: cs !< module control structure

1690 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &

1691 intent(inout) :: s_t !< Salinity at the top edge of each layer [S ~> ppt]

1692 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &

1693 intent(inout) :: s_b !< Salinity at the bottom edge of each layer [S ~> ppt]

1694 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &

1695 intent(inout) :: t_t !< Temperature at the top edge of each layer [C ~> degC]

1696 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &

1697 intent(inout) :: t_b !< Temperature at the bottom edge of each layer [C ~> degC]

1698 type(thermo_var_ptrs), intent(in) :: tv !< thermodynamics structure

1699 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &

1700 intent(in) :: h !< layer thickness [H ~> m or kg m-2]

1701 ! Local variables

1702 integer :: i, j, k

1703 type(plm_wls) :: recon !< A PLM-WLS reconstruction

1704

1705 call recon%init(gv%ke, h_neglect=gv%H_subroundoff)

1706

1707 !$OMP parallel do default(shared) firstprivate(recon)

1708 do j = g%jsc-1,g%jec+1 ; do i = g%isc-1,g%iec+1

1709

1710 call recon%reconstruct(h(i,j,:), tv%T(i,j,:))

1711 t_t(i,j,:) = recon%ul(:)

1712 t_b(i,j,:) = recon%ur(:)

1713

1714 call recon%reconstruct(h(i,j,:), tv%S(i,j,:))

1715 s_t(i,j,:) = recon%ul(:)

1716 s_b(i,j,:) = recon%ur(:)

1717

1718 enddo ; enddo

1719

1720 call recon%destroy()

1721

1722end subroutine ts_plm_wls_edge_values

1723

1724!> Initializes regridding for the main ALE algorithm

1725subroutine ale_initregridding(G, GV, US, max_depth, param_file, mdl, regridCS)

1726 type(ocean_grid_type), intent(in) :: g !< Grid structure

1727 type(verticalgrid_type), intent(in) :: gv !< Ocean vertical grid structure

1728 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type

1729 real, intent(in) :: max_depth !< The maximum depth of the ocean [Z ~> m].

1730 type(param_file_type), intent(in) :: param_file !< parameter file

1731 character(len=*), intent(in) :: mdl !< Name of calling module

1732 type(regridding_cs), intent(out) :: regridcs !< Regridding parameters and work arrays

1733 ! Local variables

1734 character(len=30) :: coord_mode

1735

1736 call get_param(param_file, mdl, "REGRIDDING_COORDINATE_MODE", coord_mode, &

1737 "Coordinate mode for vertical regridding. "//&

1738 "Choose among the following possibilities: "//&

1739 trim(regriddingcoordinatemodedoc), &

1740 default=default_coordinate_mode, fail_if_missing=.true.)

1741

1742 call initialize_regridding(regridcs, g, gv, us, max_depth, param_file, mdl, coord_mode, '', '')

1743

1744end subroutine ale_initregridding

1745

1746!> Query the target coordinate interfaces positions

1747function ale_getcoordinate( CS )

1748 type(ale_cs), pointer :: cs !< module control structure

1749

1750 real, dimension(CS%nk+1) :: ale_getcoordinate !< The coordinate positions, in the appropriate units

1751 !! of the target coordinate, e.g. [Z ~> m] for z*,

1752 !! non-dimensional for sigma, etc.

1753 ale_getcoordinate(:) = getcoordinateinterfaces( cs%regridCS, undo_scaling=.true. )

1754

1755end function ale_getcoordinate

1756

1757

1758!> Query the target coordinate units

1759function ale_getcoordinateunits( CS )

1760 type(ale_cs), pointer :: cs !< module control structure

1761

1762 character(len=20) :: ale_getcoordinateunits

1763

1764 ale_getcoordinateunits = getcoordinateunits( cs%regridCS )

1765

1766end function ale_getcoordinateunits

1767

1768

1769!> Returns true if initial conditions should be regridded and remapped

1770logical function ale_remap_init_conds( CS )

1771 type(ale_cs), pointer :: cs !< module control structure

1772

1773 ale_remap_init_conds = .false.

1774 if (associated(cs)) ale_remap_init_conds = cs%remap_after_initialization

1775end function ale_remap_init_conds

1776

1777!> Updates the weights for time filtering the new grid generated in regridding

1778subroutine ale_update_regrid_weights( dt, CS )

1779 real, intent(in) :: dt !< Time-step used between ALE calls [T ~> s]

1780 type(ale_cs), pointer :: cs !< ALE control structure

1781 ! Local variables

1782 real :: w ! An implicit weighting estimate [nondim]

1783

1784 if (associated(cs)) then

1785 w = 0.0

1786 if (cs%regrid_time_scale > 0.0) then

1787 w = cs%regrid_time_scale / (cs%regrid_time_scale + dt)

1788 endif

1789 call set_regrid_params(cs%regridCS, old_grid_weight=w)

1790 endif

1791

1792end subroutine ale_update_regrid_weights

1793

1794!> Update the vertical grid type with ALE information.

1795!! This subroutine sets information in the verticalGrid_type to be

1796!! consistent with the use of ALE mode.

1797subroutine ale_updateverticalgridtype(CS, GV)

1798 type(ale_cs), pointer :: cs !< ALE control structure

1799 type(verticalgrid_type), pointer :: gv !< vertical grid information

1800

1801 integer :: nk

1802

1803 nk = gv%ke

1804 gv%sInterface(1:nk+1) = getcoordinateinterfaces( cs%regridCS, undo_scaling=.true. )

1805 gv%sLayer(1:nk) = 0.5*( gv%sInterface(1:nk) + gv%sInterface(2:nk+1) )

1806 gv%zAxisUnits = getcoordinateunits( cs%regridCS )

1807 gv%zAxisLongName = getcoordinateshortname( cs%regridCS )

1808 gv%direction = -1 ! Because of ferret in z* mode. Need method to set

1809 ! as function of coordinate mode.

1810

1811end subroutine ale_updateverticalgridtype

1812

1813

1814!> Write the vertical coordinate information into a file.

1815!! This subroutine writes out a file containing any available data related

1816!! to the vertical grid used by the MOM ocean model when in ALE mode.

1817subroutine ale_writecoordinatefile( CS, GV, directory )

1818 type(ale_cs), pointer :: cs !< module control structure

1819 type(verticalgrid_type), intent(in) :: gv !< ocean vertical grid structure

1820 character(len=*), intent(in) :: directory !< directory for writing grid info

1821

1822 character(len=240) :: filepath

1823

1824 filepath = trim(directory) // trim("Vertical_coordinate.nc")

1825

1826 call write_regrid_file(cs%regridCS, gv, filepath)

1827

1828end subroutine ale_writecoordinatefile

1829

1830!> Set h to coordinate values for fixed coordinate systems

1831subroutine ale_initthicknesstocoord( CS, G, GV, h, height_units )

1832 type(ale_cs), intent(inout) :: cs !< module control structure

1833 type(ocean_grid_type), intent(in) :: g !< module grid structure

1834 type(verticalgrid_type), intent(in) :: gv !< Ocean vertical grid structure

1835 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: h !< layer thickness in thickness units

1836 !! [H ~> m or kg m-2] or height units [Z ~> m]

1837 logical, optional, intent(in) :: height_units !< If present and true, the

1838 !! thicknesses are in height units

1839

1840 ! Local variables

1841 real :: scale ! A scaling value for the thicknesses [nondim] or [H Z-1 ~> nondim or kg m-3]

1842 integer :: i, j

1843

1844 scale = gv%Z_to_H

1845 if (present(height_units)) then ; if (height_units) scale = 1.0 ; endif

1846 do j = g%jsd,g%jed ; do i = g%isd,g%ied

1847 h(i,j,:) = scale * getstaticthickness( cs%regridCS, 0., max(g%meanSL(i,j)+g%bathyT(i,j), 0.0) )

1848 enddo ; enddo

1849

1850end subroutine ale_initthicknesstocoord

1851

1852end module mom_ale