MOM_remapping.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!> Provides column-wise vertical remapping functions

6module mom_remapping

7

8! Original module written by Laurent White, 2008.06.09

9

10use mom_error_handler, only : mom_error, fatal

11use mom_string_functions, only : uppercase

12use numerical_testing_type, only : testing

13use regrid_edge_values, only : edge_values_explicit_h4, edge_values_implicit_h4

14use regrid_edge_values, only : edge_values_explicit_h4cw

15use regrid_edge_values, only : edge_values_implicit_h4, edge_values_implicit_h6

16use regrid_edge_values, only : edge_slopes_implicit_h3, edge_slopes_implicit_h5

17use pcm_functions, only : pcm_reconstruction

18use plm_functions, only : plm_reconstruction, plm_boundary_extrapolation

19use ppm_functions, only : ppm_reconstruction, ppm_boundary_extrapolation

20use ppm_functions, only : ppm_monotonicity

21use pqm_functions, only : pqm_reconstruction, pqm_boundary_extrapolation_v1

22use mom_hybgen_remap, only : hybgen_plm_coefs, hybgen_ppm_coefs, hybgen_weno_coefs

23

24use recon1d_type, only : recon1d

25use recon1d_pcm, only : pcm

26use recon1d_plm_cw, only : plm_cw

27use recon1d_plm_hybgen, only : plm_hybgen

28use recon1d_plm_cwk, only : plm_cwk

29use recon1d_mplm_cwk, only : mplm_cwk

30use recon1d_emplm_cwk, only : emplm_cwk

31use recon1d_mplm_wa, only : mplm_wa

32use recon1d_emplm_wa, only : emplm_wa

33use recon1d_mplm_wa_poly, only : mplm_wa_poly

34use recon1d_emplm_wa_poly, only : emplm_wa_poly

35use recon1d_ppm_cw, only : ppm_cw

36use recon1d_ppm_hybgen, only : ppm_hybgen

37use recon1d_ppm_cwk, only : ppm_cwk

38use recon1d_eppm_cwk, only : eppm_cwk

39use recon1d_ppm_h4_2019, only : ppm_h4_2019

40use recon1d_ppm_h4_2018, only : ppm_h4_2018

41use recon1d_plm_wls, only : plm_wls

42

43implicit none ; private

44

45!> Container for remapping parameters

46type, public :: remapping_cs ; private

47 !> Determines which reconstruction to use

48 integer :: remapping_scheme = -911

49 !> Degree of polynomial reconstruction

50 integer :: degree = 0

51 !> If true, extrapolate boundaries

52 logical :: boundary_extrapolation = .true.

53 !> If true, reconstructions are checked for consistency.

54 logical :: check_reconstruction = .false.

55 !> If true, the result of remapping are checked for conservation and bounds.

56 logical :: check_remapping = .false.

57 !> If true, the intermediate values used in remapping are forced to be bounded.

58 logical :: force_bounds_in_subcell = .false.

59 !> If true, impose bounds on the remapping from sub-cells to target grid

60 logical :: force_bounds_in_target = .true.

61 !> If true, impose bounds on the remapping from non-vanished sub-cells to target grid

62 logical :: better_force_bounds_in_target = .false.

63 !> If true, calculate and use an offset when summing sub-cells to the target grid

64 logical :: offset_tgt_summation = .false.

65 !> The vintage of the expressions to use for remapping. Values below 20190101 result

66 !! in the use of older, less accurate expressions.

67 integer :: answer_date = 99991231

68 !> If true, use the OM4 version of the remapping algorithm that makes poor assumptions

69 !! about the reconstructions in top and bottom layers of the source grid

70 logical :: om4_remap_via_sub_cells = .false.

71

72 !> A negligibly small width for the purpose of cell reconstructions in the same units

73 !! as the h0 argument to remapping_core_h [H]

74 real :: h_neglect

75 !> A negligibly small width for the purpose of edge value calculations in the same units

76 !! as the h0 argument to remapping_core_h [H]

77 real :: h_neglect_edge

78

79 !> If true, do some debugging as operations proceed

80 logical :: debug = .false.

81

82 !> The instance of the actual equation of state

83 class(recon1d), pointer :: reconstruction => null()

84end type

85

86! The following routines are visible to the outside world

87public remapping_core_h, remapping_core_w

88public initialize_remapping, end_remapping, remapping_set_param, extract_member_remapping_cs

89public remapping_unit_tests, build_reconstructions_1d, average_value_ppoly

90public interpolate_column, reintegrate_column, dzfromh1h2

91

92! The following are private parameter constants

93integer, parameter :: remapping_pcm = 0 !< O(h^1) remapping scheme

94integer, parameter :: remapping_plm = 2 !< O(h^2) remapping scheme

95integer, parameter :: remapping_plm_hybgen = 3 !< O(h^2) remapping scheme

96integer, parameter :: remapping_ppm_cw =10 !< O(h^3) remapping scheme

97integer, parameter :: remapping_ppm_h4 = 4 !< O(h^3) remapping scheme

98integer, parameter :: remapping_ppm_ih4 = 5 !< O(h^3) remapping scheme

99integer, parameter :: remapping_ppm_hybgen = 6 !< O(h^3) remapping scheme

100integer, parameter :: remapping_weno_hybgen= 7 !< O(h^3) remapping scheme

101integer, parameter :: remapping_pqm_ih4ih3 = 8 !< O(h^4) remapping scheme

102integer, parameter :: remapping_pqm_ih6ih5 = 9 !< O(h^5) remapping scheme

103integer, parameter :: remapping_via_class =99 !< Scheme is controlled by Recon1d class

104

105integer, parameter :: integration_pcm = 0 !< Piecewise Constant Method

106integer, parameter :: integration_plm = 1 !< Piecewise Linear Method

107integer, parameter :: integration_ppm = 3 !< Piecewise Parabolic Method

108integer, parameter :: integration_pqm = 5 !< Piecewise Quartic Method

109

110!> Documentation for external callers

111character(len=360), public :: remappingschemesdoc = &

112 "PCM (1st-order accurate)\n"//&

113 "PLM (2nd-order accurate)\n"//&

114 "PLM_HYBGEN (2nd-order accurate)\n"//&

115 "PPM_H4 (3rd-order accurate)\n"//&

116 "PPM_IH4 (3rd-order accurate)\n"//&

117 "PPM_HYBGEN (3rd-order accurate)\n"//&

118 "WENO_HYBGEN (3rd-order accurate)\n"//&

119 "PQM_IH4IH3 (4th-order accurate)\n"//&

120 "PQM_IH6IH5 (5th-order accurate)\n"

121character(len=3), public :: remappingdefaultscheme = "PLM" !< Default remapping method

122

123contains

124

125!> Set parameters within remapping object

126subroutine remapping_set_param(CS, remapping_scheme, boundary_extrapolation, &

127 check_reconstruction, check_remapping, force_bounds_in_subcell, &

128 force_bounds_in_target, better_force_bounds_in_target, offset_tgt_summation, &

129 om4_remap_via_sub_cells, answers_2018, answer_date, nk, &

130 h_neglect, h_neglect_edge)

131 type(remapping_cs), intent(inout) :: cs !< Remapping control structure

132 character(len=*), optional, intent(in) :: remapping_scheme !< Remapping scheme to use

133 logical, optional, intent(in) :: boundary_extrapolation !< Indicate to extrapolate in boundary cells

134 logical, optional, intent(in) :: check_reconstruction !< Indicate to check reconstructions

135 logical, optional, intent(in) :: check_remapping !< Indicate to check results of remapping

136 logical, optional, intent(in) :: force_bounds_in_subcell !< Force subcells values to be bounded

137 logical, optional, intent(in) :: force_bounds_in_target !< Force target values to be bounded

138 logical, optional, intent(in) :: better_force_bounds_in_target !< Force target values to be bounded

139 logical, optional, intent(in) :: offset_tgt_summation !< Use an offset when summing sub-cells

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

141 logical, optional, intent(in) :: answers_2018 !< If true use older, less accurate expressions.

142 integer, optional, intent(in) :: answer_date !< The vintage of the expressions to use

143 real, optional, intent(in) :: h_neglect !< A negligibly small width for the purpose of cell

144 !! reconstructions in the same units as the h0 argument

145 !! to remapping_core_h [H]

146 real, optional, intent(in) :: h_neglect_edge !< A negligibly small width for the purpose of edge

147 !! value calculations in the same units as as the h0

148 !! argument to remapping_core_h [H]

149 integer, optional, intent(in) :: nk !< Number of levels to initialize reconstruction class with

150

151 if (present(remapping_scheme)) then

152 call setreconstructiontype( remapping_scheme, cs )

153 if (index(trim(remapping_scheme),'C_')>0) then

154 if (present(nk)) then

155 call cs%reconstruction%init(nk, h_neglect=h_neglect)

156 else

157 call mom_error( fatal, 'MOM_remapping, remapping_set_param: '//&

158 'Using the Recon1d class for remapping requires nk to be passed' )

159 endif

160 endif

161 endif

162 if (present(boundary_extrapolation)) then

163 cs%boundary_extrapolation = boundary_extrapolation

164 endif

165 if (present(check_reconstruction)) then

166 cs%check_reconstruction = check_reconstruction

167 endif

168 if (present(check_remapping)) then

169 cs%check_remapping = check_remapping

170 endif

171 if (present(force_bounds_in_subcell)) then

172 cs%force_bounds_in_subcell = force_bounds_in_subcell

173 endif

174 if (present(force_bounds_in_target)) then

175 cs%force_bounds_in_target = force_bounds_in_target

176 endif

177 if (present(better_force_bounds_in_target)) then

178 cs%better_force_bounds_in_target = better_force_bounds_in_target

179 endif

180 if (present(offset_tgt_summation)) then

181 cs%offset_tgt_summation = offset_tgt_summation

182 endif

183 if (present(om4_remap_via_sub_cells)) then

184 cs%om4_remap_via_sub_cells = om4_remap_via_sub_cells

185 endif

186 if (present(answers_2018)) then

187 if (answers_2018) then

188 cs%answer_date = 20181231

189 else

190 cs%answer_date = 20190101

191 endif

192 endif

193 if (present(answer_date)) then

194 cs%answer_date = answer_date

195 endif

196 if (present(h_neglect)) then

197 cs%h_neglect = h_neglect

198 endif

199 if (present(h_neglect_edge)) then

200 cs%h_neglect_edge = h_neglect_edge

201 endif

202

203end subroutine remapping_set_param

204

205subroutine extract_member_remapping_cs(CS, remapping_scheme, degree, boundary_extrapolation, check_reconstruction, &

206 check_remapping, force_bounds_in_subcell, force_bounds_in_target, &

207 better_force_bounds_in_target, offset_tgt_summation)

208 type(remapping_cs), intent(in) :: cs !< Control structure for remapping module

209 integer, optional, intent(out) :: remapping_scheme !< Determines which reconstruction scheme to use

210 integer, optional, intent(out) :: degree !< Degree of polynomial reconstruction

211 logical, optional, intent(out) :: boundary_extrapolation !< If true, extrapolate boundaries

212 logical, optional, intent(out) :: check_reconstruction !< If true, reconstructions are checked for consistency.

213 logical, optional, intent(out) :: check_remapping !< If true, the result of remapping are checked

214 !! for conservation and bounds.

215 logical, optional, intent(out) :: force_bounds_in_subcell !< If true, the intermediate values used in

216 !! remapping are forced to be bounded.

217 logical, optional, intent(out) :: force_bounds_in_target !< Force target values to be bounded

218 logical, optional, intent(out) :: better_force_bounds_in_target !< Force target values to be bounded

219 logical, optional, intent(out) :: offset_tgt_summation !< Use an offset when summing sub-cells

220

221 if (present(remapping_scheme)) remapping_scheme = cs%remapping_scheme

222 if (present(degree)) degree = cs%degree

223 if (present(boundary_extrapolation)) boundary_extrapolation = cs%boundary_extrapolation

224 if (present(check_reconstruction)) check_reconstruction = cs%check_reconstruction

225 if (present(check_remapping)) check_remapping = cs%check_remapping

226 if (present(force_bounds_in_subcell)) force_bounds_in_subcell = cs%force_bounds_in_subcell

227 if (present(force_bounds_in_target)) force_bounds_in_target = cs%force_bounds_in_target

228 if (present(better_force_bounds_in_target)) better_force_bounds_in_target = cs%better_force_bounds_in_target

229 if (present(offset_tgt_summation)) offset_tgt_summation = cs%offset_tgt_summation

230

231end subroutine extract_member_remapping_cs

232

233!> Remaps column of values u0 on grid h0 to grid h1 assuming the top edge is aligned and using the OM4

234!! reconstruction methods

235!!

236!! \todo Remove h_neglect argument by moving into remapping_CS

237!! \todo Remove PCM_cell argument by adding new method in Recon1D class

238subroutine remapping_core_h(CS, n0, h0, u0, n1, h1, u1, net_err, PCM_cell)

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

240 integer, intent(in) :: n0 !< Number of cells on source grid

241 real, dimension(n0), intent(in) :: h0 !< Cell widths on source grid [H]

242 real, dimension(n0), intent(in) :: u0 !< Cell averages on source grid [A]

243 integer, intent(in) :: n1 !< Number of cells on target grid

244 real, dimension(n1), intent(in) :: h1 !< Cell widths on target grid [H]

245 real, dimension(n1), intent(out) :: u1 !< Cell averages on target grid [A]

246 real, optional, intent(out) :: net_err !< Error in total column [A H]

247 logical, dimension(n0), optional, intent(in) :: pcm_cell !< If present, use PCM remapping for

248 !! cells in the source grid where this is true.

249 ! Local variables

250 real, dimension(n0+n1+1) :: h_sub ! Width of each each sub-cell [H]

251 real, dimension(n0+n1+1) :: uh_sub ! Integral of u*h over each sub-cell [A H]

252 real, dimension(n0+n1+1) :: u_sub ! Average of u over each sub-cell [A]

253 integer, dimension(n0+n1+1) :: isub_src ! Index of source cell for each sub-cell

254 integer, dimension(n0) :: isrc_start ! Index of first sub-cell within each source cell

255 integer, dimension(n0) :: isrc_end ! Index of last sub-cell within each source cell

256 integer, dimension(n0) :: isrc_max ! Index of thickest sub-cell within each source cell

257 real, dimension(n0) :: h0_eff ! Effective thickness of source cells [H]

258 integer, dimension(n1) :: itgt_start ! Index of first sub-cell within each target cell

259 integer, dimension(n1) :: itgt_end ! Index of last sub-cell within each target cell

260 ! For error checking/debugging

261 real :: u02_err ! Integrated reconstruction error estimates [H A]

262 real, dimension(n0,2) :: ppoly_r_e ! Edge value of polynomial [A]

263 real, dimension(n0,2) :: ppoly_r_s ! Edge slope of polynomial [A H-1]

264 real, dimension(n0,CS%degree+1) :: ppoly_r_coefs ! Coefficients of polynomial reconstructions [A]

265 real :: uh_err ! A bound on the error in the sum of u*h, as estimated by the remapping code [A H]

266 integer :: imethod ! An integer indicating the integration method used

267

268 ! Calculate sub-layer thicknesses and indices connecting sub-layers to source and target grids

269 ! Sets: h_sub, h0_eff, isrc_start, isrc_end, isrc_max, isub_src, itgt_start, itgt_end

270 call intersect_src_tgt_grids(n0, h0, n1, h1, h_sub, h0_eff, &

271 isrc_start, isrc_end, isrc_max, itgt_start, itgt_end, isub_src)

272

273 if (cs%remapping_scheme == remapping_via_class) then

274

275! if (CS%debug) call CS%reconstruction%set_debug() ! Sets an internal flag

276

277 call cs%reconstruction%reconstruct(h0, u0)

278

279 ! Adjust h_sub so that the Hallberg conservation trick works properly

280! call adjust_h_sub( n0, h0, n1, isrc_start, isrc_end, isrc_max, h_sub )

281

282 ! Loop over each sub-cell to calculate average/integral values within each sub-cell.

283 ! Uses: h_sub, isrc_start, isrc_end, isrc_max, isub_src

284 ! Sets: u_sub, uh_sub

285 call cs%reconstruction%remap_to_sub_grid(h0, u0, n1, h_sub, &

286 isrc_start, isrc_end, isrc_max, isub_src, &

287 u_sub, uh_sub, u02_err)

288

289 ! Loop over each target cell summing the integrals from sub-cells within the target cell.

290 ! Uses: itgt_start, itgt_end, h1, h_sub, uh_sub, u_sub

291 ! Sets: u1, uh_err

292 call remap_sub_to_tgt_grid(n0, n1, h1, h_sub, u_sub, uh_sub, itgt_start, itgt_end, &

293 cs%force_bounds_in_target, cs%offset_tgt_summation, &

294 cs%better_force_bounds_in_target, u1, uh_err)

295

296 ! Include the error remapping from source to sub-cells in the estimate of total remapping error

297 uh_err = uh_err + u02_err

298

299 else ! Uses the OM4-era reconstruction functions

300

301 call build_reconstructions_1d(cs, n0, h0, u0, ppoly_r_coefs, ppoly_r_e, ppoly_r_s, imethod, &

302 cs%h_neglect, cs%h_neglect_edge, pcm_cell, debug=cs%debug)

303

304 if (cs%check_reconstruction) call check_reconstructions_1d(n0, h0, u0, cs%degree, &

305 cs%boundary_extrapolation, ppoly_r_coefs, ppoly_r_e)

306

307 ! Loop over each sub-cell to calculate average/integral values within each sub-cell.

308 ! Uses: h_sub, h0_eff, isub_src

309 ! Sets: u_sub, uh_sub

310 if (cs%om4_remap_via_sub_cells) then ! Uses the version from OM4 with a bug at the bottom

311

312 call remap_src_to_sub_grid_om4(n0, h0, u0, ppoly_r_e, ppoly_r_coefs, n1, h_sub, &

313 h0_eff, isrc_start, isrc_end, isrc_max, isub_src, &

314 imethod, cs%force_bounds_in_subcell, u_sub, uh_sub, u02_err)

315

316 else ! i.e. if (CS%om4_remap_via_sub_cells == .false.)

317

318 call remap_src_to_sub_grid(n0, h0, u0, ppoly_r_e, ppoly_r_coefs, n1, h_sub, &

319 isrc_start, isrc_end, isrc_max, isub_src, &

320 imethod, cs%force_bounds_in_subcell, u_sub, uh_sub, u02_err)

321

322 endif

323

324 ! Loop over each target cell summing the integrals from sub-cells within the target cell.

325 ! Uses: itgt_start, itgt_end, h1, h_sub, uh_sub, u_sub

326 ! Sets: u1, uh_err

327 call remap_sub_to_tgt_grid_om4(n0, n1, h1, h_sub, u_sub, uh_sub, itgt_start, itgt_end, &

328 cs%force_bounds_in_target, u1, uh_err)

329 ! Include the error remapping from source to sub-cells in the estimate of total remapping error

330 uh_err = uh_err + u02_err

331

332 if (cs%check_remapping) call check_remapped_values(n0, h0, u0, ppoly_r_e, cs%degree, ppoly_r_coefs, &

333 n1, h1, u1, imethod, uh_err, "remapping_core_h")

334

335 endif

336

337 if (present(net_err)) net_err = uh_err

338

339end subroutine remapping_core_h

340

341!> Remaps column of values u0 on grid h0 to implied grid h1

342!! where the interfaces of h1 differ from those of h0 by dx.

343subroutine remapping_core_w( CS, n0, h0, u0, n1, dx, u1)

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

345 integer, intent(in) :: n0 !< Number of cells on source grid

346 real, dimension(n0), intent(in) :: h0 !< Cell widths on source grid [H]

347 real, dimension(n0), intent(in) :: u0 !< Cell averages on source grid [A]

348 integer, intent(in) :: n1 !< Number of cells on target grid

349 real, dimension(n1+1), intent(in) :: dx !< Cell widths on target grid [H]

350 real, dimension(n1), intent(out) :: u1 !< Cell averages on target grid [A]

351

352 ! Local variables

353 real, dimension(n0+n1+1) :: h_sub ! Width of each each sub-cell [H]

354 real, dimension(n0+n1+1) :: uh_sub ! Integral of u*h over each sub-cell [A H]

355 real, dimension(n0+n1+1) :: u_sub ! Average of u over each sub-cell [A]

356 integer, dimension(n0+n1+1) :: isub_src ! Index of source cell for each sub-cell

357 integer, dimension(n0) :: isrc_start ! Index of first sub-cell within each source cell

358 integer, dimension(n0) :: isrc_end ! Index of last sub-cell within each source cell

359 integer, dimension(n0) :: isrc_max ! Index of thickest sub-cell within each source cell

360 real, dimension(n0) :: h0_eff ! Effective thickness of source cells [H]

361 integer, dimension(n1) :: itgt_start ! Index of first sub-cell within each target cell

362 integer, dimension(n1) :: itgt_end ! Index of last sub-cell within each target cell

363 ! For error checking/debugging

364 real :: u02_err ! Integrated reconstruction error estimates [H A]

365 real, dimension(n0,2) :: ppoly_r_e ! Edge value of polynomial [A]

366 real, dimension(n0,2) :: ppoly_r_s ! Edge slope of polynomial [A H-1]

367 real, dimension(n0,CS%degree+1) :: ppoly_r_coefs ! Coefficients of polynomial reconstructions [A]

368 real, dimension(n1) :: h1 !< Cell widths on target grid [H]

369 real :: uh_err ! A bound on the error in the sum of u*h, as estimated by the remapping code [A H]

370 integer :: imethod ! An integer indicating the integration method used

371 integer :: k

372

373 call build_reconstructions_1d( cs, n0, h0, u0, ppoly_r_coefs, ppoly_r_e, ppoly_r_s, imethod,&

374 cs%h_neglect, cs%h_neglect_edge )

375

376 if (cs%check_reconstruction) call check_reconstructions_1d(n0, h0, u0, cs%degree, &

377 cs%boundary_extrapolation, ppoly_r_coefs, ppoly_r_e)

378

379 ! This is a temporary step prior to switching to remapping_core_h()

380 do k = 1, n1

381 if (k<=n0) then

382 h1(k) = max( 0., h0(k) + ( dx(k+1) - dx(k) ) )

383 else

384 h1(k) = max( 0., dx(k+1) - dx(k) )

385 endif

386 enddo

387

388 ! Calculate sub-layer thicknesses and indices connecting sub-layers to source and target grids

389 call intersect_src_tgt_grids( n0, h0, n1, h1, h_sub, h0_eff, &

390 isrc_start, isrc_end, isrc_max, itgt_start, itgt_end, isub_src )

391

392 ! Loop over each sub-cell to calculate average/integral values within each sub-cell.

393 ! Uses: h_sub, h0_eff, isub_src

394 ! Sets: u_sub, uh_sub

395 call remap_src_to_sub_grid_om4(n0, h0, u0, ppoly_r_e, ppoly_r_coefs, n1, h_sub, &

396 h0_eff, isrc_start, isrc_end, isrc_max, isub_src, &

397 imethod, cs%force_bounds_in_subcell, u_sub, uh_sub, u02_err)

398

399 ! Loop over each target cell summing the integrals from sub-cells within the target cell.

400 ! Uses: itgt_start, itgt_end, h1, h_sub, uh_sub, u_sub

401 ! Sets: u1, uh_err

402 call remap_sub_to_tgt_grid_om4(n0, n1, h1, h_sub, u_sub, uh_sub, itgt_start, itgt_end, &

403 cs%force_bounds_in_target, u1, uh_err)

404

405 ! Include the error remapping from source to sub-cells in the estimate of total remapping error

406 uh_err = uh_err + u02_err

407

408 if (cs%check_remapping) call check_remapped_values(n0, h0, u0, ppoly_r_e, cs%degree, ppoly_r_coefs, &

409 n1, h1, u1, imethod, uh_err, "remapping_core_w")

410

411end subroutine remapping_core_w

412

413!> Creates polynomial reconstructions of u0 on the source grid h0.

414subroutine build_reconstructions_1d( CS, n0, h0, u0, ppoly_r_coefs, &

415 ppoly_r_E, ppoly_r_S, iMethod, h_neglect, &

416 h_neglect_edge, PCM_cell, debug )

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

418 integer, intent(in) :: n0 !< Number of cells on source grid

419 real, dimension(n0), intent(in) :: h0 !< Cell widths on source grid [H]

420 real, dimension(n0), intent(in) :: u0 !< Cell averages on source grid [A]

421 real, dimension(n0,CS%degree+1), &

422 intent(out) :: ppoly_r_coefs !< Coefficients of polynomial [A]

423 real, dimension(n0,2), intent(out) :: ppoly_r_e !< Edge value of polynomial [A]

424 real, dimension(n0,2), intent(out) :: ppoly_r_s !< Edge slope of polynomial [A H-1]

425 integer, intent(out) :: imethod !< Integration method

426 real, intent(in) :: h_neglect !< A negligibly small width for the

427 !! purpose of cell reconstructions

428 !! in the same units as h0 [H]

429 real, optional, intent(in) :: h_neglect_edge !< A negligibly small width for the purpose

430 !! of edge value calculations in the same units as h0 [H].

431 !! The default is h_neglect.

432 logical, optional, intent(in) :: pcm_cell(n0) !< If present, use PCM remapping for

433 !! cells from the source grid where this is true.

434 logical, optional, intent(in) :: debug !< If true, enable debugging

435

436 ! Local variables

437 real :: h_neg_edge ! A negligibly small width for the purpose of edge value

438 ! calculations in the same units as h0 [H]

439 integer :: local_remapping_scheme

440 integer :: k, n

441 logical :: deb ! Do debugging

442

443 deb = .false. ; if (present(debug)) deb = debug

444

445 h_neg_edge = h_neglect ; if (present(h_neglect_edge)) h_neg_edge = h_neglect_edge

446

447 ! Reset polynomial

448 ppoly_r_e(:,:) = 0.0

449 ppoly_r_s(:,:) = 0.0

450 ppoly_r_coefs(:,:) = 0.0

451 imethod = -999

452

453 local_remapping_scheme = cs%remapping_scheme

454 if (n0<=1) then

455 local_remapping_scheme = remapping_pcm

456 elseif (n0<=3) then

457 local_remapping_scheme = min( local_remapping_scheme, remapping_plm )

458 elseif (n0<=4 .and. local_remapping_scheme /= remapping_ppm_cw ) then

459 local_remapping_scheme = min( local_remapping_scheme, remapping_ppm_h4 )

460 endif

461 select case ( local_remapping_scheme )

462 case ( remapping_pcm )

463 call pcm_reconstruction( n0, u0, ppoly_r_e, ppoly_r_coefs)

464 imethod = integration_pcm

465 case ( remapping_plm )

466 call plm_reconstruction( n0, h0, u0, ppoly_r_e, ppoly_r_coefs, h_neglect )

467 if ( cs%boundary_extrapolation ) then

468 call plm_boundary_extrapolation( n0, h0, u0, ppoly_r_e, ppoly_r_coefs, h_neglect)

469 endif

470 imethod = integration_plm

471 case ( remapping_plm_hybgen )

472 call hybgen_plm_coefs(u0, h0, ppoly_r_coefs(:,2), n0, 1, h_neglect)

473 do k=1,n0

474 ppoly_r_e(k,1) = u0(k) - 0.5 * ppoly_r_coefs(k,2) ! Left edge value of cell k

475 ppoly_r_e(k,2) = u0(k) + 0.5 * ppoly_r_coefs(k,2) ! Right edge value of cell k

476 ppoly_r_coefs(k,1) = ppoly_r_e(k,1)

477 enddo

478 if ( cs%boundary_extrapolation ) &

479 call plm_boundary_extrapolation( n0, h0, u0, ppoly_r_e, ppoly_r_coefs, h_neglect )

480 imethod = integration_plm

481 case ( remapping_ppm_cw )

482 ! identical to REMAPPING_PPM_HYBGEN

483 call edge_values_explicit_h4cw( n0, h0, u0, ppoly_r_e, h_neg_edge )

484 call ppm_monotonicity( n0, u0, ppoly_r_e )

485 call ppm_reconstruction( n0, h0, u0, ppoly_r_e, ppoly_r_coefs, h_neglect, answer_date=cs%answer_date )

486 if ( cs%boundary_extrapolation ) then

487 call ppm_boundary_extrapolation( n0, h0, u0, ppoly_r_e, ppoly_r_coefs, h_neglect )

488 endif

489 imethod = integration_ppm

490 case ( remapping_ppm_h4 )

491 call edge_values_explicit_h4( n0, h0, u0, ppoly_r_e, h_neg_edge, answer_date=cs%answer_date )

492 call ppm_reconstruction( n0, h0, u0, ppoly_r_e, ppoly_r_coefs, h_neglect, answer_date=cs%answer_date )

493 if ( cs%boundary_extrapolation ) then

494 call ppm_boundary_extrapolation( n0, h0, u0, ppoly_r_e, ppoly_r_coefs, h_neglect )

495 endif

496 imethod = integration_ppm

497 case ( remapping_ppm_ih4 )

498 call edge_values_implicit_h4( n0, h0, u0, ppoly_r_e, h_neg_edge, answer_date=cs%answer_date )

499 call ppm_reconstruction( n0, h0, u0, ppoly_r_e, ppoly_r_coefs, h_neglect, answer_date=cs%answer_date )

500 if ( cs%boundary_extrapolation ) then

501 call ppm_boundary_extrapolation( n0, h0, u0, ppoly_r_e, ppoly_r_coefs, h_neglect )

502 endif

503 imethod = integration_ppm

504 case ( remapping_ppm_hybgen )

505 call hybgen_ppm_coefs(u0, h0, ppoly_r_e, n0, 1, h_neglect)

506 call ppm_reconstruction( n0, h0, u0, ppoly_r_e, ppoly_r_coefs, h_neglect, answer_date=99991231 )

507 if ( cs%boundary_extrapolation ) &

508 call ppm_boundary_extrapolation( n0, h0, u0, ppoly_r_e, ppoly_r_coefs, h_neglect )

509 imethod = integration_ppm

510 case ( remapping_weno_hybgen )

511 call hybgen_weno_coefs(u0, h0, ppoly_r_e, n0, 1, h_neglect)

512 call ppm_reconstruction( n0, h0, u0, ppoly_r_e, ppoly_r_coefs, h_neglect, answer_date=99991231 )

513 if ( cs%boundary_extrapolation ) &

514 call ppm_boundary_extrapolation( n0, h0, u0, ppoly_r_e, ppoly_r_coefs, h_neglect )

515 imethod = integration_ppm

516 case ( remapping_pqm_ih4ih3 )

517 call edge_values_implicit_h4( n0, h0, u0, ppoly_r_e, h_neg_edge, answer_date=cs%answer_date )

518 call edge_slopes_implicit_h3( n0, h0, u0, ppoly_r_s, h_neglect, answer_date=cs%answer_date )

519 call pqm_reconstruction( n0, h0, u0, ppoly_r_e, ppoly_r_s, ppoly_r_coefs, h_neglect, &

520 answer_date=cs%answer_date )

521 if ( cs%boundary_extrapolation ) then

522 call pqm_boundary_extrapolation_v1( n0, h0, u0, ppoly_r_e, ppoly_r_s, &

523 ppoly_r_coefs, h_neglect )

524 endif

525 imethod = integration_pqm

526 case ( remapping_pqm_ih6ih5 )

527 call edge_values_implicit_h6( n0, h0, u0, ppoly_r_e, h_neg_edge, answer_date=cs%answer_date )

528 call edge_slopes_implicit_h5( n0, h0, u0, ppoly_r_s, h_neglect, answer_date=cs%answer_date )

529 call pqm_reconstruction( n0, h0, u0, ppoly_r_e, ppoly_r_s, ppoly_r_coefs, h_neglect, &

530 answer_date=cs%answer_date )

531 if ( cs%boundary_extrapolation ) then

532 call pqm_boundary_extrapolation_v1( n0, h0, u0, ppoly_r_e, ppoly_r_s, &

533 ppoly_r_coefs, h_neglect )

534 endif

535 imethod = integration_pqm

536 case ( remapping_via_class )

537 call mom_error( fatal, 'MOM_remapping, build_reconstructions_1d: '//&

538 'Should not reach this point if using Recon1d class for remapping' )

539 case default

540 call mom_error( fatal, 'MOM_remapping, build_reconstructions_1d: '//&

541 'The selected remapping method is invalid' )

542 end select

543

544 if (present(pcm_cell)) then

545 ! Change the coefficients to those for the piecewise constant method in indicated cells.

546 do k=1,n0 ; if (pcm_cell(k)) then

547 ppoly_r_coefs(k,1) = u0(k)

548 ppoly_r_e(k,1:2) = u0(k)

549 ppoly_r_s(k,1:2) = 0.0

550 do n=2,cs%degree+1 ; ppoly_r_coefs(k,n) = 0.0 ; enddo

551 endif ; enddo

552 endif

553

554end subroutine build_reconstructions_1d

555

556!> Checks that edge values and reconstructions satisfy bounds

557subroutine check_reconstructions_1d(n0, h0, u0, deg, boundary_extrapolation, &

558 ppoly_r_coefs, ppoly_r_E)

559 integer, intent(in) :: n0 !< Number of cells on source grid

560 real, dimension(n0), intent(in) :: h0 !< Cell widths on source grid [H]

561 real, dimension(n0), intent(in) :: u0 !< Cell averages on source grid [A]

562 integer, intent(in) :: deg !< Degree of polynomial reconstruction

563 logical, intent(in) :: boundary_extrapolation !< Extrapolate at boundaries if true

564 real, dimension(n0,deg+1),intent(in) :: ppoly_r_coefs !< Coefficients of polynomial [A]

565 real, dimension(n0,2), intent(in) :: ppoly_r_E !< Edge value of polynomial [A]

566 ! Local variables

567 integer :: i0, n

568 real :: u_l, u_c, u_r ! Cell averages [A]

569 real :: u_min, u_max ! Cell extrema [A]

570 logical :: problem_detected

571

572 problem_detected = .false.

573 do i0 = 1, n0

574 u_l = u0(max(1,i0-1))

575 u_c = u0(i0)

576 u_r = u0(min(n0,i0+1))

577 if (i0 > 1 .or. .not. boundary_extrapolation) then

578 u_min = min(u_l, u_c)

579 u_max = max(u_l, u_c)

580 if (ppoly_r_e(i0,1) < u_min) then

581 write(0,'(a,I0,5(1x,a,1pe24.16))') 'Left edge undershoot at ',i0,'u(i0-1)=',u_l,'u(i0)=',u_c, &

582 'edge=',ppoly_r_e(i0,1),'err=',ppoly_r_e(i0,1)-u_min

583 problem_detected = .true.

584 endif

585 if (ppoly_r_e(i0,1) > u_max) then

586 write(0,'(a,I0,5(1x,a,1pe24.16))') 'Left edge overshoot at ',i0,'u(i0-1)=',u_l,'u(i0)=',u_c, &

587 'edge=',ppoly_r_e(i0,1),'err=',ppoly_r_e(i0,1)-u_max

588 problem_detected = .true.

589 endif

590 endif

591 if (i0 < n0 .or. .not. boundary_extrapolation) then

592 u_min = min(u_c, u_r)

593 u_max = max(u_c, u_r)

594 if (ppoly_r_e(i0,2) < u_min) then

595 write(0,'(a,I0,5(1x,a,1pe24.16))') 'Right edge undershoot at ',i0,'u(i0)=',u_c,'u(i0+1)=',u_r, &

596 'edge=',ppoly_r_e(i0,2),'err=',ppoly_r_e(i0,2)-u_min

597 problem_detected = .true.

598 endif

599 if (ppoly_r_e(i0,2) > u_max) then

600 write(0,'(a,I0,5(1x,a,1pe24.16))') 'Right edge overshoot at ',i0,'u(i0)=',u_c,'u(i0+1)=',u_r, &

601 'edge=',ppoly_r_e(i0,2),'err=',ppoly_r_e(i0,2)-u_max

602 problem_detected = .true.

603 endif

604 endif

605 if (i0 > 1) then

606 if ( (u_c-u_l)*(ppoly_r_e(i0,1)-ppoly_r_e(i0-1,2)) < 0.) then

607 write(0,'(a,I0,5(1x,a,1pe24.16))') 'Non-monotonic edges at',i0,'u(i0-1)=',u_l,'u(i0)=',u_c, &

608 'right edge=',ppoly_r_e(i0-1,2),'left edge=',ppoly_r_e(i0,1)

609 write(0,'(5(a,1pe24.16,1x))') 'u(i0)-u(i0-1)',u_c-u_l,'edge diff=',ppoly_r_e(i0,1)-ppoly_r_e(i0-1,2)

610 problem_detected = .true.

611 endif

612 endif

613 if (problem_detected) then

614 write(0,'(a,1p9e24.16)') 'Polynomial coeffs:',ppoly_r_coefs(i0,:)

615 write(0,'(3(a,1pe24.16,1x))') 'u_l=',u_l,'u_c=',u_c,'u_r=',u_r

616 write(0,'(a4,10a24)') 'i0','h0(i0)','u0(i0)','left edge','right edge','Polynomial coefficients'

617 do n = 1, n0

618 write(0,'(I0,1p10e24.16)') n,h0(n),u0(n),ppoly_r_e(n,1),ppoly_r_e(n,2),ppoly_r_coefs(n,:)

619 enddo

620 call mom_error(fatal, 'MOM_remapping, check_reconstructions_1d: '// &

621 'Edge values or polynomial coefficients were inconsistent!')

622 endif

623 enddo

624

625end subroutine check_reconstructions_1d

626

627!> Returns the intersection of source and targets grids along with and auxiliary lists or indices.

628!!

629!! For source grid with thicknesses h0(1:n0) and target grid with thicknesses h1(1:n1) the intersection

630!! or "subgrid" has thicknesses h_sub(1:n0+n1+1).

631!! h0 and h1 must have the same units. h_sub will return with the same units as h0 and h1.

632!!

633!! Notes on the algorithm:

634!! Internally, grids are defined by the interfaces (although we describe grids via thicknesses for accuracy).

635!! The intersection or union of two grids is thus defined by the super set of both lists of interfaces.

636!! Because both source and target grids can contain vanished cells, we do not eliminate repeated

637!! interfaces from the union.

638!! That is, the total number of interfaces of the sub-cells is equal to the total numer of interfaces of

639!! the source grid (n0+1) plus the total number of interfaces of the target grid (n1+1), i.e. n0+n1+2.

640!! Whenever target and source interfaces align, then the retention of identical interfaces leads to a

641!! vanished subcell.

642!! The remapping uses a common point of reference to the left (top) so there is always a vanished subcell

643!! at the left (top).

644!! If the total column thicknesses are the same, then the right (bottom) interfaces are also aligned and

645!! so the last subcell will also be vanished.

646subroutine intersect_src_tgt_grids( n0, h0, n1, h1, h_sub, h0_eff, &

647 isrc_start, isrc_end, isrc_max, itgt_start, itgt_end, isub_src )

648 integer, intent(in) :: n0 !< Number of cells in source grid

649 real, intent(in) :: h0(n0) !< Source grid widths (size n0) [H]

650 integer, intent(in) :: n1 !< Number of cells in target grid

651 real, intent(in) :: h1(n1) !< Target grid widths (size n1) [H]

652 real, intent(out) :: h_sub(n0+n1+1) !< Overlapping sub-cell thicknesses, h_sub [H]

653 real, intent(out) :: h0_eff(n0) !< Effective thickness of source cells [H]

654 integer, intent(out) :: isrc_start(n0) !< Index of first sub-cell within each source cell

655 integer, intent(out) :: isrc_end(n0) !< Index of last sub-cell within each source cell

656 integer, intent(out) :: isrc_max(n0) !< Index of thickest sub-cell within each source cell

657 integer, intent(out) :: itgt_start(n1) !< Index of first sub-cell within each target cell

658 integer, intent(out) :: itgt_end(n1) !< Index of last sub-cell within each target cell

659 integer, intent(out) :: isub_src(n0+n1+1) !< Index of source cell for each sub-cell

660 ! Local variables

661 integer :: i_sub ! Index of sub-cell

662 integer :: i0 ! Index into h0(1:n0), source column

663 integer :: i1 ! Index into h1(1:n1), target column

664 integer :: i_start0 ! Used to record which sub-cells map to source cells

665 integer :: i_start1 ! Used to record which sub-cells map to target cells

666 integer :: i_max ! Used to record which sub-cell is the largest contribution of a source cell

667 real :: dh_max ! Used to record which sub-cell is the largest contribution of a source cell [H]

668 real :: h0_supply, h1_supply ! The amount of width available for constructing sub-cells [H]

669 real :: dh ! The width of the sub-cell [H]

670 real :: dh0_eff ! Running sum of source cell thickness [H]

671 ! For error checking/debugging

672 logical :: src_has_volume !< True if h0 has not been consumed

673 logical :: tgt_has_volume !< True if h1 has not been consumed

674

675 ! Initialize algorithm

676 h0_supply = h0(1)

677 h1_supply = h1(1)

678 src_has_volume = .true.

679 tgt_has_volume = .true.

680 i0 = 1 ; i1 = 1

681 i_start0 = 1 ; i_start1 = 1

682 i_max = 1

683 dh_max = 0.

684 dh0_eff = 0.

685

686 ! First sub-cell is always vanished

687 h_sub(1) = 0.

688 isrc_start(1) = 1

689 isrc_end(1) = 1

690 isrc_max(1) = 1

691 isub_src(1) = 1

692

693 ! Loop over each sub-cell to calculate intersections with source and target grids

694 do i_sub = 2, n0+n1+1

695

696 ! This is the width of the sub-cell, determined by which ever column has the least

697 ! supply available to consume.

698 dh = min(h0_supply, h1_supply)

699

700 ! This is the running sum of the source cell thickness. After summing over each

701 ! sub-cell, the sum of sub-cell thickness might differ from the original source

702 ! cell thickness due to round off.

703 dh0_eff = dh0_eff + min(dh, h0_supply)

704

705 ! Record the source index (i0) that this sub-cell integral belongs to. This

706 ! is needed to index the reconstruction coefficients for the source cell

707 ! used in the integrals of the sub-cell width.

708 isub_src(i_sub) = i0

709 h_sub(i_sub) = dh

710

711 ! For recording the largest sub-cell within a source cell.

712 if (dh >= dh_max) then

713 i_max = i_sub

714 dh_max = dh

715 endif

716

717 ! Which ever column (source or target) has the least width left to consume determined

718 ! the width, dh, of sub-cell i_sub in the expression for dh above.

719 if (h0_supply <= h1_supply .and. src_has_volume) then

720 ! h0_supply is smaller than h1_supply) so we consume h0_supply and increment the

721 ! source cell index.

722 h1_supply = h1_supply - dh ! Although this is a difference the result will

723 ! be non-negative because of the conditional.

724 ! Record the sub-cell start/end index that span the source cell i0.

725 isrc_start(i0) = i_start0

726 isrc_end(i0) = i_sub

727 i_start0 = i_sub + 1

728 ! Record the sub-cell that is the largest fraction of the source cell.

729 isrc_max(i0) = i_max

730 i_max = i_sub + 1

731 dh_max = 0.

732 ! Record the source cell thickness found by summing the sub-cell thicknesses.

733 h0_eff(i0) = dh0_eff

734 ! Move the source index.

735 if (i0 < n0) then

736 i0 = i0 + 1

737 h0_supply = h0(i0)

738 dh0_eff = 0.

739 else

740 h0_supply = 0.

741 src_has_volume = .false.

742 endif

743 elseif (h0_supply >= h1_supply .and. tgt_has_volume) then

744 ! h1_supply is smaller than h0_supply) so we consume h1_supply and increment the

745 ! target cell index.

746 h0_supply = h0_supply - dh ! Although this is a difference the result will

747 ! be non-negative because of the conditional.

748 ! Record the sub-cell start/end index that span the target cell i1.

749 itgt_start(i1) = i_start1

750 itgt_end(i1) = i_sub

751 i_start1 = i_sub + 1

752 ! Move the target index.

753 if (i1 < n1) then

754 i1 = i1 + 1

755 h1_supply = h1(i1)

756 else

757 h1_supply = 0.

758 tgt_has_volume = .false.

759 endif

760 elseif (src_has_volume) then

761 ! We ran out of target volume but still have source cells to consume

762 h_sub(i_sub) = h0_supply

763 ! Record the sub-cell start/end index that span the source cell i0.

764 isrc_start(i0) = i_start0

765 isrc_end(i0) = i_sub

766 i_start0 = i_sub + 1

767 ! Record the sub-cell that is the largest fraction of the source cell.

768 isrc_max(i0) = i_max

769 i_max = i_sub + 1

770 dh_max = 0.

771 ! Record the source cell thickness found by summing the sub-cell thicknesses.

772 h0_eff(i0) = dh0_eff

773 if (i0 < n0) then

774 i0 = i0 + 1

775 h0_supply = h0(i0)

776 dh0_eff = 0.

777 else

778 h0_supply = 0.

779 src_has_volume = .false.

780 endif

781 elseif (tgt_has_volume) then

782 ! We ran out of source volume but still have target cells to consume

783 h_sub(i_sub) = h1_supply

784 ! Record the sub-cell start/end index that span the target cell i1.

785 itgt_start(i1) = i_start1

786 itgt_end(i1) = i_sub

787 i_start1 = i_sub + 1

788 ! Move the target index.

789 if (i1 < n1) then

790 i1 = i1 + 1

791 h1_supply = h1(i1)

792 else

793 h1_supply = 0.

794 tgt_has_volume = .false.

795 endif

796 else

797 stop 'intersect_src_tgt_grids: THIS SHOULD NEVER HAPPEN!'

798 endif

799

800 enddo

801

802end subroutine intersect_src_tgt_grids

803

804!> Adjust h_sub to ensure accurate conservation

805!!

806!! Loop over each source cell substituting the thickest sub-cell (within the source cell) with the

807!! residual of the source cell thickness minus the sum of other sub-cells

808!! aka a genius algorithm for accurate conservation when remapping from Robert Hallberg (\@Hallberg-NOAA).

809!subroutine adjust_h_sub( n0, h0, n1, isrc_start, isrc_end, isrc_max, h_sub )

810! integer, intent(in) :: n0 !< Number of cells in source grid

811! real, intent(in) :: h0(n0) !< Source grid widths (size n0) [H]

812! integer, intent(in) :: n1 !< Number of cells in target grid

813! integer, intent(in) :: isrc_start(n0) !< Index of first sub-cell within each source cell

814! integer, intent(in) :: isrc_end(n0) !< Index of last sub-cell within each source cell

815! integer, intent(in) :: isrc_max(n0) !< Index of thickest sub-cell within each source cell

816! real, intent(inout) :: h_sub(n0+n1+1) !< Overlapping sub-cell thicknesses, h_sub [H]

817! ! Local variables

818! integer :: i_sub ! Index of sub-cell

819! integer :: i0 ! Index into h0(1:n0), source column

820! integer :: i_max ! Used to record which sub-cell is the largest contribution of a source cell

821! real :: dh_max ! Used to record which sub-cell is the largest contribution of a source cell [H]

822! real :: dh ! The width of the sub-cell [H]

823! integer :: i0_last_thick_cell ! Last h0 cell with finite thickness

824!

825! i0_last_thick_cell = 0

826! do i0 = 1, n0

827! if (h0(i0)>0.) i0_last_thick_cell = i0

828! enddo

829!

830! do i0 = 1, i0_last_thick_cell

831! i_max = isrc_max(i0)

832! dh_max = h_sub(i_max)

833! if (dh_max > 0.) then

834! ! dh will be the sum of sub-cell thicknesses within the source cell except for the thickest sub-cell.

835! dh = 0.

836! do i_sub = isrc_start(i0), isrc_end(i0)

837! if (i_sub /= i_max) dh = dh + h_sub(i_sub)

838! enddo

839! h_sub(i_max) = h0(i0) - dh

840! endif

841! enddo

842!

843!end subroutine adjust_h_sub

844

845!> Remaps column of n0 values u0 on grid h0 to subgrid h_sub

846!!

847!! This includes an error for the scenario where the source grid is much thicker than

848!! the target grid and extrapolation is needed.

849subroutine remap_src_to_sub_grid_om4(n0, h0, u0, ppoly0_E, ppoly0_coefs, n1, h_sub, &

850 h0_eff, isrc_start, isrc_end, isrc_max, isub_src, &

851 method, force_bounds_in_subcell, u_sub, uh_sub, u02_err)

852 integer, intent(in) :: n0 !< Number of cells in source grid

853 real, intent(in) :: h0(n0) !< Source grid widths (size n0) [H]

854 real, intent(in) :: u0(n0) !< Source grid widths (size n0) [H]

855 real, intent(in) :: ppoly0_E(n0,2) !< Edge value of polynomial [A]

856 real, intent(in) :: ppoly0_coefs(:,:) !< Coefficients of polynomial [A]

857 integer, intent(in) :: n1 !< Number of cells in target grid

858 real, intent(in) :: h_sub(n0+n1+1) !< Overlapping sub-cell thicknesses, h_sub [H]

859 real, intent(in) :: h0_eff(n0) !< Effective thickness of source cells [H]

860 integer, intent(in) :: isrc_start(n0) !< Index of first sub-cell within each source cell

861 integer, intent(in) :: isrc_end(n0) !< Index of last sub-cell within each source cell

862 integer, intent(in) :: isrc_max(n0) !< Index of thickest sub-cell within each source cell

863 integer, intent(in) :: isub_src(n0+n1+1) !< Index of source cell for each sub-cell

864 integer, intent(in) :: method !< Remapping scheme to use

865 logical, intent(in) :: force_bounds_in_subcell !< Force sub-cell values to be bounded

866 real, intent(out) :: u_sub(n0+n1+1) !< Sub-cell cell averages (size n1) [A]

867 real, intent(out) :: uh_sub(n0+n1+1) !< Sub-cell cell integrals (size n1) [A H]

868 real, intent(out) :: u02_err !< Integrated reconstruction error estimates [A H]

869 ! Local variables

870 integer :: i_sub ! Index of sub-cell

871 integer :: i0 ! Index into h0(1:n0), source column

872 integer :: i_max ! Used to record which sub-cell is the largest contribution of a source cell

873 real :: dh_max ! Used to record which sub-cell is the largest contribution of a source cell [H]

874 real :: xa, xb ! Non-dimensional position within a source cell (0..1) [nondim]

875 real :: dh ! The width of the sub-cell [H]

876 real :: duh ! The total amount of accumulated stuff (u*h) [A H]

877 real :: dh0_eff ! Running sum of source cell thickness [H]

878 real :: u0_min(n0), u0_max(n0) !< Min/max of u0 for each source cell [A]

879 ! For error checking/debugging

880 logical, parameter :: adjust_thickest_subcell = .true. ! To fix round-off conservation issues

881 integer :: i0_last_thick_cell

882 real :: u_orig ! The original value of the reconstruction in a cell [A]

883

884 i0_last_thick_cell = 0

885 do i0 = 1, n0

886 u0_min(i0) = min(ppoly0_e(i0,1), ppoly0_e(i0,2))

887 u0_max(i0) = max(ppoly0_e(i0,1), ppoly0_e(i0,2))

888 if (h0(i0)>0.) i0_last_thick_cell = i0

889 enddo

890

891 ! Loop over each sub-cell to calculate average/integral values within each sub-cell.

892 ! Uses: h_sub, isub_src, h0_eff

893 ! Sets: u_sub, uh_sub

894 xa = 0.

895 dh0_eff = 0.

896 uh_sub(1) = 0.

897 u_sub(1) = ppoly0_e(1,1)

898 u02_err = 0.

899 do i_sub = 2, n0+n1

900

901 ! Sub-cell thickness from loop above

902 dh = h_sub(i_sub)

903

904 ! Source cell

905 i0 = isub_src(i_sub)

906

907 ! Evaluate average and integral for sub-cell i_sub.

908 ! Integral is over distance dh but expressed in terms of non-dimensional

909 ! positions with source cell from xa to xb (0 <= xa <= xb <= 1).

910 dh0_eff = dh0_eff + dh ! Cumulative thickness within the source cell

911 if (h0_eff(i0)>0.) then

912 xb = dh0_eff / h0_eff(i0) ! This expression yields xa <= xb <= 1.0

913 xb = min(1., xb) ! This is only needed when the total target column is wider than the source column

914 u_sub(i_sub) = average_value_ppoly( n0, u0, ppoly0_e, ppoly0_coefs, method, i0, xa, xb)

915 else ! Vanished cell

916 xb = 1.

917 u_sub(i_sub) = u0(i0)

918 endif

919 if (force_bounds_in_subcell) then

920 ! These next two lines should not be needed but when using PQM we found roundoff

921 ! can lead to overshoots. These lines sweep issues under the rug which need to be

922 ! properly .. later. -AJA

923 u_orig = u_sub(i_sub)

924 u_sub(i_sub) = max( u_sub(i_sub), u0_min(i0) )

925 u_sub(i_sub) = min( u_sub(i_sub), u0_max(i0) )

926 u02_err = u02_err + dh*abs( u_sub(i_sub) - u_orig )

927 endif

928 uh_sub(i_sub) = dh * u_sub(i_sub)

929

930 if (isub_src(i_sub+1) /= i0) then

931 ! If the next sub-cell is in a different source cell, reset the position counters

932 dh0_eff = 0.

933 xa = 0.

934 else

935 xa = xb ! Next integral will start at end of last

936 endif

937

938 enddo

939 u_sub(n0+n1+1) = ppoly0_e(n0,2) ! This value is only needed when total target column

940 uh_sub(n0+n1+1) = ppoly0_e(n0,2) * h_sub(n0+n1+1) ! is wider than the source column

941

942 if (adjust_thickest_subcell) then

943 ! Loop over each source cell substituting the integral/average for the thickest sub-cell (within

944 ! the source cell) with the residual of the source cell integral minus the other sub-cell integrals

945 ! aka a genius algorithm for accurate conservation when remapping from Robert Hallberg (\@Hallberg-NOAA).

946 ! Uses: i0_last_thick_cell, isrc_max, h_sub, isrc_start, isrc_end, uh_sub, u0, h0

947 ! Updates: uh_sub, u_sub

948 do i0 = 1, i0_last_thick_cell

949 i_max = isrc_max(i0)

950 dh_max = h_sub(i_max)

951 if (dh_max > 0.) then

952 ! duh will be the sum of sub-cell integrals within the source cell except for the thickest sub-cell.

953 duh = 0.

954 do i_sub = isrc_start(i0), isrc_end(i0)

955 if (i_sub /= i_max) duh = duh + uh_sub(i_sub)

956 enddo

957 uh_sub(i_max) = u0(i0)*h0(i0) - duh

958 u02_err = u02_err + max( abs(uh_sub(i_max)), abs(u0(i0)*h0(i0)), abs(duh) )

959 endif

960 enddo

961 endif

962

963end subroutine remap_src_to_sub_grid_om4

964

965!> Remaps column of n0 values u0 on grid h0 to subgrid h_sub

966subroutine remap_src_to_sub_grid(n0, h0, u0, ppoly0_E, ppoly0_coefs, n1, h_sub, &

967 isrc_start, isrc_end, isrc_max, isub_src, &

968 method, force_bounds_in_subcell, u_sub, uh_sub, u02_err)

969 integer, intent(in) :: n0 !< Number of cells in source grid

970 real, intent(in) :: h0(n0) !< Source grid widths (size n0) [H]

971 real, intent(in) :: u0(n0) !< Source grid widths (size n0) [H]

972 real, intent(in) :: ppoly0_E(n0,2) !< Edge value of polynomial [A]

973 real, intent(in) :: ppoly0_coefs(:,:) !< Coefficients of polynomial [A]

974 integer, intent(in) :: n1 !< Number of cells in target grid

975 real, intent(in) :: h_sub(n0+n1+1) !< Overlapping sub-cell thicknesses, h_sub [H]

976 integer, intent(in) :: isrc_start(n0) !< Index of first sub-cell within each source cell

977 integer, intent(in) :: isrc_end(n0) !< Index of last sub-cell within each source cell

978 integer, intent(in) :: isrc_max(n0) !< Index of thickest sub-cell within each source cell

979 integer, intent(in) :: isub_src(n0+n1+1) !< Index of source cell for each sub-cell

980 integer, intent(in) :: method !< Remapping scheme to use

981 logical, intent(in) :: force_bounds_in_subcell !< Force sub-cell values to be bounded

982 real, intent(out) :: u_sub(n0+n1+1) !< Sub-cell cell averages (size n1) [A]

983 real, intent(out) :: uh_sub(n0+n1+1) !< Sub-cell cell integrals (size n1) [A H]

984 real, intent(out) :: u02_err !< Integrated reconstruction error estimates [A H]

985 ! Local variables

986 integer :: i_sub ! Index of sub-cell

987 integer :: i0 ! Index into h0(1:n0), source column

988 integer :: i_max ! Used to record which sub-cell is the largest contribution of a source cell

989 real :: dh_max ! Used to record which sub-cell is the largest contribution of a source cell [H]

990 real :: xa, xb ! Non-dimensional position within a source cell (0..1) [nondim]

991 real :: dh ! The width of the sub-cell [H]

992 real :: duh ! The total amount of accumulated stuff (u*h) [A H]

993 real :: dh0_eff ! Running sum of source cell thickness [H]

994 real :: u0_min(n0), u0_max(n0) ! Min/max of u0 for each source cell [A]

995 ! For error checking/debugging

996 logical, parameter :: adjust_thickest_subcell = .true. ! To fix round-off conservation issues

997 integer :: i0_last_thick_cell

998 real :: u_orig ! The original value of the reconstruction in a cell [A]

999

1000 i0_last_thick_cell = 0

1001 do i0 = 1, n0

1002 u0_min(i0) = min(ppoly0_e(i0,1), ppoly0_e(i0,2))

1003 u0_max(i0) = max(ppoly0_e(i0,1), ppoly0_e(i0,2))

1004 if (h0(i0)>0.) i0_last_thick_cell = i0

1005 enddo

1006

1007 ! Loop over each sub-cell to calculate average/integral values within each sub-cell.

1008 ! Uses: h_sub, isub_src, h0_eff

1009 ! Sets: u_sub, uh_sub

1010 xa = 0.

1011 dh0_eff = 0.

1012 u02_err = 0.

1013 do i_sub = 1, n0+n1

1014

1015 ! Sub-cell thickness from loop above

1016 dh = h_sub(i_sub)

1017

1018 ! Source cell

1019 i0 = isub_src(i_sub)

1020

1021 ! Evaluate average and integral for sub-cell i_sub.

1022 ! Integral is over distance dh but expressed in terms of non-dimensional

1023 ! positions with source cell from xa to xb (0 <= xa <= xb <= 1).

1024 dh0_eff = dh0_eff + dh ! Cumulative thickness within the source cell

1025 if (h0(i0)>0.) then

1026 xb = dh0_eff / h0(i0) ! This expression yields xa <= xb <= 1.0

1027 xb = min(1., xb) ! This is only needed when the total target column is wider than the source column

1028 u_sub(i_sub) = average_value_ppoly( n0, u0, ppoly0_e, ppoly0_coefs, method, i0, xa, xb)

1029 else ! Vanished cell

1030 xb = 1.

1031 u_sub(i_sub) = u0(i0)

1032 endif

1033 if (force_bounds_in_subcell) then

1034 ! These next two lines should not be needed but when using PQM we found roundoff

1035 ! can lead to overshoots. These lines sweep issues under the rug which need to be

1036 ! properly .. later. -AJA

1037 u_orig = u_sub(i_sub)

1038 u_sub(i_sub) = max( u_sub(i_sub), u0_min(i0) )

1039 u_sub(i_sub) = min( u_sub(i_sub), u0_max(i0) )

1040 u02_err = u02_err + dh*abs( u_sub(i_sub) - u_orig )

1041 endif

1042 uh_sub(i_sub) = dh * u_sub(i_sub)

1043

1044 if (isub_src(i_sub+1) /= i0) then

1045 ! If the next sub-cell is in a different source cell, reset the position counters

1046 dh0_eff = 0.

1047 xa = 0.

1048 else

1049 xa = xb ! Next integral will start at end of last

1050 endif

1051

1052 enddo

1053 i_sub = n0+n1+1

1054 ! Sub-cell thickness from loop above

1055 dh = h_sub(i_sub)

1056 ! Source cell

1057 i0 = isub_src(i_sub)

1058

1059 ! Evaluate average and integral for sub-cell i_sub.

1060 ! Integral is over distance dh but expressed in terms of non-dimensional

1061 ! positions with source cell from xa to xb (0 <= xa <= xb <= 1).

1062 dh0_eff = dh0_eff + dh ! Cumulative thickness within the source cell

1063 if (h0(i0)>0.) then

1064 xb = dh0_eff / h0(i0) ! This expression yields xa <= xb <= 1.0

1065 xb = min(1., xb) ! This is only needed when the total target column is wider than the source column

1066 u_sub(i_sub) = average_value_ppoly( n0, u0, ppoly0_e, ppoly0_coefs, method, i0, xa, xb)

1067 else ! Vanished cell

1068 xb = 1.

1069 u_sub(i_sub) = u0(i0)

1070 endif

1071 if (force_bounds_in_subcell) then

1072 ! These next two lines should not be needed but when using PQM we found roundoff

1073 ! can lead to overshoots. These lines sweep issues under the rug which need to be

1074 ! properly .. later. -AJA

1075 u_orig = u_sub(i_sub)

1076 u_sub(i_sub) = max( u_sub(i_sub), u0_min(i0) )

1077 u_sub(i_sub) = min( u_sub(i_sub), u0_max(i0) )

1078 u02_err = u02_err + dh*abs( u_sub(i_sub) - u_orig )

1079 endif

1080 uh_sub(i_sub) = dh * u_sub(i_sub)

1081

1082 if (adjust_thickest_subcell) then

1083 ! Loop over each source cell substituting the integral/average for the thickest sub-cell (within

1084 ! the source cell) with the residual of the source cell integral minus the other sub-cell integrals

1085 ! aka a genius algorithm for accurate conservation when remapping from Robert Hallberg (\@Hallberg-NOAA).

1086 ! Uses: i0_last_thick_cell, isrc_max, h_sub, isrc_start, isrc_end, uh_sub, u0, h0

1087 ! Updates: uh_sub

1088 do i0 = 1, i0_last_thick_cell

1089 i_max = isrc_max(i0)

1090 dh_max = h_sub(i_max)

1091 if (dh_max > 0.) then

1092 ! duh will be the sum of sub-cell integrals within the source cell except for the thickest sub-cell.

1093 duh = 0.

1094 do i_sub = isrc_start(i0), isrc_end(i0)

1095 if (i_sub /= i_max) duh = duh + uh_sub(i_sub)

1096 enddo

1097 uh_sub(i_max) = u0(i0)*h0(i0) - duh

1098 u02_err = u02_err + max( abs(uh_sub(i_max)), abs(u0(i0)*h0(i0)), abs(duh) )

1099 endif

1100 enddo

1101 endif

1102

1103end subroutine remap_src_to_sub_grid

1104

1105!> Remaps column of n0+n1+1 values usub on sub-grid h_sub to targets on grid h1

1106!! using the OM4-era algorithm

1107subroutine remap_sub_to_tgt_grid_om4(n0, n1, h1, h_sub, u_sub, uh_sub, &

1108 itgt_start, itgt_end, force_bounds_in_target, u1, uh_err)

1109 integer, intent(in) :: n0 !< Number of cells in source grid

1110 integer, intent(in) :: n1 !< Number of cells in target grid

1111 real, intent(in) :: h1(n1) !< Target grid widths (size n1) [H]

1112 real, intent(in) :: h_sub(n0+n1+1) !< Overlapping sub-cell thicknesses, h_sub [H]

1113 real, intent(in) :: u_sub(n0+n1+1) !< Sub-cell cell averages (size n1) [A]

1114 real, intent(in) :: uh_sub(n0+n1+1) !< Sub-cell cell integrals (size n1) [A H]

1115 integer, intent(in) :: itgt_start(n1) !< Index of first sub-cell within each target cell

1116 integer, intent(in) :: itgt_end(n1) !< Index of last sub-cell within each target cell

1117 logical, intent(in) :: force_bounds_in_target !< Force sub-cell values to be bounded

1118 real, intent(out) :: u1(n1) !< Target cell averages (size n1) [A]

1119 real, intent(out) :: uh_err !< Estimate of bound on error in sum of u*h [A H]

1120 ! Local variables

1121 integer :: i1 ! tgt loop index

1122 integer :: i_sub ! index to sub-layer

1123 real :: dh ! The width of the sub-cell [H]

1124 real :: duh ! The total amount of accumulated stuff (u*h) [A H]

1125 real :: u1min, u1max ! Minimum and maximum values of reconstructions [A]

1126 real :: u_orig ! The original value of the reconstruction in a cell prior to bounding [A]

1127

1128 u1min = 0. ! Not necessary, but avoids an overzealous compiler ...

1129 u1max = 0. ! ... warning about uninitialized variables

1130

1131 ! Loop over each target cell summing the integrals from sub-cells within the target cell.

1132 ! Uses: itgt_start, itgt_end, h_sub, uh_sub, u_sub

1133 ! Sets: u1, uh_err

1134 uh_err = 0.

1135 do i1 = 1, n1

1136 if (h1(i1) > 0.) then

1137 duh = 0. ; dh = 0.

1138 i_sub = itgt_start(i1)

1139 if (force_bounds_in_target) then

1140 u1min = u_sub(i_sub)

1141 u1max = u_sub(i_sub)

1142 endif

1143 do i_sub = itgt_start(i1), itgt_end(i1)

1144 if (force_bounds_in_target) then

1145 u1min = min(u1min, u_sub(i_sub))

1146 u1max = max(u1max, u_sub(i_sub))

1147 endif

1148 dh = dh + h_sub(i_sub)

1149 duh = duh + uh_sub(i_sub)

1150 ! This accumulates the contribution to the error bound for the sum of u*h

1151 uh_err = uh_err + max(abs(duh),abs(uh_sub(i_sub)))*epsilon(duh)

1152 enddo

1153 u1(i1) = duh / dh

1154 ! This is the contribution from the division to the error bound for the sum of u*h

1155 uh_err = uh_err + abs(duh)*epsilon(duh)

1156 if (force_bounds_in_target) then

1157 u_orig = u1(i1)

1158 u1(i1) = max(u1min, min(u1max, u1(i1)))

1159 ! Adjusting to be bounded contributes to the error for the sum of u*h

1160 uh_err = uh_err + dh*abs( u1(i1)-u_orig )

1161 endif

1162 else

1163 u1(i1) = u_sub(itgt_start(i1))

1164 endif

1165 enddo

1166

1167end subroutine remap_sub_to_tgt_grid_om4

1168

1169!> Remaps column of n0+n1+1 values usub on sub-grid h_sub to targets on grid h1

1170subroutine remap_sub_to_tgt_grid(n0, n1, h1, h_sub, u_sub, uh_sub, &

1171 itgt_start, itgt_end, force_bounds_in_target, &

1172 better_force_bounds_in_target, offset_summation, u1, uh_err)

1173 integer, intent(in) :: n0 !< Number of cells in source grid

1174 integer, intent(in) :: n1 !< Number of cells in target grid

1175 real, intent(in) :: h1(n1) !< Target grid widths (size n1) [H]

1176 real, intent(in) :: h_sub(n0+n1+1) !< Overlapping sub-cell thicknesses, h_sub [H]

1177 real, intent(in) :: u_sub(n0+n1+1) !< Sub-cell cell averages (size n1) [A]

1178 real, intent(in) :: uh_sub(n0+n1+1) !< Sub-cell cell integrals (size n1) [A H]

1179 integer, intent(in) :: itgt_start(n1) !< Index of first sub-cell within each target cell

1180 integer, intent(in) :: itgt_end(n1) !< Index of last sub-cell within each target cell

1181 logical, intent(in) :: force_bounds_in_target !< Force sub-cell values to be bounded

1182 logical, intent(in) :: better_force_bounds_in_target !< Force sub-cell values to be bounded

1183 logical, intent(in) :: offset_summation !< Offset values in summation for accuracy

1184 real, intent(out) :: u1(n1) !< Target cell averages (size n1) [A]

1185 real, intent(out) :: uh_err !< Estimate of bound on error in sum of u*h [A H]

1186 ! Local variables

1187 integer :: i1 ! tgt loop index

1188 integer :: i_sub ! index to sub-layer

1189 real :: dh ! The width of the sub-cell [H]

1190 real :: duh ! The total amount of accumulated stuff (u*h) [A H]

1191 real :: u1min, u1max ! Minimum and maximum values of reconstructions [A]

1192 real :: u_orig ! The original value of the reconstruction in a cell prior to bounding [A]

1193 real :: u_ref ! A value to offest the summation to gain accuracy [A]

1194 real :: h_max ! Thickest cell encountered [H]

1195

1196 u1min = 0. ! Not necessary, but avoids an overzealous compiler ...

1197 u1max = 0. ! ... warning about uninitialized variables

1198 u_ref = 0. ! An offset of 0. should do no harm

1199 h_max = 0.

1200

1201 ! Loop over each target cell summing the integrals from sub-cells within the target cell.

1202 ! Uses: itgt_start, itgt_end, h_sub, uh_sub, u_sub

1203 ! Sets: u1, uh_err

1204 uh_err = 0.

1205 do i1 = 1, n1

1206 if (h1(i1) > 0.) then

1207 duh = 0. ; dh = 0.

1208 i_sub = itgt_start(i1)

1209 if (force_bounds_in_target) then

1210 u1min = u_sub(i_sub)

1211 u1max = u_sub(i_sub)

1212 endif

1213 if (offset_summation) then

1214 u_ref = 0. ! An offset of 0. should do no harm

1215 h_max = 0.

1216 do i_sub = itgt_start(i1), itgt_end(i1)

1217 if (h_sub(i_sub) > h_max) then

1218 u_ref = u_sub(i_sub)

1219 h_max = h_sub(i_sub)

1220 endif

1221 enddo

1222 endif

1223 do i_sub = itgt_start(i1), itgt_end(i1)

1224 if (force_bounds_in_target .or. better_force_bounds_in_target .and. h_sub(i_sub)>0.) then

1225 u1min = min(u1min, u_sub(i_sub))

1226 u1max = max(u1max, u_sub(i_sub))

1227 endif

1228 dh = dh + h_sub(i_sub)

1229 ! Ideally u_ref would be already be substracted in uh_sub

1230 duh = duh + ( uh_sub(i_sub) - h_sub(i_sub) * u_ref )

1231 ! This accumulates the contribution to the error bound for the sum of u*h

1232 uh_err = uh_err + max(abs(duh),abs(uh_sub(i_sub)))*epsilon(duh)

1233 enddo

1234 u1(i1) = duh / dh + u_ref

1235 ! This is the contribution from the division to the error bound for the sum of u*h

1236 uh_err = uh_err + abs(duh)*epsilon(duh)

1237 if (force_bounds_in_target) then

1238 u_orig = u1(i1)

1239 u1(i1) = max(u1min, min(u1max, u1(i1)))

1240 ! Adjusting to be bounded contributes to the error for the sum of u*h

1241 uh_err = uh_err + dh*abs( u1(i1)-u_orig )

1242 endif

1243 else

1244 u1(i1) = u_sub(itgt_start(i1))

1245 endif

1246 enddo

1247

1248end subroutine remap_sub_to_tgt_grid

1249

1250!> Linearly interpolate interface data, u_src, from grid h_src to a grid h_dest

1251subroutine interpolate_column(nsrc, h_src, u_src, ndest, h_dest, u_dest, mask_edges)

1252 integer, intent(in) :: nsrc !< Number of source cells

1253 real, dimension(nsrc), intent(in) :: h_src !< Thickness of source cells [H]

1254 real, dimension(nsrc+1), intent(in) :: u_src !< Values at source cell interfaces [A]

1255 integer, intent(in) :: ndest !< Number of destination cells

1256 real, dimension(ndest), intent(in) :: h_dest !< Thickness of destination cells [H]

1257 real, dimension(ndest+1), intent(inout) :: u_dest !< Interpolated value at destination cell interfaces [A]

1258 logical, intent(in) :: mask_edges !< If true, mask the values outside of massless

1259 !! layers at the top and bottom of the column.

1260

1261 ! Local variables

1262 real :: x_dest ! Relative position of target interface [H]

1263 real :: dh ! Source cell thickness [H]

1264 real :: frac_pos(ndest+1) ! Fractional position of the destination interface

1265 ! within the source layer [nondim], 0 <= frac_pos <= 1.

1266 integer :: k_src(ndest+1) ! Source grid layer index of destination interface, 1 <= k_src <= ndest.

1267 integer :: ks, k_dest ! Index of cell in src and dest columns

1268

1269 ! The following forces the "do while" loop to do one cycle that will set u1, u2, dh.

1270 ks = 0

1271 dh = 0.

1272 x_dest = 0.

1273

1274 ! Find the layer index and fractional position of the interfaces of the target

1275 ! grid on the source grid.

1276 do k_dest=1,ndest+1

1277 do while (dh<=x_dest .and. ks<nsrc)

1278 ! Move positions pointers forward until the interval 0 .. dh spans x_dest.

1279 x_dest = x_dest - dh

1280 ks = ks + 1

1281 dh = h_src(ks) ! Thickness of layer ks

1282 enddo

1283 k_src(k_dest) = ks

1284

1285 if (dh>0.) then

1286 frac_pos(k_dest) = max(0., min(1., x_dest / dh)) ! Weight of u2

1287 else ! For a vanished source layer we need to do something reasonable...

1288 frac_pos(k_dest) = 0.5

1289 endif

1290

1291 if (k_dest <= ndest) then

1292 x_dest = x_dest + h_dest(k_dest) ! Position of interface k_dest+1

1293 endif

1294 enddo

1295

1296 do k_dest=1,ndest+1

1297 ! Linear interpolation between surrounding edge values.

1298 ks = k_src(k_dest)

1299 u_dest(k_dest) = (1.0 - frac_pos(k_dest)) * u_src(ks) + frac_pos(k_dest) * u_src(ks+1)

1300 enddo

1301

1302 if (mask_edges) then

1303 ! Mask vanished layers at the surface which would be under an ice-shelf.

1304 ! When the layer k_dest is vanished and all layers above are also vanished,

1305 ! the k_dest interface value should be missing.

1306 do k_dest=1,ndest

1307 if (h_dest(k_dest) > 0.) exit

1308 u_dest(k_dest) = 0.0

1309 enddo

1310

1311 ! Mask interfaces below vanished layers at the bottom

1312 do k_dest=ndest,1,-1

1313 if (h_dest(k_dest) > 0.) exit

1314 u_dest(k_dest+1) = 0.0

1315 enddo

1316 endif

1317

1318end subroutine interpolate_column

1319

1320!> Conservatively calculate integrated data, uh_dest, on grid h_dest, from layer-integrated data, uh_src, on grid h_src

1321subroutine reintegrate_column(nsrc, h_src, uh_src, ndest, h_dest, uh_dest)

1322 integer, intent(in) :: nsrc !< Number of source cells

1323 real, dimension(nsrc), intent(in) :: h_src !< Thickness of source cells [H]

1324 real, dimension(nsrc), intent(in) :: uh_src !< Values at source cell interfaces [A H]

1325 integer, intent(in) :: ndest !< Number of destination cells

1326 real, dimension(ndest), intent(in) :: h_dest !< Thickness of destination cells [H]

1327 real, dimension(ndest), intent(inout) :: uh_dest !< Interpolated value at destination cell interfaces [A H]

1328

1329 ! Local variables

1330 real :: h_src_rem, h_dest_rem, dh ! Incremental thicknesses [H]

1331 real :: uh_src_rem, duh ! Incremental amounts of stuff [A H]

1332 integer :: k_src, k_dest ! Index of cell in src and dest columns

1333 logical :: src_ran_out

1334

1335 uh_dest(:) = 0.0

1336

1337 k_src = 0

1338 k_dest = 0

1339 h_dest_rem = 0.

1340 h_src_rem = 0.

1341 uh_src_rem = 0.

1342 src_ran_out = .false.

1343

1344 do while(.true.)

1345 if (h_src_rem==0. .and. k_src<nsrc) then

1346 ! Supply is empty so move to the next source cell

1347 k_src = k_src + 1

1348 h_src_rem = h_src(k_src)

1349 uh_src_rem = uh_src(k_src)

1350 if (h_src_rem==0.) cycle

1351 endif

1352 if (h_dest_rem==0. .and. k_dest<ndest) then

1353 ! Sink has no capacity so move to the next destination cell

1354 k_dest = k_dest + 1

1355 h_dest_rem = h_dest(k_dest)

1356 uh_dest(k_dest) = 0.

1357 if (h_dest_rem==0.) cycle

1358 endif

1359 if (k_src==nsrc .and. h_src_rem==0.) then

1360 if (src_ran_out) exit ! This is the second time implying there is no more src

1361 src_ran_out = .true.

1362 cycle

1363 endif

1364 duh = 0.

1365 if (h_src_rem<h_dest_rem) then

1366 ! The source cell is fully within the destination cell

1367 dh = h_src_rem

1368 if (dh>0.) duh = uh_src_rem

1369 h_src_rem = 0.

1370 uh_src_rem = 0.

1371 h_dest_rem = max(0., h_dest_rem - dh)

1372 elseif (h_src_rem>h_dest_rem) then

1373 ! Only part of the source cell can be used up

1374 dh = h_dest_rem

1375 duh = (dh / h_src_rem) * uh_src_rem

1376 h_src_rem = max(0., h_src_rem - dh)

1377 uh_src_rem = uh_src_rem - duh

1378 h_dest_rem = 0.

1379 else ! h_src_rem==h_dest_rem

1380 ! The source cell exactly fits the destination cell

1381 duh = uh_src_rem

1382 h_src_rem = 0.

1383 uh_src_rem = 0.

1384 h_dest_rem = 0.

1385 endif

1386 uh_dest(k_dest) = uh_dest(k_dest) + duh

1387 if (k_dest==ndest .and. (k_src==nsrc .or. h_dest_rem==0.)) exit

1388 enddo

1389

1390end subroutine reintegrate_column

1391

1392!> Returns the average value of a reconstruction within a single source cell, i0,

1393!! between the non-dimensional positions xa and xb (xa<=xb) with dimensional

1394!! separation dh.

1395real function average_value_ppoly( n0, u0, ppoly0_E, ppoly0_coefs, method, i0, xa, xb)

1396 integer, intent(in) :: n0 !< Number of cells in source grid

1397 real, intent(in) :: u0(n0) !< Cell means [A]

1398 real, intent(in) :: ppoly0_e(n0,2) !< Edge value of polynomial [A]

1399 real, intent(in) :: ppoly0_coefs(:,:) !< Coefficients of polynomial [A]

1400 integer, intent(in) :: method !< Remapping scheme to use

1401 integer, intent(in) :: i0 !< Source cell index

1402 real, intent(in) :: xa !< Non-dimensional start position within source cell [nondim]

1403 real, intent(in) :: xb !< Non-dimensional end position within source cell [nondim]

1404 ! Local variables

1405 real :: u_ave ! The average value of the polynomial over the specified range [A]

1406 real :: xapxb ! A sum of fracional positions [nondim]

1407 real :: mx, ya, yb, my ! Various fractional positions [nondim]

1408 real :: xa_2, xb_2 ! Squared fractional positions [nondim]

1409 real :: xa2pxb2, xa2b2ab, ya2b2ab ! Sums of squared fractional positions [nondim]

1410 real :: a_l, a_r, u_c, a_c ! Values of the polynomial at various locations [A]

1411 real, parameter :: r_3 = 1.0/3.0 ! Used in evaluation of integrated polynomials [nondim]

1412

1413 u_ave = 0. ! Avoids warnings about "potentially unset values"; u_ave is always calculated for legitimate schemes

1414 if (xb > xa) then

1415 select case ( method )

1416 case ( integration_pcm )

1417 u_ave = u0(i0)

1418 case ( integration_plm )

1419 u_ave = ( &

1420 ppoly0_coefs(i0,1) &

1421 + ppoly0_coefs(i0,2) * 0.5 * ( xb + xa ) )

1422 case ( integration_ppm )

1423 mx = 0.5 * ( xa + xb )

1424 a_l = ppoly0_e(i0, 1)

1425 a_r = ppoly0_e(i0, 2)

1426 u_c = u0(i0)

1427 a_c = 0.5 * ( ( u_c - a_l ) + ( u_c - a_r ) ) ! a_6 / 6

1428 if (mx<0.5) then

1429 ! This integration of the PPM reconstruction is expressed in distances from the left edge

1430 xa2b2ab = (xa*xa+xb*xb)+xa*xb

1431 u_ave = a_l + ( ( a_r - a_l ) * mx &

1432 + a_c * ( 3. * ( xb + xa ) - 2.*xa2b2ab ) )

1433 else

1434 ! This integration of the PPM reconstruction is expressed in distances from the right edge

1435 ya = 1. - xa

1436 yb = 1. - xb

1437 my = 0.5 * ( ya + yb )

1438 ya2b2ab = (ya*ya+yb*yb)+ya*yb

1439 u_ave = a_r + ( ( a_l - a_r ) * my &

1440 + a_c * ( 3. * ( yb + ya ) - 2.*ya2b2ab ) )

1441 endif

1442 case ( integration_pqm )

1443 xa_2 = xa*xa

1444 xb_2 = xb*xb

1445 xa2pxb2 = xa_2 + xb_2

1446 xapxb = xa + xb

1447 u_ave = ( &

1448 ppoly0_coefs(i0,1) &

1449 + ( ppoly0_coefs(i0,2) * 0.5 * ( xapxb ) &

1450 + ( ppoly0_coefs(i0,3) * r_3 * ( xa2pxb2 + xa*xb ) &

1451 + ( ppoly0_coefs(i0,4) * 0.25* ( xa2pxb2 * xapxb ) &

1452 + ppoly0_coefs(i0,5) * 0.2 * ( ( xb*xb_2 + xa*xa_2 ) * xapxb + xa_2*xb_2 ) ) ) ) )

1453 case default

1454 call mom_error( fatal,'The selected integration method is invalid' )

1455 end select

1456 else ! dh == 0.

1457 select case ( method )

1458 case ( integration_pcm )

1459 u_ave = ppoly0_coefs(i0,1)

1460 case ( integration_plm )

1461 !u_ave = ppoly0_coefs(i0,1) &

1462 ! + xa * ppoly0_coefs(i0,2)

1463 a_l = ppoly0_e(i0, 1)

1464 a_r = ppoly0_e(i0, 2)

1465 ya = 1. - xa

1466 if (xa < 0.5) then

1467 u_ave = a_l + xa * ( a_r - a_l )

1468 else

1469 u_ave = a_r + ya * ( a_l - a_r )

1470 endif

1471 case ( integration_ppm )

1472 !u_ave = ppoly0_coefs(i0,1) &

1473 ! + xa * ( ppoly0_coefs(i0,2) &

1474 ! + xa * ppoly0_coefs(i0,3) )

1475 a_l = ppoly0_e(i0, 1)

1476 a_r = ppoly0_e(i0, 2)

1477 u_c = u0(i0)

1478 a_c = 3. * ( ( u_c - a_l ) + ( u_c - a_r ) ) ! a_6

1479 ya = 1. - xa

1480 if (xa < 0.5) then

1481 u_ave = a_l + xa * ( ( a_r - a_l ) + a_c * ya )

1482 else

1483 u_ave = a_r + ya * ( ( a_l - a_r ) + a_c * xa )

1484 endif

1485 case ( integration_pqm )

1486 u_ave = ppoly0_coefs(i0,1) &

1487 + xa * ( ppoly0_coefs(i0,2) &

1488 + xa * ( ppoly0_coefs(i0,3) &

1489 + xa * ( ppoly0_coefs(i0,4) &

1490 + xa * ppoly0_coefs(i0,5) ) ) )

1491 case default

1492 u_ave = 0.

1493 call mom_error( fatal,'The selected integration method is invalid' )

1494 end select

1495 endif

1496 average_value_ppoly = u_ave

1497

1498end function average_value_ppoly

1499

1500!> This subroutine checks for sufficient consistence in the extrema and total amounts on the old

1501!! and new grids.

1502subroutine check_remapped_values(n0, h0, u0, ppoly_r_E, deg, ppoly_r_coefs, &

1503 n1, h1, u1, iMethod, uh_err, caller)

1504 integer, intent(in) :: n0 !< Number of cells on source grid

1505 real, dimension(n0), intent(in) :: h0 !< Cell widths on source grid [H]

1506 real, dimension(n0), intent(in) :: u0 !< Cell averages on source grid [A]

1507 real, dimension(n0,2), intent(in) :: ppoly_r_E !< Edge values of polynomial fits [A]

1508 integer, intent(in) :: deg !< Degree of the piecewise polynomial reconstrution

1509 real, dimension(n0,deg+1), intent(in) :: ppoly_r_coefs !< Coefficients of the piecewise

1510 !! polynomial reconstructions [A]

1511 integer, intent(in) :: n1 !< Number of cells on target grid

1512 real, dimension(n1), intent(in) :: h1 !< Cell widths on target grid [H]

1513 real, dimension(n1), intent(in) :: u1 !< Cell averages on target grid [A]

1514 integer, intent(in) :: iMethod !< An integer indicating the integration method used

1515 real, intent(in) :: uh_err !< A bound on the error in the sum of u*h as

1516 !! estimated by the remapping code [H A]

1517 character(len=*), intent(in) :: caller !< The name of the calling routine.

1518

1519 ! Local variables

1520 real :: h0tot, h0err ! Sum of source cell widths and round-off error in this sum [H]

1521 real :: h1tot, h1err ! Sum of target cell widths and round-off error in this sum [H]

1522 real :: u0tot, u0err ! Integrated values on the source grid and round-off error in this sum [H A]

1523 real :: u1tot, u1err ! Integrated values on the target grid and round-off error in this sum [H A]

1524 real :: u0min, u0max, u1min, u1max ! Extrema of values on the two grids [A]

1525 integer :: k

1526

1527 ! Check errors and bounds

1528 call measure_input_bounds( n0, h0, u0, ppoly_r_e, h0tot, h0err, u0tot, u0err, u0min, u0max )

1529 call measure_output_bounds( n1, h1, u1, h1tot, h1err, u1tot, u1err, u1min, u1max )

1530

1531 if (imethod<5) return ! We except PQM until we've debugged it

1532

1533 if ( (abs(u1tot-u0tot)>(u0err+u1err)+uh_err .and. abs(h1tot-h0tot)<h0err+h1err) &

1534 .or. (u1min<u0min .or. u1max>u0max) ) then

1535 write(0,*) 'iMethod = ',imethod

1536 write(0,*) 'H: h0tot=',h0tot,'h1tot=',h1tot,'dh=',h1tot-h0tot,'h0err=',h0err,'h1err=',h1err

1537 if (abs(h1tot-h0tot)>h0err+h1err) &

1538 write(0,*) 'H non-conservation difference=',h1tot-h0tot,'allowed err=',h0err+h1err,' <-----!'

1539 write(0,*) 'UH: u0tot=',u0tot,'u1tot=',u1tot,'duh=',u1tot-u0tot,'u0err=',u0err,'u1err=',u1err,'uh_err=',uh_err

1540 if (abs(u1tot-u0tot)>(u0err+u1err)+uh_err) &

1541 write(0,*) 'U non-conservation difference=',u1tot-u0tot,'allowed err=',u0err+u1err+uh_err,' <-----!'

1542 write(0,*) 'U: u0min=',u0min,'u1min=',u1min

1543 if (u1min<u0min) write(0,*) 'U minimum overshoot=',u1min-u0min,' <-----!'

1544 write(0,*) 'U: u0max=',u0max,'u1max=',u1max

1545 if (u1max>u0max) write(0,*) 'U maximum overshoot=',u1max-u0max,' <-----!'

1546 write(0,'(a3,6a24)') 'k','h0','left edge','u0','right edge','h1','u1'

1547 do k = 1, max(n0,n1)

1548 if (k<=min(n0,n1)) then

1549 write(0,'(i3,1p6e24.16)') k,h0(k),ppoly_r_e(k,1),u0(k),ppoly_r_e(k,2),h1(k),u1(k)

1550 elseif (k>n0) then

1551 write(0,'(i3,96x,1p2e24.16)') k,h1(k),u1(k)

1552 else

1553 write(0,'(i3,1p4e24.16)') k,h0(k),ppoly_r_e(k,1),u0(k),ppoly_r_e(k,2)

1554 endif

1555 enddo

1556 write(0,'(a3,2a24)') 'k','u0','Polynomial coefficients'

1557 do k = 1, n0

1558 write(0,'(i3,1p6e24.16)') k,u0(k),ppoly_r_coefs(k,:)

1559 enddo

1560 call mom_error( fatal, 'MOM_remapping, '//trim(caller)//': '//&

1561 'Remapping result is inconsistent!' )

1562 endif

1563

1564end subroutine check_remapped_values

1565

1566!> Measure totals and bounds on source grid

1567subroutine measure_input_bounds( n0, h0, u0, edge_values, h0tot, h0err, u0tot, u0err, u0min, u0max )

1568 integer, intent(in) :: n0 !< Number of cells on source grid

1569 real, dimension(n0), intent(in) :: h0 !< Cell widths on source grid [H]

1570 real, dimension(n0), intent(in) :: u0 !< Cell averages on source grid [A]

1571 real, dimension(n0,2), intent(in) :: edge_values !< Cell edge values on source grid [A]

1572 real, intent(out) :: h0tot !< Sum of cell widths [H]

1573 real, intent(out) :: h0err !< Magnitude of round-off error in h0tot [H]

1574 real, intent(out) :: u0tot !< Sum of cell widths times values [H A]

1575 real, intent(out) :: u0err !< Magnitude of round-off error in u0tot [H A]

1576 real, intent(out) :: u0min !< Minimum value in reconstructions of u0 [A]

1577 real, intent(out) :: u0max !< Maximum value in reconstructions of u0 [A]

1578 ! Local variables

1579 real :: eps ! The smallest representable fraction of a number [nondim]

1580 integer :: k

1581

1582 eps = epsilon(h0(1))

1583 h0tot = h0(1)

1584 h0err = 0.

1585 u0tot = h0(1) * u0(1)

1586 u0err = 0.

1587 u0min = min( edge_values(1,1), edge_values(1,2) )

1588 u0max = max( edge_values(1,1), edge_values(1,2) )

1589 do k = 2, n0

1590 h0tot = h0tot + h0(k)

1591 h0err = h0err + eps * max(h0tot, h0(k))

1592 u0tot = u0tot + h0(k) * u0(k)

1593 u0err = u0err + eps * max(abs(u0tot), abs(h0(k) * u0(k)))

1594 u0min = min( u0min, edge_values(k,1), edge_values(k,2) )

1595 u0max = max( u0max, edge_values(k,1), edge_values(k,2) )

1596 enddo

1597

1598end subroutine measure_input_bounds

1599

1600!> Measure totals and bounds on destination grid

1601subroutine measure_output_bounds( n1, h1, u1, h1tot, h1err, u1tot, u1err, u1min, u1max )

1602 integer, intent(in) :: n1 !< Number of cells on destination grid

1603 real, dimension(n1), intent(in) :: h1 !< Cell widths on destination grid [H]

1604 real, dimension(n1), intent(in) :: u1 !< Cell averages on destination grid [A]

1605 real, intent(out) :: h1tot !< Sum of cell widths [H]

1606 real, intent(out) :: h1err !< Magnitude of round-off error in h1tot [H]

1607 real, intent(out) :: u1tot !< Sum of cell widths times values [H A]

1608 real, intent(out) :: u1err !< Magnitude of round-off error in u1tot [H A]

1609 real, intent(out) :: u1min !< Minimum value in reconstructions of u1 [A]

1610 real, intent(out) :: u1max !< Maximum value in reconstructions of u1 [A]

1611 ! Local variables

1612 real :: eps ! The smallest representable fraction of a number [nondim]

1613 integer :: k

1614

1615 eps = epsilon(h1(1))

1616 h1tot = h1(1)

1617 h1err = 0.

1618 u1tot = h1(1) * u1(1)

1619 u1err = 0.

1620 u1min = u1(1)

1621 u1max = u1(1)

1622 do k = 2, n1

1623 h1tot = h1tot + h1(k)

1624 h1err = h1err + eps * max(h1tot, h1(k))

1625 u1tot = u1tot + h1(k) * u1(k)

1626 u1err = u1err + eps * max(abs(u1tot), abs(h1(k) * u1(k)))

1627 u1min = min(u1min, u1(k))

1628 u1max = max(u1max, u1(k))

1629 enddo

1630

1631end subroutine measure_output_bounds

1632

1633!> Calculates the change in interface positions based on h1 and h2

1634subroutine dzfromh1h2( n1, h1, n2, h2, dx )

1635 integer, intent(in) :: n1 !< Number of cells on source grid

1636 real, dimension(:), intent(in) :: h1 !< Cell widths of source grid (size n1) [H]

1637 integer, intent(in) :: n2 !< Number of cells on target grid

1638 real, dimension(:), intent(in) :: h2 !< Cell widths of target grid (size n2) [H]

1639 real, dimension(:), intent(out) :: dx !< Change in interface position (size n2+1) [H]

1640 ! Local variables

1641 integer :: k

1642 real :: x1, x2 ! Interface positions [H]

1643

1644 x1 = 0.

1645 x2 = 0.

1646 dx(1) = 0.

1647 do k = 1, max(n1,n2)

1648 if (k <= n1) x1 = x1 + h1(k) ! Interface k+1, right of source cell k

1649 if (k <= n2) then

1650 x2 = x2 + h2(k) ! Interface k+1, right of target cell k

1651 dx(k+1) = x2 - x1 ! Change of interface k+1, target - source

1652 endif

1653 enddo

1654

1655end subroutine dzfromh1h2

1656

1657!> Constructor for remapping control structure

1658subroutine initialize_remapping( CS, remapping_scheme, boundary_extrapolation, &

1659 check_reconstruction, check_remapping, force_bounds_in_subcell, &

1660 force_bounds_in_target, better_force_bounds_in_target, offset_tgt_summation, &

1661 om4_remap_via_sub_cells, answers_2018, answer_date, nk, &

1662 h_neglect, h_neglect_edge)

1663 ! Arguments

1664 type(remapping_cs), intent(inout) :: cs !< Remapping control structure

1665 character(len=*), intent(in) :: remapping_scheme !< Remapping scheme to use

1666 logical, optional, intent(in) :: boundary_extrapolation !< Indicate to extrapolate in boundary cells

1667 logical, optional, intent(in) :: check_reconstruction !< Indicate to check reconstructions

1668 logical, optional, intent(in) :: check_remapping !< Indicate to check results of remapping

1669 logical, optional, intent(in) :: force_bounds_in_subcell !< Force subcells values to be bounded

1670 logical, optional, intent(in) :: force_bounds_in_target !< Force target values to be bounded

1671 logical, optional, intent(in) :: better_force_bounds_in_target !< Force target values to be bounded

1672 logical, optional, intent(in) :: offset_tgt_summation !< Use an offset when summing sub-cells

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

1674 logical, optional, intent(in) :: answers_2018 !< If true use older, less accurate expressions.

1675 integer, optional, intent(in) :: answer_date !< The vintage of the expressions to use

1676 real, optional, intent(in) :: h_neglect !< A negligibly small width for the purpose of cell

1677 !! reconstructions in the same units as h0 [H]

1678 real, optional, intent(in) :: h_neglect_edge !< A negligibly small width for the purpose of edge

1679 !! value calculations in the same units as h0 [H].

1680 integer, optional, intent(in) :: nk !< Number of levels to initialize reconstruction class with

1681

1682 ! Note that remapping_scheme is mandatory for initialize_remapping()

1683 call remapping_set_param(cs, &

1684 remapping_scheme=remapping_scheme, &

1685 boundary_extrapolation=boundary_extrapolation, &

1686 check_reconstruction=check_reconstruction, &

1687 check_remapping=check_remapping, &

1688 force_bounds_in_subcell=force_bounds_in_subcell, &

1689 om4_remap_via_sub_cells=om4_remap_via_sub_cells, &

1690 force_bounds_in_target=force_bounds_in_target, &

1691 better_force_bounds_in_target=better_force_bounds_in_target, &

1692 offset_tgt_summation=offset_tgt_summation, &

1693 answers_2018=answers_2018, &

1694 answer_date=answer_date, &

1695 nk=nk, &

1696 h_neglect=h_neglect, &

1697 h_neglect_edge=h_neglect_edge)

1698

1699end subroutine initialize_remapping

1700

1701!> Changes the method of reconstruction

1702!! Use this routine to parse a string parameter specifying the reconstruction

1703!! and re-allocates work arrays appropriately. It is called from

1704!! initialize_remapping but can be called from an external module too.

1705subroutine setreconstructiontype(string,CS)

1706 character(len=*), intent(in) :: string !< String to parse for method

1707 type(remapping_cs), intent(inout) :: CS !< Remapping control structure

1708 ! Local variables

1709 integer :: degree

1710 degree = -99

1711 if (associated(cs%reconstruction)) then

1712 ! We have a choice of being careless and allowing easy re-use (e.g. when testing)

1713 cs%remapping_scheme = -911

1714 call cs%reconstruction%destroy()

1715 deallocate( cs%reconstruction )

1716 ! or being careful and make sure we've properly clean up...

1717 ! call MOM_error(FATAL, "setReconstructionType: "//&

1718 ! "Recon1d type is already associated when initializing.")

1719 endif

1720 select case ( uppercase(trim(string)) )

1721 case ("PCM")

1722 cs%remapping_scheme = remapping_pcm

1723 degree = 0

1724 case ("PLM")

1725 cs%remapping_scheme = remapping_plm

1726 degree = 1

1727 case ("PLM_HYBGEN")

1728 cs%remapping_scheme = remapping_plm_hybgen

1729 degree = 1

1730 case ("PPM_CW")

1731 cs%remapping_scheme = remapping_ppm_cw

1732 degree = 2

1733 case ("PPM_H4")

1734 cs%remapping_scheme = remapping_ppm_h4

1735 degree = 2

1736 case ("PPM_IH4")

1737 cs%remapping_scheme = remapping_ppm_ih4

1738 degree = 2

1739 case ("PPM_HYBGEN")

1740 cs%remapping_scheme = remapping_ppm_hybgen

1741 degree = 2

1742 case ("WENO_HYBGEN")

1743 cs%remapping_scheme = remapping_weno_hybgen

1744 degree = 2

1745 case ("PQM_IH4IH3")

1746 cs%remapping_scheme = remapping_pqm_ih4ih3

1747 degree = 4

1748 case ("PQM_IH6IH5")

1749 cs%remapping_scheme = remapping_pqm_ih6ih5

1750 degree = 4

1751 case ("C_PCM")

1752 allocate( pcm :: cs%reconstruction )

1753 cs%remapping_scheme = remapping_via_class

1754 case ("C_PLM_CW")

1755 allocate( plm_cw :: cs%reconstruction )

1756 cs%remapping_scheme = remapping_via_class

1757 case ("C_PLM_HYBGEN")

1758 allocate( plm_hybgen :: cs%reconstruction )

1759 cs%remapping_scheme = remapping_via_class

1760 case ("C_MPLM_WA")

1761 allocate( mplm_wa :: cs%reconstruction )

1762 cs%remapping_scheme = remapping_via_class

1763 case ("C_EMPLM_WA")

1764 allocate( emplm_wa :: cs%reconstruction )

1765 cs%remapping_scheme = remapping_via_class

1766 case ("C_MPLM_WA_POLY")

1767 allocate( mplm_wa_poly :: cs%reconstruction )

1768 cs%remapping_scheme = remapping_via_class

1769 case ("C_EMPLM_WA_POLY")

1770 allocate( emplm_wa_poly :: cs%reconstruction )

1771 cs%remapping_scheme = remapping_via_class

1772 case ("C_PLM_CWK")

1773 allocate( plm_cwk :: cs%reconstruction )

1774 cs%remapping_scheme = remapping_via_class

1775 case ("C_MPLM_CWK")

1776 allocate( mplm_cwk :: cs%reconstruction )

1777 cs%remapping_scheme = remapping_via_class

1778 case ("C_EMPLM_CWK")

1779 allocate( emplm_cwk :: cs%reconstruction )

1780 cs%remapping_scheme = remapping_via_class

1781 case ("C_PPM_CW")

1782 allocate( ppm_cw :: cs%reconstruction )

1783 cs%remapping_scheme = remapping_via_class

1784 case ("C_PPM_HYBGEN")

1785 allocate( ppm_hybgen :: cs%reconstruction )

1786 cs%remapping_scheme = remapping_via_class

1787 case ("C_PPM_CWK")

1788 allocate( ppm_cwk :: cs%reconstruction )

1789 cs%remapping_scheme = remapping_via_class

1790 case ("C_EPPM_CWK")

1791 allocate( eppm_cwk :: cs%reconstruction )

1792 cs%remapping_scheme = remapping_via_class

1793 case ("C_PPM_H4_2019")

1794 allocate( ppm_h4_2019 :: cs%reconstruction )

1795 cs%remapping_scheme = remapping_via_class

1796 case ("C_PPM_H4_2018")

1797 allocate( ppm_h4_2018 :: cs%reconstruction )

1798 cs%remapping_scheme = remapping_via_class

1799 case ("C_PLM_WLS")

1800 allocate( plm_wls :: cs%reconstruction )

1801 cs%remapping_scheme = remapping_via_class

1802 case default

1803 call mom_error(fatal, "setReconstructionType: "//&

1804 "Unrecognized choice for REMAPPING_SCHEME ("//trim(string)//").")

1805 end select

1806

1807 cs%degree = degree

1808

1809end subroutine setreconstructiontype

1810

1811!> Destrcutor for remapping control structure

1812subroutine end_remapping(CS)

1813 type(remapping_cs), intent(inout) :: cs !< Remapping control structure

1814

1815 cs%degree = 0

1816

1817end subroutine end_remapping

1818

1819!> Test if interpolate_column() produces the wrong answer

1820subroutine test_interp(test, msg, nsrc, h_src, u_src, ndest, h_dest, u_true)

1821 type(testing), intent(inout) :: test !< Unit testing convenience functions

1822 character(len=*), intent(in) :: msg !< Message to label test

1823 integer, intent(in) :: nsrc !< Number of source cells

1824 real, dimension(nsrc), intent(in) :: h_src !< Thickness of source cells [H]

1825 real, dimension(nsrc+1), intent(in) :: u_src !< Values at source cell interfaces [A]

1826 integer, intent(in) :: ndest !< Number of destination cells

1827 real, dimension(ndest), intent(in) :: h_dest !< Thickness of destination cells [H]

1828 real, dimension(ndest+1), intent(in) :: u_true !< Correct value at destination cell interfaces [A]

1829 ! Local variables

1830 real, dimension(ndest+1) :: u_dest ! Interpolated value at destination cell interfaces [A]

1831

1832 ! Interpolate from src to dest

1833 call interpolate_column(nsrc, h_src, u_src, ndest, h_dest, u_dest, .true.)

1834 call test%real_arr(ndest, u_dest, u_true, msg)

1835end subroutine test_interp

1836

1837!> Test if reintegrate_column() produces the wrong answer

1838subroutine test_reintegrate(test, msg, nsrc, h_src, uh_src, ndest, h_dest, uh_true)

1839 type(testing), intent(inout) :: test !< Unit testing convenience functions

1840 character(len=*), intent(in) :: msg !< Message to label test

1841 integer, intent(in) :: nsrc !< Number of source cells

1842 real, dimension(nsrc), intent(in) :: h_src !< Thickness of source cells [H]

1843 real, dimension(nsrc), intent(in) :: uh_src !< Values of source cell stuff [A H]

1844 integer, intent(in) :: ndest !< Number of destination cells

1845 real, dimension(ndest), intent(in) :: h_dest !< Thickness of destination cells [H]

1846 real, dimension(ndest), intent(in) :: uh_true !< Correct value of destination cell stuff [A H]

1847 ! Local variables

1848 real, dimension(ndest) :: uh_dest ! Reintegrated value on destination cells [A H]

1849

1850 ! Interpolate from src to dest

1851 call reintegrate_column(nsrc, h_src, uh_src, ndest, h_dest, uh_dest)

1852 call test%real_arr(ndest, uh_dest, uh_true, msg)

1853

1854end subroutine test_reintegrate

1855

1856!> Test class-based remapping for internal consistency on random data

1857subroutine test_recon_consistency(test, scheme, n0, niter, h_neglect)

1858 type(testing), intent(inout) :: test !< Unit testing convenience functions

1859 character(len=*), intent(in) :: scheme !< Name of scheme to use

1860 integer, intent(in) :: n0 !< Number of source cells

1861 integer, intent(in) :: niter !< Number of randomized columns to try

1862 real, intent(in) :: h_neglect !< A negligibly small width used in cell reconstructions [H]

1863 ! Local

1864 type(remapping_cs) :: remapCS !< Remapping control structure

1865 real :: h0(n0) ! Source grid [H but really nondim]

1866 real :: u0(n0) ! Source values [A]

1867 logical :: error ! Indicates a divergence

1868 integer :: iter ! Loop counter

1869 integer :: seed_size ! Number of integers used by seed

1870 integer, allocatable :: seed(:) ! Random number seed

1871 character(len=16) :: label ! Generated label

1872

1873 call initialize_remapping(remapcs, scheme, nk=n0, h_neglect=h_neglect, &

1874 force_bounds_in_subcell=.false. )

1875

1876 call random_seed(size=seed_size)

1877 allocate( seed(seed_size) )

1878 seed(:) = 102030405 ! Repeatable sequences

1879 call random_seed(put=seed)

1880

1881 error = .false.

1882 do iter = 1, niter

1883 call random_number( h0 ) ! In range 0-1

1884 h0(:) = max(0., h0(:) - 0.05) ! Make 5% of values equal to zero

1885 call random_number( u0 ) ! In range 0-1

1886

1887 call remapcs%reconstruction%reconstruct(h0, u0)

1888 if ( remapcs%reconstruction%check_reconstruction(h0, u0) ) then

1889 if ( .not. error ) then ! Only dump first error

1890 print *,'iter=',iter

1891 print *,'h0',h0

1892 print *,'u0',u0

1893 endif

1894 error = .true.

1895 endif

1896

1897 enddo

1898

1899 write(label,'(I0)') niter

1900 call test%test( error, trim(label)//' consistency tests of '//scheme )

1901

1902 call remapcs%reconstruction%destroy()

1903

1904end subroutine test_recon_consistency

1905

1906!> Test that remapping a uniform field remains uniform

1907subroutine test_preserve_uniform(test, scheme, n0, niter, h_neglect)

1908 type(testing), intent(inout) :: test !< Unit testing convenience functions

1909 character(len=*), intent(in) :: scheme !< Name of scheme to use

1910 integer, intent(in) :: n0 !< Number of source cells

1911 integer, intent(in) :: niter !< Number of randomized columns to try

1912 real, intent(in) :: h_neglect !< A negligibly small width used in cell reconstructions [H]

1913 ! Local

1914 type(remapping_cs) :: remapCS !< Remapping control structure

1915 real :: h0(n0), h1(n0) ! Source and target grids [H but really nondim]

1916 real :: u0(n0), u1(n0) ! Source and target values [A]

1917 logical :: error ! Indicates a divergence

1918 integer :: iter ! Loop counter

1919 integer :: seed_size ! Number of integers used by seed

1920 integer, allocatable :: seed(:) ! Random number seed

1921 character(len=16) :: label ! Generated label

1922

1923 call initialize_remapping(remapcs, scheme, nk=n0, h_neglect=h_neglect, &

1924 force_bounds_in_subcell=.true., &

1925 force_bounds_in_target=.true., &

1926 better_force_bounds_in_target=.true., &

1927 offset_tgt_summation=.false., &

1928 om4_remap_via_sub_cells=.false.)

1929

1930 call random_seed(size=seed_size)

1931 allocate( seed(seed_size) )

1932 seed(:) = 102030405 ! Repeatable sequences

1933 call random_seed(put=seed)

1934

1935 error = .false.

1936 do iter = 1, niter

1937 call random_number( h0 ) ! In range 0-1

1938 h0(:) = max(0., h0(:) - 0.05) ! Make 5% of values equal to zero

1939 call random_number( h1 ) ! In range 0-1

1940 h1(:) = max(0., h1(:) - 0.05) ! Make 5% of values equal to zero

1941 call random_number( u0(1) ) ! In range 0-1

1942 u0(:) = u0(1) ! Make u0 uniform

1943

1944 call remapping_core_h( remapcs, n0, h0, u0, n0, h1, u1 )

1945 if ( maxval( abs( u1(:) - u0(1) ) ) > 0. ) then

1946 if ( .not. error ) then ! Only dump first error

1947 print *,'iter=',iter

1948 print *,'u0(1)',u0(1)

1949 print *,'u1',u1

1950 print *,'u1-u0(1)',u1 - u0(1)

1951 endif

1952 error = .true.

1953 endif

1954

1955 enddo

1956

1957 write(label,'(I0)') niter

1958 call test%test( error, trim(label)//' uniformity tests of '//scheme )

1959

1960end subroutine test_preserve_uniform

1961

1962!> Test that remapping to the same grid preserves answers

1963!!

1964!! Notes:

1965!! 1) this test is currently imperfect since occasionally we see round-off

1966!! implying that ( A * B ) / A != B

1967!! 2) this test does not work for vanished layers

1968subroutine test_unchanged_grid(test, scheme, n0, niter, h_neglect)

1969 type(testing), intent(inout) :: test !< Unit testing convenience functions

1970 character(len=*), intent(in) :: scheme !< Name of scheme to use

1971 integer, intent(in) :: n0 !< Number of source cells

1972 integer, intent(in) :: niter !< Number of randomized columns to try

1973 real, intent(in) :: h_neglect !< A negligibly small width used in cell reconstructions [H]

1974 ! Local

1975 type(remapping_cs) :: remapCS !< Remapping control structure

1976 real :: h0(n0), h1(n0) ! Source and target grids [H but really nondim]

1977 real :: u0(n0), u1(n0) ! Source and target values [A]

1978 logical :: error ! Indicates a divergence

1979 integer :: iter ! Loop counter

1980 character(len=16) :: label ! Generated label

1981

1982 call initialize_remapping(remapcs, scheme, nk=n0, h_neglect=h_neglect, &

1983 force_bounds_in_subcell=.true., &

1984 force_bounds_in_target=.false., &

1985 better_force_bounds_in_target=.true., &

1986 offset_tgt_summation=.true., &

1987 om4_remap_via_sub_cells=.false.)

1988

1989 error = .false.

1990 do iter = 1, niter

1991 call random_number( h0 ) ! In range 0-1

1992 h0(:) = max(0., h0(:) - 0.00) ! Note we do NOT test with vanished layers

1993 h1(:) = h0(:) ! Exact copy

1994 call random_number( u0 ) ! In range 0-1

1995

1996 call remapping_core_h( remapcs, n0, h0, u0, n0, h1, u1 )

1997 if ( maxval( abs( u1(:) - u0(:) ) ) > epsilon(h0(1)) * maxval( abs( u0 ) ) ) then

1998 if ( .not. error ) then ! Only dump first error

1999 print *,'iter=',iter

2000 print *,'h0',h0

2001 print *,'u0',u0

2002 print *,'u1',u1

2003 print *,'u1-u0',u1 - u0

2004 endif

2005 error = .true.

2006 endif

2007

2008 enddo

2009

2010 write(label,'(I0)') niter

2011 call test%test( error, trim(label)//' unchanged grid tests of '//scheme )

2012

2013 call remapcs%reconstruction%destroy()

2014

2015end subroutine test_unchanged_grid

2016

2017!> Test class-based remapping bitwise reproduces original implementation

2018subroutine compare_two_schemes(test, CS1, CS2, n0, n1, niter, msg)

2019 type(testing), intent(inout) :: test !< Unit testing convenience functions

2020 type(remapping_cs), intent(inout) :: CS1 !< Remapping control structure configured for

2021 !! original implementation

2022 type(remapping_cs), intent(inout) :: CS2 !< Remapping control structure configured for

2023 !! class-based implementation

2024 integer, intent(in) :: n0 !< Number of source cells

2025 integer, intent(in) :: n1 !< Number of destination cells

2026 integer, intent(in) :: niter !< Number of randomized columns to try

2027 character(len=*), intent(in) :: msg !< Message to label test

2028 ! Local

2029 real :: h0(n0), h1(n1) ! Source and target grids [H but really nondim]

2030 real :: u0(n0), u1(n1), u2(n1) ! Source and two target values [A]

2031 logical :: error ! Indicates a divergence

2032 integer :: iter ! Loop counter

2033 integer :: seed_size ! Number of integers used by seed

2034 integer, allocatable :: seed(:) ! Random number seed

2035 character(len=16) :: label ! Generated label

2036

2037 call random_seed(size=seed_size)

2038 allocate( seed(seed_size) )

2039 seed(:) = 102030405 ! Repeatable sequences

2040 call random_seed(put=seed)

2041

2042 error = .false.

2043 do iter = 1, niter

2044 call random_number( h0 ) ! In range 0-1

2045 h0(:) = max(0., h0(:) - 0.00) ! Make 5% of values equal to zero

2046 h0(:) = h0(:) / sum( h0 ) ! Approximately normalize to total depth of 1

2047 call random_number(h1) ! In range 0-1

2048 h1(:) = max(0., h1(:) - 0.00) ! Make 5% of values equal to zero

2049 h1(:) = h1(:) / sum( h1 ) ! Approximately normalize to total depth of 1

2050 call random_number( u0 ) ! In range 0-1

2051

2052 call remapping_core_h( cs1, n0, h0, u0, n1, h1, u1 )

2053 call remapping_core_h( cs2, n0, h0, u0, n1, h1, u2 )

2054 error = sum( abs( u2(:) - u1(:) ) ) > 0.

2055 if (error) then

2056 print *,'iter=',iter

2057 print *,'h1',h1

2058 print *,'h0',h0

2059 print *,'u0',u0

2060 print *,'u1',u1

2061 print *,'u2',u2

2062 print *,'e',u2-u1

2063 ! CS1%debug = .true.

2064 ! call remapping_core_h( CS1, n0, h0, u0, n1, h1, u1 )

2065 ! CS2%debug = .true.

2066 ! call remapping_core_h( CS2, n0, h0, u0, n1, h1, u2 )

2067 exit

2068 endif

2069 enddo

2070

2071 write(label,'(I0)') niter

2072 call test%test( error, trim(label)//' comparisons of '//msg )

2073

2074end subroutine compare_two_schemes

2075

2076!> Runs unit tests on remapping functions.

2077!! Should only be called from a single/root thread

2078!! Returns True if a test fails, otherwise False

2079logical function remapping_unit_tests(verbose, num_comp_samp)

2080 logical, intent(in) :: verbose !< If true, write results to stdout

2081 integer, optional, intent(in) :: num_comp_samp !< If present, number of samples to

2082 !! try comparing class-based cade against OM4 code

2083 ! Local variables

2084 integer :: n0, n1, n2

2085 real, allocatable :: h0(:), h1(:), h2(:) ! Thicknesses for test columns [H]

2086 real, allocatable :: u0(:), u1(:), u2(:) ! Values for test profiles [A]

2087 real, allocatable :: dx1(:) ! Change in interface position [H]

2088 type(remapping_cs) :: cs, cs2 !< Remapping control structures

2089 real, allocatable, dimension(:,:) :: ppoly0_e ! Edge values of polynomials [A]

2090 real, allocatable, dimension(:,:) :: ppoly0_s ! Edge slopes of polynomials [A H-1]

2091 real, allocatable, dimension(:,:) :: ppoly0_coefs ! Coefficients of polynomials [A]

2092 real, allocatable, dimension(:) :: h_sub, h0_eff ! Subgrid and effective source thicknesses [H]

2093 real, allocatable, dimension(:) :: u_sub, uh_sub ! Subgrid values and totals [A, A H]

2094 real :: u02_err ! Error in remaping [A]

2095 integer, allocatable, dimension(:) :: isrc_start, isrc_end, isrc_max, itgt_start, itgt_end, isub_src ! Indices

2096 integer :: answer_date ! The vintage of the expressions to test

2097 real :: err ! Errors in the remapped thicknesses [H] or values [A]

2098 real :: h_neglect, h_neglect_edge ! Tiny thicknesses used in remapping [H]

2099 integer :: seed_size ! Number of integers used by seed

2100 integer, allocatable :: seed(:) ! Random number seed

2101 type(testing) :: test ! Unit testing convenience functions

2102 integer :: om4 ! Loop parameter, 0 or 1

2103 integer :: ntests ! Number of iterations when brute force testing

2104 character(len=4) :: om4_tag ! Generated label

2105 type(pcm) :: pcm

2106 type(plm_cw) :: plm_cw

2107 type(plm_hybgen) :: plm_hybgen

2108 type(mplm_wa) :: mplm_wa

2109 type(emplm_wa) :: emplm_wa

2110 type(mplm_wa_poly) :: mplm_wa_poly

2111 type(emplm_wa_poly) :: emplm_wa_poly

2112 type(plm_cwk) :: plm_cwk

2113 type(mplm_cwk) :: mplm_cwk

2114 type(emplm_cwk) :: emplm_cwk

2115 type(ppm_h4_2019) :: ppm_h4_2019

2116 type(ppm_h4_2018) :: ppm_h4_2018

2117 type(ppm_cw) :: ppm_cw

2118 type(ppm_hybgen) :: ppm_hybgen

2119 type(ppm_cwk) :: ppm_cwk

2120 type(eppm_cwk) :: eppm_cwk

2121 type(plm_wls) :: plm_wls

2122

2123 call test%set( verbose=verbose ) ! Sets the verbosity flag in test

2124! call test%set( stop_instantly=.true. ) ! While debugging

2125

2126 answer_date = 20190101 ! 20181231

2127 h_neglect = 1.0e-30

2128 h_neglect_edge = h_neglect ; if (answer_date < 20190101) h_neglect_edge = 1.0e-10

2129

2130 if (verbose) write(test%stdout,*) ' ===== MOM_remapping: remapping_unit_tests ================='

2131

2132 if (verbose) write(test%stdout,*) ' - - - - - 1st generation tests - - - - -'

2133

2134 call initialize_remapping(cs, 'PPM_H4', answer_date=answer_date, &

2135 h_neglect=h_neglect, h_neglect_edge=h_neglect_edge)

2136

2137 ! Profile 0: 4 layers of thickness 0.75 and total depth 3, with du/dz=8

2138 n0 = 4

2139 allocate( h0(n0), u0(n0) )

2140 h0 = (/0.75, 0.75, 0.75, 0.75/)

2141 u0 = (/9., 3., -3., -9./)

2142

2143 ! Profile 1: 3 layers of thickness 1.0 and total depth 3

2144 n1 = 3

2145 allocate( h1(n1), u1(n1), dx1(n1+1) )

2146 h1 = (/1.0, 1.0, 1.0/)

2147

2148 ! Profile 2: 6 layers of thickness 0.5 and total depth 3

2149 n2 = 6

2150 allocate( h2(n2), u2(n2) )

2151 h2 = (/0.5, 0.5, 0.5, 0.5, 0.5, 0.5/)

2152

2153 ! Mapping u1 from h1 to h2

2154 call dzfromh1h2( n0, h0, n1, h1, dx1 )

2155 call remapping_core_w( cs, n0, h0, u0, n1, dx1, u1 )

2156 call test%real_arr(3, u1, (/8.,0.,-8./), 'remapping_core_w() PPM_H4')

2157

2158 allocate(ppoly0_e(n0,2), ppoly0_s(n0,2), ppoly0_coefs(n0,cs%degree+1))

2159 ppoly0_e(:,:) = 0.0

2160 ppoly0_s(:,:) = 0.0

2161 ppoly0_coefs(:,:) = 0.0

2162

2163 call initialize_remapping(cs, 'PPM_H4', force_bounds_in_subcell=.false., answer_date=answer_date)

2164

2165 call remapping_core_h( cs, n0, h0, u0, n2, h2, u2, net_err=err )

2166 call test%real_arr(6, u2, (/10.,6.,2.,-2.,-6.,-10./), 'remapping_core_h() 2')

2167

2168 call remapping_core_h( cs, n0, h0, u0, 6, (/.125,.125,.125,.125,.125,.125/), u2, net_err=err )

2169 call test%real_arr(6, u2, (/11.5,10.5,9.5,8.5,7.5,6.5/), 'remapping_core_h() 3')

2170

2171 call remapping_core_h( cs, n0, h0, u0, 3, (/2.25,1.5,1./), u2, net_err=err )

2172 call test%real_arr(3, u2, (/3.,-10.5,-12./), 'remapping_core_h() 4')

2173

2174 deallocate(h0, u0, h1, u1, h2, u2, ppoly0_e, ppoly0_s, ppoly0_coefs)

2175 call end_remapping(cs)

2176

2177 ! ===============================================

2178 ! This section tests the reconstruction functions

2179 ! ===============================================

2180 if (verbose) write(test%stdout,*) ' - - - - - reconstruction tests - - - - -'

2181

2182 allocate( ppoly0_coefs(5,6), ppoly0_e(5,2), ppoly0_s(5,2), u2(2) )

2183

2184 call pcm_reconstruction(3, (/1.,2.,4./), &

2185 ppoly0_e(1:3,:), ppoly0_coefs(1:3,:) )

2186 call test%real_arr(3, ppoly0_e(:,1), (/1.,2.,4./), 'PCM: left edges')

2187 call test%real_arr(3, ppoly0_e(:,2), (/1.,2.,4./), 'PCM: right edges')

2188 call test%real_arr(3, ppoly0_coefs(:,1), (/1.,2.,4./), 'PCM: P0')

2189

2190 call plm_reconstruction(3, (/1.,1.,1./), (/1.,3.,5./), &

2191 ppoly0_e(1:3,:), ppoly0_coefs(1:3,:), h_neglect )

2192 call test%real_arr(3, ppoly0_e(:,1), (/1.,2.,5./), 'Unlim PLM: left edges')

2193 call test%real_arr(3, ppoly0_e(:,2), (/1.,4.,5./), 'Unlim PLM: right edges')

2194 call test%real_arr(3, ppoly0_coefs(:,1), (/1.,2.,5./), 'Unlim PLM: P0')

2195 call test%real_arr(3, ppoly0_coefs(:,2), (/0.,2.,0./), 'Unlim PLM: P1')

2196

2197 call plm_reconstruction(3, (/1.,1.,1./), (/1.,2.,7./), &

2198 ppoly0_e(1:3,:), ppoly0_coefs(1:3,:), h_neglect )

2199 call test%real_arr(3, ppoly0_e(:,1), (/1.,1.,7./), 'Left lim PLM: left edges')

2200 call test%real_arr(3, ppoly0_e(:,2), (/1.,3.,7./), 'Left lim PLM: right edges')

2201 call test%real_arr(3, ppoly0_coefs(:,1), (/1.,1.,7./), 'Left lim PLM: P0')

2202 call test%real_arr(3, ppoly0_coefs(:,2), (/0.,2.,0./), 'Left lim PLM: P1')

2203

2204 call plm_reconstruction(3, (/1.,1.,1./), (/1.,6.,7./), &

2205 ppoly0_e(1:3,:), ppoly0_coefs(1:3,:), h_neglect )

2206 call test%real_arr(3, ppoly0_e(:,1), (/1.,5.,7./), 'Right lim PLM: left edges')

2207 call test%real_arr(3, ppoly0_e(:,2), (/1.,7.,7./), 'Right lim PLM: right edges')

2208 call test%real_arr(3, ppoly0_coefs(:,1), (/1.,5.,7./), 'Right lim PLM: P0')

2209 call test%real_arr(3, ppoly0_coefs(:,2), (/0.,2.,0./), 'Right lim PLM: P1')

2210

2211 call plm_reconstruction(3, (/1.,2.,3./), (/1.,4.,9./), &

2212 ppoly0_e(1:3,:), ppoly0_coefs(1:3,:), h_neglect )

2213 call test%real_arr(3, ppoly0_e(:,1), (/1.,2.,9./), 'Non-uniform line PLM: left edges')

2214 call test%real_arr(3, ppoly0_e(:,2), (/1.,6.,9./), 'Non-uniform line PLM: right edges')

2215 call test%real_arr(3, ppoly0_coefs(:,1), (/1.,2.,9./), 'Non-uniform line PLM: P0')

2216 call test%real_arr(3, ppoly0_coefs(:,2), (/0.,4.,0./), 'Non-uniform line PLM: P1')

2217

2218 call edge_values_explicit_h4(5, (/1.,1.,1.,1.,1./), (/1.,3.,5.,7.,9./), &

2219 ppoly0_e, h_neglect=1e-10, answer_date=answer_date )

2220 ! The next two tests currently fail due to roundoff, but pass when given a reasonable tolerance.

2221 call test%real_arr(5, ppoly0_e(:,1), (/0.,2.,4.,6.,8./), 'Line H4: left edges', tol=8.0e-15)

2222 call test%real_arr(5, ppoly0_e(:,2), (/2.,4.,6.,8.,10./), 'Line H4: right edges', tol=1.0e-14)

2223

2224 ppoly0_e(:,1) = (/0.,2.,4.,6.,8./)

2225 ppoly0_e(:,2) = (/2.,4.,6.,8.,10./)

2226 call ppm_reconstruction(5, (/1.,1.,1.,1.,1./), (/1.,3.,5.,7.,9./), ppoly0_e(1:5,:), &

2227 ppoly0_coefs(1:5,:), h_neglect, answer_date=answer_date )

2228 call test%real_arr(5, ppoly0_coefs(:,1), (/1.,2.,4.,6.,9./), 'Line PPM: P0')

2229 call test%real_arr(5, ppoly0_coefs(:,2), (/0.,2.,2.,2.,0./), 'Line PPM: P1')

2230 call test%real_arr(5, ppoly0_coefs(:,3), (/0.,0.,0.,0.,0./), 'Line PPM: P2')

2231

2232 call edge_values_explicit_h4( 5, (/1.,1.,1.,1.,1./), (/1.,1.,7.,19.,37./), ppoly0_e, &

2233 h_neglect=1e-10, answer_date=answer_date )

2234 ! The next two tests are now passing when answer_date >= 20190101, but otherwise only work to roundoff.

2235 call test%real_arr(5, ppoly0_e(:,1), (/3.,0.,3.,12.,27./), 'Parabola H4: left edges', tol=2.7e-14)

2236 call test%real_arr(5, ppoly0_e(:,2), (/0.,3.,12.,27.,48./), 'Parabola H4: right edges', tol=4.8e-14)

2237 ppoly0_e(:,1) = (/0.,0.,3.,12.,27./)

2238 ppoly0_e(:,2) = (/0.,3.,12.,27.,48./)

2239 call ppm_reconstruction(5, (/1.,1.,1.,1.,1./), (/0.,1.,7.,19.,37./), ppoly0_e(1:5,:), &

2240 ppoly0_coefs(1:5,:), h_neglect, answer_date=answer_date )

2241 call test%real_arr(5, ppoly0_e(:,1), (/0.,0.,3.,12.,37./), 'Parabola PPM: left edges')

2242 call test%real_arr(5, ppoly0_e(:,2), (/0.,3.,12.,27.,37./), 'Parabola PPM: right edges')

2243 call test%real_arr(5, ppoly0_coefs(:,1), (/0.,0.,3.,12.,37./), 'Parabola PPM: P0')

2244 call test%real_arr(5, ppoly0_coefs(:,2), (/0.,0.,6.,12.,0./), 'Parabola PPM: P1')

2245 call test%real_arr(5, ppoly0_coefs(:,3), (/0.,3.,3.,3.,0./), 'Parabola PPM: P2')

2246

2247 ppoly0_e(:,1) = (/0.,0.,6.,10.,15./)

2248 ppoly0_e(:,2) = (/0.,6.,12.,17.,15./)

2249 call ppm_reconstruction(5, (/1.,1.,1.,1.,1./), (/0.,5.,7.,16.,15./), ppoly0_e(1:5,:), &

2250 ppoly0_coefs(1:5,:), h_neglect, answer_date=answer_date )

2251 call test%real_arr(5, ppoly0_e(:,1), (/0.,3.,6.,16.,15./), 'Limits PPM: left edges')

2252 call test%real_arr(5, ppoly0_e(:,2), (/0.,6.,9.,16.,15./), 'Limits PPM: right edges')

2253 call test%real_arr(5, ppoly0_coefs(:,1), (/0.,3.,6.,16.,15./), 'Limits PPM: P0')

2254 call test%real_arr(5, ppoly0_coefs(:,2), (/0.,6.,0.,0.,0./), 'Limits PPM: P1')

2255 call test%real_arr(5, ppoly0_coefs(:,3), (/0.,-3.,3.,0.,0./), 'Limits PPM: P2')

2256

2257 deallocate(ppoly0_e, ppoly0_s, ppoly0_coefs, u2)

2258

2259 ! ==============================================================

2260 ! This section tests the components of remapping_core_h()

2261 ! ==============================================================

2262

2263 if (verbose) write(test%stdout,*) ' - - - - - remapping algororithm tests - - - - -'

2264

2265 ! Test 1: n0=2, n1=3 Maps uniform grids with one extra target layer and no implicitly-vanished interior sub-layers

2266 ! h_src = | 3 | 3 |

2267 ! h_tgt = | 2 | 2 | 2 |

2268 ! h_sub = |0| 2 | 1 | 1 | 2 |0|

2269 ! isrc_start |1 | 4 |

2270 ! isrc_end | 3 | 5 |

2271 ! isrc_max | 2 | 5 |

2272 ! itgt_start |1 | 3 | 5 |

2273 ! itgt_end | 2 | 4 | 6|

2274 ! isub_src |1| 1 | 1 | 2 | 2 |2|

2275 allocate( h_sub(6), h0_eff(2), isrc_start(2), isrc_end(2), isrc_max(2), itgt_start(3), itgt_end(3), isub_src(6) )

2276 call intersect_src_tgt_grids( 2, (/3., 3./), & ! n0, h0

2277 3, (/2., 2., 2./), & ! n1, h1

2278 h_sub, h0_eff, &

2279 isrc_start, isrc_end, isrc_max, itgt_start, itgt_end, isub_src )

2280 if (verbose) write(test%stdout,*) "intersect_src_tgt_grids test 1: n0=2, n1=3"

2281 if (verbose) write(test%stdout,*) " h_src = | 3 | 3 |"

2282 if (verbose) write(test%stdout,*) " h_tgt = | 2 | 2 | 2 |"

2283 call test%real_arr(6, h_sub, (/0.,2.,1.,1.,2.,0./), 'h_sub')

2284 call test%real_arr(2, h0_eff, (/3.,3./), 'h0_eff')

2285 call test%int_arr(2, isrc_start, (/1,4/), 'isrc_start')

2286 call test%int_arr(2, isrc_end, (/3,5/), 'isrc_end')

2287 call test%int_arr(2, isrc_max, (/2,5/), 'isrc_max')

2288 call test%int_arr(3, itgt_start, (/1,3,5/), 'itgt_start')

2289 call test%int_arr(3, itgt_end, (/2,4,6/), 'itgt_end')

2290 call test%int_arr(6, isub_src, (/1,1,1,2,2,2/), 'isub_src')

2291 deallocate( h_sub, h0_eff, isrc_start, isrc_end, isrc_max, itgt_start, itgt_end, isub_src )

2292

2293 ! Test 2: n0=3, n1=2 Reverses "test 1" with more source than target layers

2294 ! h_src = | 2 | 2 | 2 |

2295 ! h_tgt = | 3 | 3 |

2296 ! h_sub = |0| 2 | 1 | 1 | 2 |0|

2297 ! isrc_start |1 | 3 | 5 |

2298 ! isrc_end | 2 | 4 | 5 |

2299 ! isrc_max | 2 | 4 | 5 |

2300 ! itgt_start |1 | 4 |

2301 ! itgt_end | 3 | 6|

2302 ! isub_src |1| 1 | 2 | 2 | 3 |3|

2303 allocate( h_sub(6), h0_eff(3), isrc_start(3), isrc_end(3), isrc_max(3), itgt_start(2), itgt_end(2), isub_src(6) )

2304 call intersect_src_tgt_grids( 3, (/2., 2., 2./), & ! n0, h0

2305 2, (/3., 3./), & ! n1, h1

2306 h_sub, h0_eff, &

2307 isrc_start, isrc_end, isrc_max, itgt_start, itgt_end, isub_src )

2308 if (verbose) write(test%stdout,*) "intersect_src_tgt_grids test 2: n0=3, n1=2"

2309 if (verbose) write(test%stdout,*) " h_src = | 2 | 2 | 2 |"

2310 if (verbose) write(test%stdout,*) " h_tgt = | 3 | 3 |"

2311 call test%real_arr(6, h_sub, (/0.,2.,1.,1.,2.,0./), 'h_sub')

2312 call test%real_arr(3, h0_eff, (/2.,2.,2./), 'h0_eff')

2313 call test%int_arr(3, isrc_start, (/1,3,5/), 'isrc_start')

2314 call test%int_arr(3, isrc_end, (/2,4,5/), 'isrc_end')

2315 call test%int_arr(3, isrc_max, (/2,4,5/), 'isrc_max')

2316 call test%int_arr(2, itgt_start, (/1,4/), 'itgt_start')

2317 call test%int_arr(2, itgt_end, (/3,6/), 'itgt_end')

2318 call test%int_arr(6, isub_src, (/1,1,2,2,3,3/), 'isub_src')

2319 deallocate( h_sub, h0_eff, isrc_start, isrc_end, isrc_max, itgt_start, itgt_end, isub_src )

2320

2321 ! Test 3: n0=2, n1=3 With aligned interfaces that lead to implicitly-vanished interior sub-layers

2322 n0 = 2 ; n1 = 3

2323 allocate( h0_eff(n0), isrc_start(n0), isrc_end(n0), isrc_max(n0), h0(n0), u0(n0) )

2324 allocate( itgt_start(n1), itgt_end(n1), h1(n1), u1(n1) )

2325 allocate( h_sub(n0+n1+1), isub_src(n0+n1+1) )

2326 u0 = (/ 2. , 5. /)

2327 h0 = (/ 2. , 4. /)

2328 h1 = (/ 2. , 2. , 2. /)

2329 ! h_src = |<- 2 ->|<- 4 ->|

2330 ! h_tgt = |<- 2 ->|<- 2 ->|<- 2 ->|

2331 ! h_sub = |0|<- 2 ->|0|<- 2 ->|<- 2 ->|0|

2332 ! isrc_start |1 |3 |

2333 ! isrc_end | 2 | 5 |

2334 ! isrc_max | 2 | 5 |

2335 ! itgt_start |1 | 4 | 5 |

2336 ! itgt_end | 3| 4 | 6|

2337 ! isub_src |1| 1 |2| 2 | 2 |2|

2338 call intersect_src_tgt_grids( n0, h0, n1, h1, h_sub, h0_eff, &

2339 isrc_start, isrc_end, isrc_max, itgt_start, itgt_end, isub_src )

2340 if (verbose) write(test%stdout,*) "intersect_src_tgt_grids test 3: n0=2, n1=3"

2341 if (verbose) write(test%stdout,*) " h_src = | 2 | 4 |"

2342 if (verbose) write(test%stdout,*) " h_tgt = | 2 | 2 | 2 |"

2343 call test%real_arr(6, h_sub, (/0.,2.,0.,2.,2.,0./), 'h_sub')

2344 call test%real_arr(2, h0_eff, (/2.,4./), 'h0_eff')

2345 call test%int_arr(2, isrc_start, (/1,3/), 'isrc_start')

2346 call test%int_arr(2, isrc_end, (/2,5/), 'isrc_end')

2347 call test%int_arr(2, isrc_max, (/2,5/), 'isrc_max')

2348 call test%int_arr(3, itgt_start, (/1,4,5/), 'itgt_start')

2349 call test%int_arr(3, itgt_end, (/3,4,6/), 'itgt_end')

2350 call test%int_arr(6, isub_src, (/1,1,2,2,2,2/), 'isub_src')

2351 allocate(ppoly0_coefs(n0,2), ppoly0_e(n0,2), ppoly0_s(n0,2))

2352 ! h_src = |<- 2 ->|<- 4 ->|

2353 ! h_sub = |0|<- 2 ->|0|<- 2 ->|<- 2 ->|0|

2354 ! u_src = | 2 | 5 |

2355 ! edge = |1 3|3 7|

2356 ! u_sub = |1| 2 |3| 4 | 6 |7|

2357 call plm_reconstruction(n0, h0, u0, ppoly0_e, ppoly0_coefs, h_neglect )

2358 call plm_boundary_extrapolation(n0, h0, u0, ppoly0_e, ppoly0_coefs, h_neglect)

2359 allocate(u_sub(n0+n1+1), uh_sub(n0+n1+1))

2360 call remap_src_to_sub_grid_om4(n0, h0, u0, ppoly0_e, ppoly0_coefs, &

2361 n1, h_sub, h0_eff, isrc_start, isrc_end, isrc_max, isub_src, &

2362 integration_plm, .false., u_sub, uh_sub, u02_err)

2363 call test%real_arr(6, u_sub, (/1.,2.,3.,4.,6.,7./), 'u_sub om4')

2364 call remap_src_to_sub_grid(n0, h0, u0, ppoly0_e, ppoly0_coefs, &

2365 n1, h_sub, isrc_start, isrc_end, isrc_max, isub_src, &

2366 integration_plm, .false., u_sub, uh_sub, u02_err)

2367 call test%real_arr(6, u_sub, (/1.,2.,3.,4.,6.,7./), 'u_sub')

2368 ! h_sub = |0|<- 2 ->|0|<- 2 ->|<- 2 ->|0|

2369 ! u_sub = |1| 2 |3| 4 | 6 |7|

2370 ! h_tgt = |<- 2 ->|<- 2 ->|<- 2 ->|

2371 ! u_tgt = | 2 | 4 | 6 |

2372 call remap_sub_to_tgt_grid(n0, n1, h1, h_sub, u_sub, uh_sub, itgt_start, itgt_end, &

2373 .false., .false., .false., u1, u02_err)

2374 call test%real_arr(3, u1, (/2.,4.,6./), 'u1')

2375 call remap_sub_to_tgt_grid(n0, n1, h1, h_sub, u_sub, uh_sub, itgt_start, itgt_end, &

2376 .true., .false., .false., u1, u02_err)

2377 call test%real_arr(3, u1, (/2.,4.,6./), 'u1.b')

2378 deallocate( ppoly0_coefs, ppoly0_e, ppoly0_s, u_sub, uh_sub, h0, u0, h1, u1)

2379 deallocate( h_sub, h0_eff, isrc_start, isrc_end, isrc_max, itgt_start, itgt_end, isub_src )

2380

2381 ! Test 4: n0=2, n1=3 Incomplete target column, sum(h_tgt)<sum(h_src), useful for diagnostics

2382 n0 = 2 ; n1 = 3

2383 allocate( h0_eff(n0), isrc_start(n0), isrc_end(n0), isrc_max(n0), h0(n0), u0(n0) )

2384 allocate( itgt_start(n1), itgt_end(n1), h1(n1), u1(n1) )

2385 allocate( h_sub(n0+n1+1), isub_src(n0+n1+1) )

2386 u0 = (/ 2. , 5. /)

2387 h0 = (/ 2. , 4. /)

2388 h1 = (/ 2. , 2. , 1. /)

2389 ! h_src = |<- 2 ->|<- 4 ->|

2390 ! h_tgt = |<- 2 ->|<- 2 ->|< 1 >|

2391 ! h_sub = |0|<- 2 ->|0|<- 2 ->|< 1 >|< 1 >|

2392 ! isrc_start |1 |3 |

2393 ! isrc_end | 2 | 6 |

2394 ! isrc_max | 2 | 4 |

2395 ! itgt_start |1 | 4 | 5 |

2396 ! itgt_end | 3| 4 | 5 |

2397 ! isub_src |1| 1 |2| 2 | 2 | 2 |

2398 call intersect_src_tgt_grids( n0, h0, n1, h1, h_sub, h0_eff, &

2399 isrc_start, isrc_end, isrc_max, itgt_start, itgt_end, isub_src )

2400 if (verbose) write(test%stdout,*) "intersect_src_tgt_grids test 4: n0=2, n1=3"

2401 if (verbose) write(test%stdout,*) " h_src = | 2 | 4 |"

2402 if (verbose) write(test%stdout,*) " h_tgt = | 2 | 2 | 1 |"

2403 call test%real_arr(6, h_sub, (/0.,2.,0.,2.,1.,1./), 'h_sub')

2404 call test%real_arr(2, h0_eff, (/2.,3./), 'h0_eff')

2405 call test%int_arr(2, isrc_start, (/1,3/), 'isrc_start')

2406 call test%int_arr(2, isrc_end, (/2,6/), 'isrc_end')

2407 call test%int_arr(2, isrc_max, (/2,4/), 'isrc_max')

2408 call test%int_arr(3, itgt_start, (/1,4,5/), 'itgt_start')

2409 call test%int_arr(3, itgt_end, (/3,4,5/), 'itgt_end')

2410 call test%int_arr(6, isub_src, (/1,1,2,2,2,2/), 'isub_src')

2411 allocate(ppoly0_coefs(n0,2), ppoly0_e(n0,2), ppoly0_s(n0,2))

2412 ! h_src = |<- 2 ->|<- 4 ->|

2413 ! h_sub = |0|<- 2 ->|0|<- 2 ->|<-1->|<-1->|

2414 ! u_src = | 2 | 5 |

2415 ! edge = |1 3|3 7|

2416 ! u_sub = |1| 2 |3| 4 | 5.5 | 6.5 |

2417 call plm_reconstruction(n0, h0, u0, ppoly0_e, ppoly0_coefs, h_neglect )

2418 call plm_boundary_extrapolation(n0, h0, u0, ppoly0_e, ppoly0_coefs, h_neglect)

2419 allocate(u_sub(n0+n1+1), uh_sub(n0+n1+1))

2420 call remap_src_to_sub_grid(2, (/2.,4./), (/2.,5./), ppoly0_e, ppoly0_coefs, &

2421 3, h_sub, isrc_start, isrc_end, isrc_max, isub_src, &

2422 integration_plm, .false., u_sub, uh_sub, u02_err)

2423 call test%real_arr(6, u_sub, (/1.,2.,3.,4.,5.5,6.5/), 'u_sub')

2424 deallocate( ppoly0_coefs, ppoly0_e, ppoly0_s, u_sub, uh_sub, h0, u0, h1, u1)

2425 deallocate( h_sub, h0_eff, isrc_start, isrc_end, isrc_max, itgt_start, itgt_end, isub_src )

2426

2427 ! Test 5: n0=3, n1=2 Target column exceeds source column, sum(h_tgt)>sum(h_src), useful for diagnostics

2428 n0 = 3 ; n1 = 2

2429 allocate( h0_eff(n0), isrc_start(n0), isrc_end(n0), isrc_max(n0), h0(n0), u0(n0) )

2430 allocate( itgt_start(n1), itgt_end(n1), h1(n1), u1(n1) )

2431 allocate( h_sub(n0+n1+1), isub_src(n0+n1+1) )

2432 u0 = (/ 2. , 4. , 5.5 /)

2433 h0 = (/ 2. , 2. , 1. /)

2434 h1 = (/ 2. , 4. /)

2435 ! h_src = |<- 2 ->|<- 2 ->|< 1 >|

2436 ! h_tgt = |<- 2 ->|<- 4 ->|

2437 ! h_sub = |0|<- 2 ->|0|<- 2 ->|< 1 >|< 1 >|

2438 ! isrc_start |1 |3 | 5 |

2439 ! isrc_end | 2 | 4 | 5 |

2440 ! isrc_max | 2 | 4 | 5 |

2441 ! itgt_start |1 | 4 |

2442 ! itgt_end | 3| 6 |

2443 ! isub_src |1| 1 |2| 2 | 3 | 3 |

2444 call intersect_src_tgt_grids( n0, h0, n1, h1, h_sub, h0_eff, &

2445 isrc_start, isrc_end, isrc_max, itgt_start, itgt_end, isub_src )

2446 if (verbose) write(test%stdout,*) "intersect_src_tgt_grids test 5: n0=3, n1=2"

2447 if (verbose) write(test%stdout,*) " h_src = | 2 | 2 | 1 |"

2448 if (verbose) write(test%stdout,*) " h_tgt = | 2 | 4 |"

2449 call test%real_arr(6, h_sub, (/0.,2.,0.,2.,1.,1./), 'h_sub')

2450 call test%real_arr(3, h0_eff, (/2.,2.,1./), 'h0_eff')

2451 call test%int_arr(3, isrc_start, (/1,3,5/), 'isrc_start')

2452 call test%int_arr(3, isrc_end, (/2,4,5/), 'isrc_end')

2453 call test%int_arr(3, isrc_max, (/2,4,5/), 'isrc_max')

2454 call test%int_arr(2, itgt_start, (/1,4/), 'itgt_start')

2455 call test%int_arr(2, itgt_end, (/3,6/), 'itgt_end')

2456 call test%int_arr(6, isub_src, (/1,1,2,2,3,3/), 'isub_src')

2457 allocate(ppoly0_coefs(n0,2), ppoly0_e(n0,2), ppoly0_s(n0,2))

2458 ! h_src = |<- 2 ->|<- 2 ->|< 1 >|

2459 ! h_sub = |0|<- 2 ->|0|<- 2 ->|< 1 >|< 1 >|

2460 ! u_src = | 2 | 4 | 5.5 |

2461 ! edge = |1 3|3 5|5 6|

2462 ! u_sub = |1| 2 |3| 4 | 5.5 | 6 |

2463 call plm_reconstruction(n0, h0, u0, ppoly0_e, ppoly0_coefs, h_neglect )

2464 call plm_boundary_extrapolation(n0, h0, u0, ppoly0_e, ppoly0_coefs, h_neglect)

2465 allocate(u_sub(n0+n1+1), uh_sub(n0+n1+1))

2466 call remap_src_to_sub_grid_om4(n0, h0, u0, ppoly0_e, ppoly0_coefs, &

2467 n1, h_sub, h0_eff, isrc_start, isrc_end, isrc_max, isub_src, &

2468 integration_plm, .false., u_sub, uh_sub, u02_err)

2469 call test%real_arr(6, u_sub, (/1.,2.,3.,4.,5.5,6./), 'u_sub om4')

2470 call remap_src_to_sub_grid(n0, h0, u0, ppoly0_e, ppoly0_coefs, &

2471 n1, h_sub, isrc_start, isrc_end, isrc_max, isub_src, &

2472 integration_plm, .false., u_sub, uh_sub, u02_err)

2473 call test%real_arr(6, u_sub, (/1.,2.,3.,4.,5.5,6./), 'u_sub')

2474 ! h_sub = |0|<- 2 ->|0|<- 2 ->|< 1 >|< 1 >|

2475 ! u_sub = |1| 2 |3| 4 | 5.5 | 6 |

2476 ! h_tgt = |<- 2 ->|<- 4 ->|

2477 ! u_tgt = | 2 | 4 7/8 |

2478 call remap_sub_to_tgt_grid(n0, n1, h1, h_sub, u_sub, uh_sub, itgt_start, itgt_end, &

2479 .false., .false., .false., u1, u02_err)

2480 call test%real_arr(2, u1, (/2.,4.875/), 'u1')

2481 call remap_sub_to_tgt_grid(n0, n1, h1, h_sub, u_sub, uh_sub, itgt_start, itgt_end, &

2482 .true., .false., .false., u1, u02_err)

2483 call test%real_arr(2, u1, (/2.,4.875/), 'u1.b')

2484 deallocate( ppoly0_coefs, ppoly0_e, ppoly0_s, u_sub, uh_sub, h0, u0, h1, u1)

2485 deallocate( h_sub, h0_eff, isrc_start, isrc_end, isrc_max, itgt_start, itgt_end, isub_src )

2486

2487 ! Test 6: n0=3, n1=5 Source and targets with vanished layers

2488 n0 = 3 ; n1 = 5

2489 allocate( h0_eff(n0), isrc_start(n0), isrc_end(n0), isrc_max(n0), h0(n0), u0(n0) )

2490 allocate( itgt_start(n1), itgt_end(n1), h1(n1), u1(n1) )

2491 allocate( h_sub(n0+n1+1), isub_src(n0+n1+1) )

2492 u0 = (/ 2. ,3., 4. /)

2493 h0 = (/ 2. ,0., 2. /)

2494 h1 = (/ 1. ,0., 1. ,0., 2. /)

2495 ! h_src = |<- 2 ->|0|<- 2 ->|

2496 ! h_tgt = |<- 1 ->|0|<- 1 ->|0|<- 2 ->|

2497 ! h_sub = |0|< 1 ->|0|< 1 >|0|0|0|<- 2 ->|0|

2498 ! isrc_start |1 |5|6 |

2499 ! isrc_end | 4 |5| 8 |

2500 ! isrc_max | 4 |5| 8 |

2501 ! itgt_start |1 |3| 4 |7| 8 |

2502 ! itgt_end | 2 |3| 6|7| 9|

2503 ! isub_src |1| 1 |1| 1 |2|3|3| 3 |3|

2504 call intersect_src_tgt_grids( n0, h0, n1, h1, h_sub, h0_eff, &

2505 isrc_start, isrc_end, isrc_max, itgt_start, itgt_end, isub_src )

2506 if (verbose) write(test%stdout,*) "intersect_src_tgt_grids test 6: n0=3, n1=5"

2507 if (verbose) write(test%stdout,*) " h_src = | 2 |0| 2 |"

2508 if (verbose) write(test%stdout,*) " h_tgt = | 1 |0| 1 |0| 2 |"

2509 call test%real_arr(9, h_sub, (/0.,1.,0.,1.,0.,0.,0.,2.,0./), 'h_sub')

2510 call test%real_arr(3, h0_eff, (/2.,0.,2./), 'h0_eff')

2511 call test%int_arr(3, isrc_start, (/1,5,6/), 'isrc_start')

2512 call test%int_arr(3, isrc_end, (/4,5,8/), 'isrc_end')

2513 call test%int_arr(3, isrc_max, (/4,5,8/), 'isrc_max')

2514 call test%int_arr(5, itgt_start, (/1,3,4,7,8/), 'itgt_start')

2515 call test%int_arr(5, itgt_end, (/2,3,6,7,9/), 'itgt_end')

2516 call test%int_arr(9, isub_src, (/1,1,1,1,2,3,3,3,3/), 'isub_src')

2517 allocate(ppoly0_coefs(n0,2), ppoly0_e(n0,2), ppoly0_s(n0,2))

2518 ! h_src = |<- 2 ->|0|<- 2 ->|

2519 ! h_sub = |0|< 1 ->|0|< 1 >|0|0|0|<- 2 ->|0|

2520 ! u_src = | 2 |3| 4 |

2521 ! edge = |1 3|3|3 5|

2522 ! u_sub = |1| 1.5 |2| 2.5 |3|3|3| 4 |5|

2523 call plm_reconstruction(n0, h0, u0, ppoly0_e, ppoly0_coefs, h_neglect )

2524 call plm_boundary_extrapolation(n0, h0, u0, ppoly0_e, ppoly0_coefs, h_neglect)

2525 allocate(u_sub(n0+n1+1), uh_sub(n0+n1+1))

2526 call remap_src_to_sub_grid_om4(n0, h0, u0, ppoly0_e, ppoly0_coefs, &

2527 n1, h_sub, h0_eff, isrc_start, isrc_end, isrc_max, isub_src, &

2528 integration_plm, .false., u_sub, uh_sub, u02_err)

2529 call test%real_arr(9, u_sub, (/1.,1.5,2.,2.5,3.,3.,3.,4.,5./), 'u_sub om4')

2530 call remap_src_to_sub_grid(n0, h0, u0, ppoly0_e, ppoly0_coefs, &

2531 n1, h_sub, isrc_start, isrc_end, isrc_max, isub_src, &

2532 integration_plm, .false., u_sub, uh_sub, u02_err)

2533 call test%real_arr(9, u_sub, (/1.,1.5,2.,2.5,3.,3.,3.,4.,5./), 'u_sub')

2534 ! h_sub = |0|< 1 ->|0|< 1 >|0|0|0|<- 2 ->|0|

2535 ! u_sub = |1| 1.5 |2| 2.5 |3|3|3| 4 |5|

2536 ! h_tgt = |<- 1 ->|0|<- 1 ->|0|<- 2 ->|

2537 ! u_tgt = | 1.5 |2| 2.5 |3| 4 |

2538 call remap_sub_to_tgt_grid(n0, n1, h1, h_sub, u_sub, uh_sub, itgt_start, itgt_end, &

2539 .false., .false., .false., u1, u02_err)

2540 call test%real_arr(5, u1, (/1.5,2.,2.5,3.,4./), 'u1')

2541 call remap_sub_to_tgt_grid(n0, n1, h1, h_sub, u_sub, uh_sub, itgt_start, itgt_end, &

2542 .true., .false., .false., u1, u02_err)

2543 call test%real_arr(5, u1, (/1.5,2.,2.5,3.,4./), 'u1.b')

2544 deallocate( ppoly0_coefs, ppoly0_e, ppoly0_s, u_sub, uh_sub, h0, u0, h1, u1)

2545 deallocate( h_sub, h0_eff, isrc_start, isrc_end, isrc_max, itgt_start, itgt_end, isub_src )

2546

2547 ! ============================================================

2548 ! This section tests remapping_core_h()

2549 ! ============================================================

2550 if (verbose) write(test%stdout,*) '- - - - - - - - - - remapping_core_h() tests - - - - - - - - -'

2551

2552 allocate(u2(2))

2553

2554 call initialize_remapping(cs, 'PLM', force_bounds_in_subcell=.false., answer_date=answer_date)

2555

2556 ! Remapping to just the two interior layers yields the same values as u_src(2:3)

2557 call remapping_core_h(cs, 4, (/0.,1.,1.,0./), (/5.,4.,2.,1./), 2, (/1.,1./), u2)

2558 call test%real_arr(2, u2, (/4.,2./), 'PLM: remapped h=0110->h=11 om4')

2559 call remapping_core_h(cs, 4, (/0.,1.,1.,0./), (/5.,4.,2.,1./), 2, (/1.,1./), u2)

2560 call test%real_arr(2, u2, (/4.,2./), 'PLM: remapped h=0110->h=11')

2561

2562 ! Remapping to two layers that are deeper. For the bottom layer of thickness 4,

2563 ! the first 1/4 has average 2, the remaining 3/4 has the bottom edge value or 1

2564 ! yield ing and average or 1.25

2565 call remapping_core_h(cs, 4, (/0.,1.,1.,0./), (/5.,4.,2.,1./), 2, (/1.,4./), u2)

2566 call test%real_arr(2, u2, (/4.,1.25/), 'PLM: remapped h=0110->h=14 om4')

2567 call remapping_core_h(cs, 4, (/0.,1.,1.,0./), (/5.,4.,2.,1./), 2, (/1.,4./), u2)

2568 call test%real_arr(2, u2, (/4.,1.25/), 'PLM: remapped h=0110->h=14')

2569

2570 ! Remapping to two layers with lowest layer not reach the bottom.

2571 ! Here, the bottom layer samples top half of source yeilding 2.5.

2572 ! Note: OM4 used the value as if the target layer was the same thickness as source.

2573 call remapping_set_param(cs, om4_remap_via_sub_cells=.true.)

2574 call remapping_core_h(cs, 4, (/0.,4.,4.,0./), (/5.,4.,2.,1./), 2, (/4.,2./), u2)

2575 call test%real_arr(2, u2, (/4.,2./), 'PLM: remapped h=0440->h=42 om4 (with known bug)')

2576 call remapping_set_param(cs, om4_remap_via_sub_cells=.false.)

2577 call remapping_core_h(cs, 4, (/0.,4.,4.,0./), (/5.,4.,2.,1./), 2, (/4.,2./), u2)

2578 call test%real_arr(2, u2, (/4.,2.5/), 'PLM: remapped h=0440->h=42')

2579

2580 ! Remapping to two layers with no layers sampling the bottom source layer

2581 ! The first layer samples the top half of u1, yielding 4.5

2582 ! The second layer samples the next quarter of u1, yielding 3.75

2583 call remapping_set_param(cs, om4_remap_via_sub_cells=.true.)

2584 call remapping_core_h(cs, 4, (/0.,5.,5.,0./), (/5.,4.,2.,1./), 2, (/2.,2./), u2)

2585 call test%real_arr(2, u2, (/4.5,3.5/), 'PLM: remapped h=0880->h=21 om4 (with known bug)')

2586 call remapping_set_param(cs, om4_remap_via_sub_cells=.false.)

2587 call remapping_core_h(cs, 4, (/0.,4.,4.,0./), (/5.,4.,2.,1./), 2, (/2.,1./), u2)

2588 call test%real_arr(2, u2, (/4.5,3.75/), 'PLM: remapped h=0440->h=21')

2589

2590 deallocate(u2)

2591

2592 ! Profile 0: 8 layers, 1x top/2x bottom vanished, and the rest with thickness 1.0, total depth 5, u(z) = 1 + z

2593 n0 = 8

2594 allocate( h0(n0), u0(n0) )

2595 h0 = (/0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0/)

2596 u0 = (/1.0, 1.5, 2.5, 3.5, 4.5, 5.5, 6.0, 6.0/)

2597 allocate( u1(8) )

2598

2599 call initialize_remapping(cs, 'PLM', answer_date=99990101, h_neglect=1.e-17, h_neglect_edge=1.e-2)

2600

2601 do om4 = 0, 1

2602 if ( om4 == 0 ) then

2603 cs%om4_remap_via_sub_cells = .false.

2604 om4_tag(:) = ' '

2605 else

2606 cs%om4_remap_via_sub_cells = .true.

2607 om4_tag(:) = ' om4'

2608 endif

2609

2610 ! Unchanged grid

2611 call remapping_core_h( cs, n0, h0, u0, 8, [0.,1.,1.,1.,1.,1.,0.,0.], u1)

2612 call test%real_arr(8, u1, (/1.0,1.5,2.5,3.5,4.5,5.5,6.0,6.0/), 'PLM: remapped h=01111100->h=01111100'//om4_tag)

2613

2614 ! Removing vanished layers (unchanged values for non-vanished layers, layer centers 0.5, 1.5, 2.5, 3.5, 4.5)

2615 call remapping_core_h( cs, n0, h0, u0, 5, [1.,1.,1.,1.,1.], u1)

2616 call test%real_arr(5, u1, (/1.5,2.5,3.5,4.5,5.5/), 'PLM: remapped h=01111100->h=11111'//om4_tag)

2617

2618 ! Remapping to variable thickness layers (layer centers 0.25, 1.0, 2.25, 4.0)

2619 call remapping_core_h( cs, n0, h0, u0, 4, [0.5,1.,1.5,2.], u1)

2620 call test%real_arr(4, u1, (/1.25,2.,3.25,5./), 'PLM: remapped h=01111100->h=h1t2'//om4_tag)

2621

2622 ! Remapping to variable thickness + vanished layers (layer centers 0.25, 1.0, 1.5, 2.25, 4.0)

2623 call remapping_core_h( cs, n0, h0, u0, 6, [0.5,1.,0.,1.5,2.,0.], u1)

2624 call test%real_arr(6, u1, (/1.25,2.,2.5,3.25,5.,6./), 'PLM: remapped h=01111100->h=h10t20'//om4_tag)

2625

2626 ! Remapping to deeper water column (layer centers 0.75, 2.25, 3., 5., 8.)

2627 call remapping_core_h( cs, n0, h0, u0, 5, [1.5,1.5,0.,4.,2.], u1)

2628 call test%real_arr(5, u1, (/1.75,3.25,4.,5.5,6./), 'PLM: remapped h=01111100->h=tt02'//om4_tag)

2629

2630 ! Remapping to slightly shorter water column (layer centers 0.5, 1.5, 2.5,, 3.5, 4.25)

2631 call remapping_core_h( cs, n0, h0, u0, 5, [1.,1.,1.,1.,0.5], u1)

2632 if ( om4 == 0 ) then

2633 call test%real_arr(5, u1, (/1.5,2.5,3.5,4.5,5.25/), 'PLM: remapped h=01111100->h=1111h')

2634 else

2635 call test%real_arr(5, u1, (/1.5,2.5,3.5,4.5,5.5/), 'PLM: remapped h=01111100->h=1111h om4 (known bug)')

2636 endif

2637

2638 ! Remapping to much shorter water column (layer centers 0.25, 0.5, 1.)

2639 call remapping_core_h( cs, n0, h0, u0, 3, [0.5,0.,1.], u1)

2640 if ( om4 == 0 ) then

2641 call test%real_arr(3, u1, (/1.25,1.5,2./), 'PLM: remapped h=01111100->h=h01')

2642 else

2643 call test%real_arr(3, u1, (/1.25,1.5,1.875/), 'PLM: remapped h=01111100->h=h01 om4 (known bug)')

2644 endif

2645

2646 enddo ! om4

2647

2648 call end_remapping(cs)

2649 deallocate( h0, u0, u1 )

2650

2651 ! ============================================================

2652 ! This section tests interpolation and reintegration functions

2653 ! ============================================================

2654 if (verbose) write(test%stdout,*) '- - - - - - - - - - interpolation tests - - - - - - - - -'

2655

2656 call test_interp(test, 'Identity: 3 layer', &

2657 3, (/1.,2.,3./), (/1.,2.,3.,4./), &

2658 3, (/1.,2.,3./), (/1.,2.,3.,4./) )

2659

2660 call test_interp(test, 'A: 3 layer to 2', &

2661 3, (/1.,1.,1./), (/1.,2.,3.,4./), &

2662 2, (/1.5,1.5/), (/1.,2.5,4./) )

2663

2664 call test_interp(test, 'B: 2 layer to 3', &

2665 2, (/1.5,1.5/), (/1.,4.,7./), &

2666 3, (/1.,1.,1./), (/1.,3.,5.,7./) )

2667

2668 call test_interp(test, 'C: 3 layer (vanished middle) to 2', &

2669 3, (/1.,0.,2./), (/1.,2.,2.,3./), &

2670 2, (/1.,2./), (/1.,2.,3./) )

2671

2672 call test_interp(test, 'D: 3 layer (deep) to 3', &

2673 3, (/1.,2.,3./), (/1.,2.,4.,7./), &

2674 2, (/2.,2./), (/1.,3.,5./) )

2675

2676 call test_interp(test, 'E: 3 layer to 3 (deep)', &

2677 3, (/1.,2.,4./), (/1.,2.,4.,8./), &

2678 3, (/2.,3.,4./), (/1.,3.,6.,8./) )

2679

2680 call test_interp(test, 'F: 3 layer to 4 with vanished top/botton', &

2681 3, (/1.,2.,4./), (/1.,2.,4.,8./), &

2682 4, (/0.,2.,5.,0./), (/0.,1.,3.,8.,0./) )

2683

2684 call test_interp(test, 'Fs: 3 layer to 4 with vanished top/botton (shallow)', &

2685 3, (/1.,2.,4./), (/1.,2.,4.,8./), &

2686 4, (/0.,2.,4.,0./), (/0.,1.,3.,7.,0./) )

2687

2688 call test_interp(test, 'Fd: 3 layer to 4 with vanished top/botton (deep)', &

2689 3, (/1.,2.,4./), (/1.,2.,4.,8./), &

2690 4, (/0.,2.,6.,0./), (/0.,1.,3.,8.,0./) )

2691

2692 if (verbose) write(test%stdout,*) ' - - - - - reintegration tests - - - - -'

2693

2694 call test_reintegrate(test, 'Identity: 3 layer', &

2695 3, (/1.,2.,3./), (/-5.,2.,1./), &

2696 3, (/1.,2.,3./), (/-5.,2.,1./) )

2697

2698 call test_reintegrate(test, 'A: 3 layer to 2', &

2699 3, (/2.,2.,2./), (/-5.,2.,1./), &

2700 2, (/3.,3./), (/-4.,2./) )

2701

2702 call test_reintegrate(test, 'A: 3 layer to 2 (deep)', &

2703 3, (/2.,2.,2./), (/-5.,2.,1./), &

2704 2, (/3.,4./), (/-4.,2./) )

2705

2706 call test_reintegrate(test, 'A: 3 layer to 2 (shallow)', &

2707 3, (/2.,2.,2./), (/-5.,2.,1./), &

2708 2, (/3.,2./), (/-4.,1.5/) )

2709

2710 call test_reintegrate(test, 'B: 3 layer to 4 with vanished top/bottom', &

2711 3, (/2.,2.,2./), (/-5.,2.,1./), &

2712 4, (/0.,3.,3.,0./), (/0.,-4.,2.,0./) )

2713

2714 call test_reintegrate(test, 'C: 3 layer to 4 with vanished top//middle/bottom', &

2715 3, (/2.,2.,2./), (/-5.,2.,1./), &

2716 5, (/0.,3.,0.,3.,0./), (/0.,-4.,0.,2.,0./) )

2717

2718 call test_reintegrate(test, 'D: 3 layer to 3 (vanished)', &

2719 3, (/2.,2.,2./), (/-5.,2.,1./), &

2720 3, (/0.,0.,0./), (/0.,0.,0./) )

2721

2722 call test_reintegrate(test, 'D: 3 layer (vanished) to 3', &

2723 3, (/0.,0.,0./), (/-5.,2.,1./), &

2724 3, (/2.,2.,2./), (/0., 0., 0./) )

2725

2726 call test_reintegrate(test, 'D: 3 layer (vanished) to 3 (vanished)', &

2727 3, (/0.,0.,0./), (/-5.,2.,1./), &

2728 3, (/0.,0.,0./), (/0.,0.,0./) )

2729

2730 call test_reintegrate(test, 'D: 3 layer (vanished) to 3 (vanished)', &

2731 3, (/0.,0.,0./), (/0.,0.,0./), &

2732 3, (/0.,0.,0./), (/0.,0.,0./) )

2733

2734 if (verbose) write(test%stdout,*) '- - - - - - - - - - Recon1d PCM tests - - - - - - - - -'

2735 call test%test( pcm%unit_tests(verbose, test%stdout, test%stderr), 'PCM unit test')

2736 call test%test( mplm_wa%unit_tests(verbose, test%stdout, test%stderr), 'MPLM_WA unit test')

2737 call test%test( emplm_wa%unit_tests(verbose, test%stdout, test%stderr), 'EMPLM_WA unit test')

2738 call test%test( mplm_wa_poly%unit_tests(verbose, test%stdout, test%stderr), 'MPLM_WA_poly unit test')

2739 call test%test( emplm_wa_poly%unit_tests(verbose, test%stdout, test%stderr), 'EMPLM_WA_poly unit test')

2740 call test%test( plm_hybgen%unit_tests(verbose, test%stdout, test%stderr), 'PLM_hybgen unit test')

2741 call test%test( plm_cw%unit_tests(verbose, test%stdout, test%stderr), 'PLM_CW unit test')

2742 call test%test( plm_cwk%unit_tests(verbose, test%stdout, test%stderr), 'PLM_CWK unit test')

2743 call test%test( mplm_cwk%unit_tests(verbose, test%stdout, test%stderr), 'MPLM_CWK unit test')

2744 call test%test( emplm_cwk%unit_tests(verbose, test%stdout, test%stderr), 'EMPLM_CWK unit test')

2745 call test%test( ppm_h4_2019%unit_tests(verbose, test%stdout, test%stderr), 'PPM_H4_2019 unit test')

2746 call test%test( ppm_h4_2018%unit_tests(verbose, test%stdout, test%stderr), 'PPM_H4_2018 unit test')

2747 call test%test( ppm_hybgen%unit_tests(verbose, test%stdout, test%stderr), 'PPM_hybgen unit test')

2748 call test%test( ppm_cw%unit_tests(verbose, test%stdout, test%stderr), 'PPM_CW unit test')

2749 call test%test( ppm_cwk%unit_tests(verbose, test%stdout, test%stderr), 'PPM_CWK unit test')

2750 call test%test( eppm_cwk%unit_tests(verbose, test%stdout, test%stderr), 'EPPM_CWK unit test')

2751 call test%test( plm_wls%unit_tests(verbose, test%stdout, test%stderr), 'PLM_WLS unit test')

2752

2753 ! Randomized, brute force tests

2754 ntests = 3000

2755 if (present(num_comp_samp)) ntests = num_comp_samp

2756

2757 call random_seed(size=seed_size)

2758 allocate( seed(seed_size) )

2759 seed(:) = 102030405

2760 call random_seed(put=seed)

2761

2762 n0 = 9

2763

2764 ! Internal consistency

2765 call test_recon_consistency(test, 'C_PCM', n0, ntests, h_neglect)

2766 call test_recon_consistency(test, 'C_PLM_CW', n0, ntests, h_neglect)

2767 call test_recon_consistency(test, 'C_PLM_HYBGEN', n0, ntests, h_neglect)

2768 call test_recon_consistency(test, 'C_MPLM_WA', n0, ntests, h_neglect)

2769 call test_recon_consistency(test, 'C_EMPLM_WA', n0, ntests, h_neglect)

2770 call test_recon_consistency(test, 'C_MPLM_WA_poly', n0, ntests, h_neglect)

2771 call test_recon_consistency(test, 'C_EMPLM_WA_poly', n0, ntests, h_neglect)

2772 call test_recon_consistency(test, 'C_PLM_CWK', n0, ntests, h_neglect)

2773 call test_recon_consistency(test, 'C_MPLM_CWK', n0, ntests, h_neglect)

2774 call test_recon_consistency(test, 'C_EMPLM_CWK', n0, ntests, h_neglect)

2775 call test_recon_consistency(test, 'C_PPM_H4_2018', n0, ntests, h_neglect)

2776 call test_recon_consistency(test, 'C_PPM_H4_2019', n0, ntests, h_neglect)

2777 call test_recon_consistency(test, 'C_PPM_HYBGEN', n0, ntests, h_neglect)

2778 call test_recon_consistency(test, 'C_PPM_CW', n0, ntests, h_neglect)

2779 call test_recon_consistency(test, 'C_PPM_CWK', n0, ntests, h_neglect)

2780 call test_recon_consistency(test, 'C_EPPM_CWK', n0, ntests, h_neglect)

2781 call test_recon_consistency(test, 'C_PLM_WLS', n0, ntests, h_neglect)

2782

2783 call test_preserve_uniform(test, 'PCM', n0, ntests, h_neglect)

2784 call test_preserve_uniform(test, 'C_PCM', n0, ntests, h_neglect)

2785! call test_preserve_uniform(test, 'PLM', n0, ntests, h_neglect) ! Fails

2786! call test_preserve_uniform(test, 'PLM_HYBGEN', n0, ntests, h_neglect) ! Fails

2787! call test_preserve_uniform(test, 'PPM_H4', n0, ntests, h_neglect) ! Fails

2788! call test_preserve_uniform(test, 'PPM_IH4', n0, ntests, h_neglect) ! Fails

2789! call test_preserve_uniform(test, 'PPM_HYBGEN', n0, ntests, h_neglect) ! Fails

2790! call test_preserve_uniform(test, 'PPM_CW', n0, ntests, h_neglect) ! Fails

2791! call test_preserve_uniform(test, 'WENO_HYBGEN', n0, ntests, h_neglect) ! Fails

2792! call test_preserve_uniform(test, 'PQM_IH4IH3', n0, ntests, h_neglect) ! Fails

2793 call test_preserve_uniform(test, 'C_PLM_CW', n0, ntests, h_neglect)

2794 call test_preserve_uniform(test, 'C_PLM_HYBGEN', n0, ntests, h_neglect)

2795 call test_preserve_uniform(test, 'C_MPLM_WA', n0, ntests, h_neglect)

2796 call test_preserve_uniform(test, 'C_EMPLM_WA', n0, ntests, h_neglect)

2797 call test_preserve_uniform(test, 'C_MPLM_WA_poly', n0, ntests, h_neglect) ! Surprised this passes -AJA

2798! call test_preserve_uniform(test, 'C_EMPLM_WA_poly', n0, ntests, h_neglect) ! This is known to fail

2799 call test_preserve_uniform(test, 'C_PLM_CWK', n0, ntests, h_neglect)

2800 call test_preserve_uniform(test, 'C_MPLM_CWK', n0, ntests, h_neglect)

2801 call test_preserve_uniform(test, 'C_EMPLM_CWK', n0, ntests, h_neglect)

2802 call test_preserve_uniform(test, 'C_PPM_H4_2019', n0, ntests, h_neglect)

2803 call test_preserve_uniform(test, 'C_PPM_H4_2018', n0, ntests, h_neglect)

2804 call test_preserve_uniform(test, 'C_PPM_HYBGEN', n0, ntests, h_neglect)

2805 call test_preserve_uniform(test, 'C_PPM_CW', n0, ntests, h_neglect)

2806 call test_preserve_uniform(test, 'C_PPM_CWK', n0, ntests, h_neglect)

2807 call test_preserve_uniform(test, 'C_EPPM_CWK', n0, ntests, h_neglect)

2808 call test_preserve_uniform(test, 'C_PLM_WLS', n0, ntests, h_neglect)

2809

2810 call test_unchanged_grid(test, 'C_PCM', n0, ntests, h_neglect)

2811 call test_unchanged_grid(test, 'C_PLM_CW', n0, ntests, h_neglect)

2812 call test_unchanged_grid(test, 'C_PLM_HYBGEN', n0, ntests, h_neglect)

2813 call test_unchanged_grid(test, 'C_PLM_CWK', n0, ntests, h_neglect)

2814 call test_unchanged_grid(test, 'C_MPLM_CWK', n0, ntests, h_neglect)

2815 call test_unchanged_grid(test, 'C_EMPLM_CWK', n0, ntests, h_neglect)

2816 call test_unchanged_grid(test, 'C_PPM_HYBGEN', n0, ntests, h_neglect)

2817 call test_unchanged_grid(test, 'C_PPM_CW', n0, ntests, h_neglect)

2818 call test_unchanged_grid(test, 'C_PPM_CWK', n0, ntests, h_neglect)

2819 call test_unchanged_grid(test, 'C_EPPM_CWK', n0, ntests, h_neglect)

2820 call test_unchanged_grid(test, 'C_PLM_WLS', n0, ntests, h_neglect)

2821

2822 ! Check that remapping to the exact same grid leaves values unchanged

2823 allocate( h0(8), u0(8) )

2824 h0 = (/0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0/)

2825 u0 = (/1.0, 1.5, 2.5, 3.5, 4.5, 5.5, 6.0, 6.0/)

2826 allocate( u1(8) )

2827 call initialize_remapping(cs, 'C_PLM_CW', nk=8)

2828 call remapping_core_h( cs, 8, h0, u0, 8, [0.,1.,1.,1.,1.,1.,0.,0.], u1 )

2829 call test%real_arr(8, u1, u0, 'remapping_core to unchanged grid with class')

2830

2831 call end_remapping(cs)

2832 deallocate( h0, u0, u1 )

2833

2834 ! Brute force test that we have bitwise identical answers with the new classes

2835 n0 = 7

2836 n1 = 4

2837

2838 ! PPM_CW and PPM_HYBGEN are identical, but are different options in build_reconstructions_1d()

2839 call initialize_remapping(cs, 'PPM_CW', answer_date=99990101, boundary_extrapolation=.false., &

2840 om4_remap_via_sub_cells=.false., force_bounds_in_subcell=.false., &

2841 h_neglect=h_neglect, h_neglect_edge=h_neglect_edge)

2842 call initialize_remapping(cs2, 'PPM_HYBGEN', answer_date=99990101, boundary_extrapolation=.false., &

2843 om4_remap_via_sub_cells=.false., force_bounds_in_subcell=.false., &

2844 h_neglect=h_neglect, h_neglect_edge=h_neglect_edge)

2845 call compare_two_schemes(test, cs, cs2, n0, n1, ntests, 'PPM_CW <-> PPM_HYBGEN')

2846

2847 ! PPM_CW <-> PPM_HYBGEN, as above but with OM4 subcells

2848 call initialize_remapping(cs, 'PPM_CW', answer_date=99990101, boundary_extrapolation=.false., &

2849 om4_remap_via_sub_cells=.true., force_bounds_in_subcell=.false., &

2850 h_neglect=h_neglect, h_neglect_edge=h_neglect_edge)

2851 call initialize_remapping(cs2, 'PPM_HYBGEN', answer_date=99990101, boundary_extrapolation=.false., &

2852 om4_remap_via_sub_cells=.true., force_bounds_in_subcell=.false., &

2853 h_neglect=h_neglect, h_neglect_edge=h_neglect_edge)

2854 call compare_two_schemes(test, cs, cs2, n0, n1, ntests, 'PPM_CW <-> PPM_HYBGEN OM4')

2855

2856 ! PPM_CW <-> PPM_HYBGEN, as above but with extrapolation

2857 call initialize_remapping(cs, 'PPM_CW', answer_date=99990101, boundary_extrapolation=.true., &

2858 om4_remap_via_sub_cells=.false., force_bounds_in_subcell=.false., &

2859 h_neglect=h_neglect, h_neglect_edge=h_neglect_edge)

2860 call initialize_remapping(cs2, 'PPM_HYBGEN', answer_date=99990101, boundary_extrapolation=.true., &

2861 om4_remap_via_sub_cells=.false., force_bounds_in_subcell=.false., &

2862 h_neglect=h_neglect, h_neglect_edge=h_neglect_edge)

2863 call compare_two_schemes(test, cs, cs2, n0, n1, ntests, 'PPM_CW <-> PPM_HYBGEN Extrap')

2864

2865 ! PPM_CW <-> PPM_HYBGEN, as above but with OM4 subcells and subcell bounds

2866 call initialize_remapping(cs, 'PPM_CW', answer_date=99990101, boundary_extrapolation=.false., &

2867 om4_remap_via_sub_cells=.true., force_bounds_in_subcell=.true., &

2868 h_neglect=h_neglect, h_neglect_edge=h_neglect_edge)

2869 call initialize_remapping(cs2, 'PPM_HYBGEN', answer_date=99990101, boundary_extrapolation=.false., &

2870 om4_remap_via_sub_cells=.true., force_bounds_in_subcell=.true., &

2871 h_neglect=h_neglect, h_neglect_edge=h_neglect_edge)

2872 call compare_two_schemes(test, cs, cs2, n0, n1, ntests, 'PPM_CW <-> PPM_HYBGEN')

2873

2874 ! PCM <-> C_PCM

2875 call initialize_remapping(cs, 'PCM', answer_date=99990101, om4_remap_via_sub_cells=.false., &

2876 force_bounds_in_subcell=.false.)

2877 call initialize_remapping(cs2, 'C_PCM', nk=n0, h_neglect=h_neglect)

2878 call compare_two_schemes(test, cs, cs2, n0, n1, ntests, 'PCM <-> C_PCM')

2879

2880 ! PLM <-> C_MPLM_WA_POLY

2881 call initialize_remapping(cs, 'PLM', answer_date=99990101, boundary_extrapolation=.false., &

2882 om4_remap_via_sub_cells=.false., force_bounds_in_subcell=.false., &

2883 h_neglect=h_neglect, h_neglect_edge=h_neglect_edge)

2884 call initialize_remapping(cs2, 'C_MPLM_WA_POLY', nk=n0, h_neglect=h_neglect)

2885 call compare_two_schemes(test, cs, cs2, n0, n1, ntests, 'PLM <-> C_MPLM_WA_poly')

2886

2887 ! PLM (with subcell bounds) <-> C_MPLM_WA_POLY

2888 call initialize_remapping(cs, 'PLM', answer_date=99990101, boundary_extrapolation=.false., &

2889 om4_remap_via_sub_cells=.false., force_bounds_in_subcell=.true., &

2890 h_neglect=h_neglect, h_neglect_edge=h_neglect_edge)

2891 call initialize_remapping(cs2, 'C_MPLM_WA_POLY', nk=n0, h_neglect=h_neglect)

2892 call compare_two_schemes(test, cs, cs2, n0, n1, ntests, 'PLM bounded <-> C_MPLM_WA_poly')

2893

2894 ! PLM + extrapolation <-> C_EMPLM_WA_POLY

2895 call initialize_remapping(cs, 'PLM', answer_date=99990101, boundary_extrapolation=.true., &

2896 om4_remap_via_sub_cells=.false., force_bounds_in_subcell=.false., &

2897 h_neglect=h_neglect, h_neglect_edge=h_neglect_edge)

2898 call initialize_remapping(cs2, 'C_EMPLM_WA_POLY', nk=n0, h_neglect=h_neglect)

2899 call compare_two_schemes(test, cs, cs2, n0, n1, ntests, 'PLM <-> C_EMPLM_WA_poly')

2900

2901 ! PLM + extrapolation (with subcell bounds) <-> C_EMPLM_WA_POLY

2902 call initialize_remapping(cs, 'PLM', answer_date=99990101, boundary_extrapolation=.true., &

2903 om4_remap_via_sub_cells=.false., force_bounds_in_subcell=.true., &

2904 h_neglect=h_neglect, h_neglect_edge=h_neglect_edge)

2905 call initialize_remapping(cs2, 'C_EMPLM_WA_POLY', nk=n0, h_neglect=h_neglect)

2906 call compare_two_schemes(test, cs, cs2, n0, n1, ntests, 'PLM bounded <-> C_EMPLM_WA_poly')

2907

2908 ! PPM_H4 (2018 answers) <-> C_PPM_H4_2018

2909 call initialize_remapping(cs, 'PPM_H4', answer_date=20180101, boundary_extrapolation=.false., &

2910 om4_remap_via_sub_cells=.false., force_bounds_in_subcell=.false., &

2911 h_neglect=h_neglect, h_neglect_edge=h_neglect_edge)

2912 call initialize_remapping(cs2, 'C_PPM_H4_2018', nk=n0, h_neglect=h_neglect)

2913 call compare_two_schemes(test, cs, cs2, n0, n1, ntests, 'PPM_H4 2018 <-> C_PPM_H4_2018')

2914

2915 ! PPM_H4 (2018 answers with subcell bounds) <-> C_PPM_H4_2018

2916 call initialize_remapping(cs, 'PPM_H4', answer_date=20180101, boundary_extrapolation=.false., &

2917 om4_remap_via_sub_cells=.false., force_bounds_in_subcell=.true., &

2918 h_neglect=h_neglect, h_neglect_edge=h_neglect_edge)

2919 call initialize_remapping(cs2, 'C_PPM_H4_2018', nk=n0, h_neglect=h_neglect)

2920 call compare_two_schemes(test, cs, cs2, n0, n1, ntests, 'PPM_H4 2018 bounded <-> C_PPM_H4_2018')

2921

2922 ! PPM_H4 (latest answers) <-> C_PPM_H4_2019

2923 call initialize_remapping(cs, 'PPM_H4', answer_date=99990101, boundary_extrapolation=.false., &

2924 om4_remap_via_sub_cells=.false., force_bounds_in_subcell=.false., &

2925 h_neglect=h_neglect, h_neglect_edge=h_neglect_edge)

2926 call initialize_remapping(cs2, 'C_PPM_H4_2019', nk=n0, h_neglect=h_neglect)

2927 call compare_two_schemes(test, cs, cs2, n0, n1, ntests, 'PPM_H4 <-> C_PPM_H4_2019')

2928

2929 ! PPM_H4 (latest answers with subcell bounds) <-> C_PPM_H4_2019

2930 call initialize_remapping(cs, 'PPM_H4', answer_date=99990101, boundary_extrapolation=.false., &

2931 om4_remap_via_sub_cells=.false., force_bounds_in_subcell=.true., &

2932 h_neglect=h_neglect, h_neglect_edge=h_neglect_edge)

2933 call initialize_remapping(cs2, 'C_PPM_H4_2019', nk=n0, h_neglect=h_neglect)

2934 call compare_two_schemes(test, cs, cs2, n0, n1, ntests, 'PPM_H4 bounded <-> C_PPM_H4_2019')

2935

2936 ! PLM_HYBGEN (latest answers with subcell bounds) <-> C_PLM_hybgen

2937 call initialize_remapping(cs, 'PLM_HYBGEN', answer_date=99990101, boundary_extrapolation=.false., &

2938 om4_remap_via_sub_cells=.false., force_bounds_in_subcell=.true., &

2939 h_neglect=h_neglect, h_neglect_edge=h_neglect_edge)

2940 call initialize_remapping(cs2, 'C_PLM_hybgen', nk=n0, h_neglect=h_neglect)

2941 call compare_two_schemes(test, cs, cs2, n0, n1, ntests, 'PLM_HYBGEN bounded <-> C_PLM_hygen')

2942

2943 ! PPM_HYBGEN (latest answers with subcell bounds) <-> C_PPM_hybgen

2944 call initialize_remapping(cs, 'PPM_HYBGEN', answer_date=99990101, boundary_extrapolation=.false., &

2945 om4_remap_via_sub_cells=.false., force_bounds_in_subcell=.true., &

2946 h_neglect=h_neglect, h_neglect_edge=h_neglect_edge)

2947 call initialize_remapping(cs2, 'C_PPM_HYBGEN', nk=n0, h_neglect=h_neglect)

2948 call compare_two_schemes(test, cs, cs2, n0, n1, ntests, 'PPM_HYBGEN bounded <-> C_PPM_HYGEN')

2949

2950 call end_remapping(cs)

2951 call end_remapping(cs2)

2952

2953 remapping_unit_tests = test%summarize('remapping_unit_tests')

2954

2955end function remapping_unit_tests

2956

2957end module mom_remapping