portedportable, not yet portedexecuted, not portableexecutable, not hit by this run
| 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 | |
| 6 | module MOM_remapping | |
| 7 | ||
| 8 | ! Original module written by Laurent White, 2008.06.09 | |
| 9 | ||
| 10 | use MOM_error_handler, only : MOM_error, FATAL | |
| 11 | use MOM_string_functions, only : uppercase | |
| 12 | use numerical_testing_type, only : testing | |
| 13 | use regrid_edge_values, only : edge_values_explicit_h4, edge_values_implicit_h4 | |
| 14 | use regrid_edge_values, only : edge_values_explicit_h4cw | |
| 15 | use regrid_edge_values, only : edge_values_implicit_h4, edge_values_implicit_h6 | |
| 16 | use regrid_edge_values, only : edge_slopes_implicit_h3, edge_slopes_implicit_h5 | |
| 17 | use PCM_functions, only : PCM_reconstruction | |
| 18 | use PLM_functions, only : PLM_reconstruction, PLM_boundary_extrapolation | |
| 19 | use PPM_functions, only : PPM_reconstruction, PPM_boundary_extrapolation | |
| 20 | use PPM_functions, only : PPM_monotonicity | |
| 21 | use PQM_functions, only : PQM_reconstruction, PQM_boundary_extrapolation_v1 | |
| 22 | use MOM_hybgen_remap, only : hybgen_plm_coefs, hybgen_ppm_coefs, hybgen_weno_coefs | |
| 23 | ||
| 24 | use Recon1d_type, only : Recon1d | |
| 25 | use Recon1d_PCM, only : PCM | |
| 26 | use Recon1d_PLM_CW, only : PLM_CW | |
| 27 | use Recon1d_PLM_hybgen, only : PLM_hybgen | |
| 28 | use Recon1d_PLM_CWK, only : PLM_CWK | |
| 29 | use Recon1d_MPLM_CWK, only : MPLM_CWK | |
| 30 | use Recon1d_EMPLM_CWK, only : EMPLM_CWK | |
| 31 | use Recon1d_MPLM_WA, only : MPLM_WA | |
| 32 | use Recon1d_EMPLM_WA, only : EMPLM_WA | |
| 33 | use Recon1d_MPLM_WA_poly, only : MPLM_WA_poly | |
| 34 | use Recon1d_EMPLM_WA_poly, only : EMPLM_WA_poly | |
| 35 | use Recon1d_PPM_CW, only : PPM_CW | |
| 36 | use Recon1d_PPM_hybgen, only : PPM_hybgen | |
| 37 | use Recon1d_PPM_CWK, only : PPM_CWK | |
| 38 | use Recon1d_EPPM_CWK, only : EPPM_CWK | |
| 39 | use Recon1d_PPM_H4_2019, only : PPM_H4_2019 | |
| 40 | use Recon1d_PPM_H4_2018, only : PPM_H4_2018 | |
| 41 | use Recon1d_PLM_WLS, only : PLM_WLS | |
| 42 | ||
| 43 | implicit none ; private | |
| 44 | ||
| 45 | !> Container for remapping parameters | |
| 46 | type, 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() | |
| 84 | end type | |
| 85 | ||
| 86 | ! The following routines are visible to the outside world | |
| 87 | public remapping_core_h, remapping_core_w | |
| 88 | public initialize_remapping, end_remapping, remapping_set_param, extract_member_remapping_CS | |
| 89 | public remapping_unit_tests, build_reconstructions_1d, average_value_ppoly | |
| 90 | public interpolate_column, reintegrate_column, dzFromH1H2 | |
| 91 | ||
| 92 | ! The following are private parameter constants | |
| 93 | integer, parameter :: REMAPPING_PCM = 0 !< O(h^1) remapping scheme | |
| 94 | integer, parameter :: REMAPPING_PLM = 2 !< O(h^2) remapping scheme | |
| 95 | integer, parameter :: REMAPPING_PLM_HYBGEN = 3 !< O(h^2) remapping scheme | |
| 96 | integer, parameter :: REMAPPING_PPM_CW =10 !< O(h^3) remapping scheme | |
| 97 | integer, parameter :: REMAPPING_PPM_H4 = 4 !< O(h^3) remapping scheme | |
| 98 | integer, parameter :: REMAPPING_PPM_IH4 = 5 !< O(h^3) remapping scheme | |
| 99 | integer, parameter :: REMAPPING_PPM_HYBGEN = 6 !< O(h^3) remapping scheme | |
| 100 | integer, parameter :: REMAPPING_WENO_HYBGEN= 7 !< O(h^3) remapping scheme | |
| 101 | integer, parameter :: REMAPPING_PQM_IH4IH3 = 8 !< O(h^4) remapping scheme | |
| 102 | integer, parameter :: REMAPPING_PQM_IH6IH5 = 9 !< O(h^5) remapping scheme | |
| 103 | integer, parameter :: REMAPPING_VIA_CLASS =99 !< Scheme is controlled by Recon1d class | |
| 104 | ||
| 105 | integer, parameter :: INTEGRATION_PCM = 0 !< Piecewise Constant Method | |
| 106 | integer, parameter :: INTEGRATION_PLM = 1 !< Piecewise Linear Method | |
| 107 | integer, parameter :: INTEGRATION_PPM = 3 !< Piecewise Parabolic Method | |
| 108 | integer, parameter :: INTEGRATION_PQM = 5 !< Piecewise Quartic Method | |
| 109 | ||
| 110 | !> Documentation for external callers | |
| 111 | character(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" | |
| 121 | character(len=3), public :: remappingDefaultScheme = "PLM" !< Default remapping method | |
| 122 | ||
| 123 | contains | |
| 124 | ||
| 125 | !> Set parameters within remapping object | |
| 126 | 6 | subroutine 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 | 6 | if (present(remapping_scheme)) then |
| 152 | 5 | call setReconstructionType( remapping_scheme, CS ) |
| 153 | 5 | if (index(trim(remapping_scheme),'C_')>0) then |
| 154 | 0 | if (present(nk)) then |
| 155 | 0 | call CS%reconstruction%init(nk, h_neglect=h_neglect) |
| 156 | else | |
| 157 | call MOM_error( FATAL, 'MOM_remapping, remapping_set_param: '//& | |
| 158 | 0 | 'Using the Recon1d class for remapping requires nk to be passed' ) |
| 159 | endif | |
| 160 | endif | |
| 161 | endif | |
| 162 | 6 | if (present(boundary_extrapolation)) then |
| 163 | 6 | CS%boundary_extrapolation = boundary_extrapolation |
| 164 | endif | |
| 165 | 6 | if (present(check_reconstruction)) then |
| 166 | 2 | CS%check_reconstruction = check_reconstruction |
| 167 | endif | |
| 168 | 6 | if (present(check_remapping)) then |
| 169 | 2 | CS%check_remapping = check_remapping |
| 170 | endif | |
| 171 | 6 | if (present(force_bounds_in_subcell)) then |
| 172 | 2 | CS%force_bounds_in_subcell = force_bounds_in_subcell |
| 173 | endif | |
| 174 | 6 | if (present(force_bounds_in_target)) then |
| 175 | 0 | CS%force_bounds_in_target = force_bounds_in_target |
| 176 | endif | |
| 177 | 6 | if (present(better_force_bounds_in_target)) then |
| 178 | 0 | CS%better_force_bounds_in_target = better_force_bounds_in_target |
| 179 | endif | |
| 180 | 6 | if (present(offset_tgt_summation)) then |
| 181 | 0 | CS%offset_tgt_summation = offset_tgt_summation |
| 182 | endif | |
| 183 | 6 | if (present(om4_remap_via_sub_cells)) then |
| 184 | 5 | CS%om4_remap_via_sub_cells = om4_remap_via_sub_cells |
| 185 | endif | |
| 186 | 6 | if (present(answers_2018)) then |
| 187 | 0 | if (answers_2018) then |
| 188 | 0 | CS%answer_date = 20181231 |
| 189 | else | |
| 190 | 0 | CS%answer_date = 20190101 |
| 191 | endif | |
| 192 | endif | |
| 193 | 6 | if (present(answer_date)) then |
| 194 | 5 | CS%answer_date = answer_date |
| 195 | endif | |
| 196 | 6 | if (present(h_neglect)) then |
| 197 | 5 | CS%h_neglect = h_neglect |
| 198 | endif | |
| 199 | 6 | if (present(h_neglect_edge)) then |
| 200 | 5 | CS%h_neglect_edge = h_neglect_edge |
| 201 | endif | |
| 202 | ||
| 203 | 6 | end subroutine remapping_set_param |
| 204 | ||
| 205 | 0 | subroutine 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 | 0 | if (present(remapping_scheme)) remapping_scheme = CS%remapping_scheme |
| 222 | 0 | if (present(degree)) degree = CS%degree |
| 223 | 0 | if (present(boundary_extrapolation)) boundary_extrapolation = CS%boundary_extrapolation |
| 224 | 0 | if (present(check_reconstruction)) check_reconstruction = CS%check_reconstruction |
| 225 | 0 | if (present(check_remapping)) check_remapping = CS%check_remapping |
| 226 | 0 | if (present(force_bounds_in_subcell)) force_bounds_in_subcell = CS%force_bounds_in_subcell |
| 227 | 0 | if (present(force_bounds_in_target)) force_bounds_in_target = CS%force_bounds_in_target |
| 228 | 0 | if (present(better_force_bounds_in_target)) better_force_bounds_in_target = CS%better_force_bounds_in_target |
| 229 | 0 | if (present(offset_tgt_summation)) offset_tgt_summation = CS%offset_tgt_summation |
| 230 | ||
| 231 | 0 | end 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 | |
| 238 | 195546 | subroutine 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 | 1253992 | real, dimension(n0+n1+1) :: h_sub ! Width of each each sub-cell [H] |
| 251 | 1253992 | real, dimension(n0+n1+1) :: uh_sub ! Integral of u*h over each sub-cell [A H] |
| 252 | 1253992 | real, dimension(n0+n1+1) :: u_sub ! Average of u over each sub-cell [A] |
| 253 | 1253992 | integer, dimension(n0+n1+1) :: isub_src ! Index of source cell for each sub-cell |
| 254 | 1253992 | integer, dimension(n0) :: isrc_start ! Index of first sub-cell within each source cell |
| 255 | 1253992 | integer, dimension(n0) :: isrc_end ! Index of last sub-cell within each source cell |
| 256 | 1253992 | integer, dimension(n0) :: isrc_max ! Index of thickest sub-cell within each source cell |
| 257 | 1253992 | real, dimension(n0) :: h0_eff ! Effective thickness of source cells [H] |
| 258 | 1253992 | integer, dimension(n1) :: itgt_start ! Index of first sub-cell within each target cell |
| 259 | 1253992 | 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 | 1253992 | real, dimension(n0,2) :: ppoly_r_E ! Edge value of polynomial [A] |
| 263 | 1253992 | real, dimension(n0,2) :: ppoly_r_S ! Edge slope of polynomial [A H-1] |
| 264 | 1253992 | 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 | 626996 | isrc_start, isrc_end, isrc_max, itgt_start, itgt_end, isub_src) |
| 272 | ||
| 273 | 626996 | if (CS%remapping_scheme == REMAPPING_VIA_CLASS) then |
| 274 | ||
| 275 | ! if (CS%debug) call CS%reconstruction%set_debug() ! Sets an internal flag | |
| 276 | ||
| 277 | 0 | 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 | 0 | 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 | 0 | 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 | 0 | 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 | 626996 | CS%h_neglect, CS%h_neglect_edge, PCM_cell, debug=CS%debug) |
| 303 | ||
| 304 | 626996 | if (CS%check_reconstruction) call check_reconstructions_1d(n0, h0, u0, CS%degree, & |
| 305 | 0 | 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 | 626996 | 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 | 626996 | 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 | 0 | 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 | 626996 | 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 | 626996 | uh_err = uh_err + u02_err |
| 331 | ||
| 332 | 626996 | if (CS%check_remapping) call check_remapped_values(n0, h0, u0, ppoly_r_E, CS%degree, ppoly_r_coefs, & |
| 333 | 0 | n1, h1, u1, iMethod, uh_err, "remapping_core_h") |
| 334 | ||
| 335 | endif | |
| 336 | ||
| 337 | 626996 | if (present(net_err)) net_err = uh_err |
| 338 | ||
| 339 | 626996 | end 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. | |
| 343 | 0 | subroutine 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 | 0 | real, dimension(n0+n1+1) :: h_sub ! Width of each each sub-cell [H] |
| 354 | 0 | real, dimension(n0+n1+1) :: uh_sub ! Integral of u*h over each sub-cell [A H] |
| 355 | 0 | real, dimension(n0+n1+1) :: u_sub ! Average of u over each sub-cell [A] |
| 356 | 0 | integer, dimension(n0+n1+1) :: isub_src ! Index of source cell for each sub-cell |
| 357 | 0 | integer, dimension(n0) :: isrc_start ! Index of first sub-cell within each source cell |
| 358 | 0 | integer, dimension(n0) :: isrc_end ! Index of last sub-cell within each source cell |
| 359 | 0 | integer, dimension(n0) :: isrc_max ! Index of thickest sub-cell within each source cell |
| 360 | 0 | real, dimension(n0) :: h0_eff ! Effective thickness of source cells [H] |
| 361 | 0 | integer, dimension(n1) :: itgt_start ! Index of first sub-cell within each target cell |
| 362 | 0 | 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 | 0 | real, dimension(n0,2) :: ppoly_r_E ! Edge value of polynomial [A] |
| 366 | 0 | real, dimension(n0,2) :: ppoly_r_S ! Edge slope of polynomial [A H-1] |
| 367 | 0 | real, dimension(n0,CS%degree+1) :: ppoly_r_coefs ! Coefficients of polynomial reconstructions [A] |
| 368 | 0 | 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 | 0 | CS%h_neglect, CS%h_neglect_edge ) |
| 375 | ||
| 376 | 0 | if (CS%check_reconstruction) call check_reconstructions_1d(n0, h0, u0, CS%degree, & |
| 377 | 0 | CS%boundary_extrapolation, ppoly_r_coefs, ppoly_r_E) |
| 378 | ||
| 379 | ! This is a temporary step prior to switching to remapping_core_h() | |
| 380 | 0 | do k = 1, n1 |
| 381 | 0 | if (k<=n0) then |
| 382 | 0 | h1(k) = max( 0., h0(k) + ( dx(k+1) - dx(k) ) ) |
| 383 | else | |
| 384 | 0 | 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 | 0 | 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 | 0 | 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 | 0 | 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 | 0 | uh_err = uh_err + u02_err |
| 407 | ||
| 408 | 0 | if (CS%check_remapping) call check_remapped_values(n0, h0, u0, ppoly_r_E, CS%degree, ppoly_r_coefs, & |
| 409 | 0 | n1, h1, u1, iMethod, uh_err, "remapping_core_w") |
| 410 | ||
| 411 | 0 | end subroutine remapping_core_w |
| 412 | ||
| 413 | !> Creates polynomial reconstructions of u0 on the source grid h0. | |
| 414 | 626996 | subroutine build_reconstructions_1d( CS, n0, h0, u0, ppoly_r_coefs, & |
| 415 | 626996 | ppoly_r_E, ppoly_r_S, iMethod, h_neglect, & |
| 416 | 195546 | 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 | 626996 | deb = .false. ; if (present(debug)) deb = debug |
| 444 | ||
| 445 | 626996 | h_neg_edge = h_neglect ; if (present(h_neglect_edge)) h_neg_edge = h_neglect_edge |
| 446 | ||
| 447 | ! Reset polynomial | |
| 448 | 95930388 | ppoly_r_E(:,:) = 0.0 |
| 449 | 95930388 | ppoly_r_S(:,:) = 0.0 |
| 450 | 118996692 | ppoly_r_coefs(:,:) = 0.0 |
| 451 | 626996 | iMethod = -999 |
| 452 | ||
| 453 | 626996 | local_remapping_scheme = CS%remapping_scheme |
| 454 | 626996 | if (n0<=1) then |
| 455 | 0 | local_remapping_scheme = REMAPPING_PCM |
| 456 | 626996 | elseif (n0<=3) then |
| 457 | 0 | local_remapping_scheme = min( local_remapping_scheme, REMAPPING_PLM ) |
| 458 | 626996 | elseif (n0<=4 .and. local_remapping_scheme /= REMAPPING_PPM_CW ) then |
| 459 | 0 | local_remapping_scheme = min( local_remapping_scheme, REMAPPING_PPM_H4 ) |
| 460 | endif | |
| 461 | 0 | select case ( local_remapping_scheme ) |
| 462 | case ( REMAPPING_PCM ) | |
| 463 | 0 | call PCM_reconstruction( n0, u0, ppoly_r_E, ppoly_r_coefs) |
| 464 | 0 | iMethod = INTEGRATION_PCM |
| 465 | case ( REMAPPING_PLM ) | |
| 466 | 323492 | call PLM_reconstruction( n0, h0, u0, ppoly_r_E, ppoly_r_coefs, h_neglect ) |
| 467 | 323492 | if ( CS%boundary_extrapolation ) then |
| 468 | 0 | call PLM_boundary_extrapolation( n0, h0, u0, ppoly_r_E, ppoly_r_coefs, h_neglect) |
| 469 | endif | |
| 470 | 323492 | iMethod = INTEGRATION_PLM |
| 471 | case ( REMAPPING_PLM_HYBGEN ) | |
| 472 | 0 | call hybgen_PLM_coefs(u0, h0, ppoly_r_coefs(:,2), n0, 1, h_neglect) |
| 473 | 0 | do k=1,n0 |
| 474 | 0 | ppoly_r_E(k,1) = u0(k) - 0.5 * ppoly_r_coefs(k,2) ! Left edge value of cell k |
| 475 | 0 | ppoly_r_E(k,2) = u0(k) + 0.5 * ppoly_r_coefs(k,2) ! Right edge value of cell k |
| 476 | 0 | ppoly_r_coefs(k,1) = ppoly_r_E(k,1) |
| 477 | enddo | |
| 478 | 0 | if ( CS%boundary_extrapolation ) & |
| 479 | 0 | call PLM_boundary_extrapolation( n0, h0, u0, ppoly_r_E, ppoly_r_coefs, h_neglect ) |
| 480 | 0 | iMethod = INTEGRATION_PLM |
| 481 | case ( REMAPPING_PPM_CW ) | |
| 482 | ! identical to REMAPPING_PPM_HYBGEN | |
| 483 | 0 | call edge_values_explicit_h4cw( n0, h0, u0, ppoly_r_E, h_neg_edge ) |
| 484 | 0 | call PPM_monotonicity( n0, u0, ppoly_r_E ) |
| 485 | 0 | call PPM_reconstruction( n0, h0, u0, ppoly_r_E, ppoly_r_coefs, h_neglect, answer_date=CS%answer_date ) |
| 486 | 0 | if ( CS%boundary_extrapolation ) then |
| 487 | 0 | call PPM_boundary_extrapolation( n0, h0, u0, ppoly_r_E, ppoly_r_coefs, h_neglect ) |
| 488 | endif | |
| 489 | 0 | iMethod = INTEGRATION_PPM |
| 490 | case ( REMAPPING_PPM_H4 ) | |
| 491 | 0 | call edge_values_explicit_h4( n0, h0, u0, ppoly_r_E, h_neg_edge, answer_date=CS%answer_date ) |
| 492 | 0 | call PPM_reconstruction( n0, h0, u0, ppoly_r_E, ppoly_r_coefs, h_neglect, answer_date=CS%answer_date ) |
| 493 | 0 | if ( CS%boundary_extrapolation ) then |
| 494 | 0 | call PPM_boundary_extrapolation( n0, h0, u0, ppoly_r_E, ppoly_r_coefs, h_neglect ) |
| 495 | endif | |
| 496 | 0 | iMethod = INTEGRATION_PPM |
| 497 | case ( REMAPPING_PPM_IH4 ) | |
| 498 | 303504 | call edge_values_implicit_h4( n0, h0, u0, ppoly_r_E, h_neg_edge, answer_date=CS%answer_date ) |
| 499 | 303504 | call PPM_reconstruction( n0, h0, u0, ppoly_r_E, ppoly_r_coefs, h_neglect, answer_date=CS%answer_date ) |
| 500 | 303504 | if ( CS%boundary_extrapolation ) then |
| 501 | 0 | call PPM_boundary_extrapolation( n0, h0, u0, ppoly_r_E, ppoly_r_coefs, h_neglect ) |
| 502 | endif | |
| 503 | 303504 | iMethod = INTEGRATION_PPM |
| 504 | case ( REMAPPING_PPM_HYBGEN ) | |
| 505 | 0 | call hybgen_PPM_coefs(u0, h0, ppoly_r_E, n0, 1, h_neglect) |
| 506 | 0 | call PPM_reconstruction( n0, h0, u0, ppoly_r_E, ppoly_r_coefs, h_neglect, answer_date=99991231 ) |
| 507 | 0 | if ( CS%boundary_extrapolation ) & |
| 508 | 0 | call PPM_boundary_extrapolation( n0, h0, u0, ppoly_r_E, ppoly_r_coefs, h_neglect ) |
| 509 | 0 | iMethod = INTEGRATION_PPM |
| 510 | case ( REMAPPING_WENO_HYBGEN ) | |
| 511 | 0 | call hybgen_weno_coefs(u0, h0, ppoly_r_E, n0, 1, h_neglect) |
| 512 | 0 | call PPM_reconstruction( n0, h0, u0, ppoly_r_E, ppoly_r_coefs, h_neglect, answer_date=99991231 ) |
| 513 | 0 | if ( CS%boundary_extrapolation ) & |
| 514 | 0 | call PPM_boundary_extrapolation( n0, h0, u0, ppoly_r_E, ppoly_r_coefs, h_neglect ) |
| 515 | 0 | iMethod = INTEGRATION_PPM |
| 516 | case ( REMAPPING_PQM_IH4IH3 ) | |
| 517 | 0 | call edge_values_implicit_h4( n0, h0, u0, ppoly_r_E, h_neg_edge, answer_date=CS%answer_date ) |
| 518 | 0 | 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 | 0 | answer_date=CS%answer_date ) |
| 521 | 0 | if ( CS%boundary_extrapolation ) then |
| 522 | call PQM_boundary_extrapolation_v1( n0, h0, u0, ppoly_r_E, ppoly_r_S, & | |
| 523 | 0 | ppoly_r_coefs, h_neglect ) |
| 524 | endif | |
| 525 | 0 | iMethod = INTEGRATION_PQM |
| 526 | case ( REMAPPING_PQM_IH6IH5 ) | |
| 527 | 0 | call edge_values_implicit_h6( n0, h0, u0, ppoly_r_E, h_neg_edge, answer_date=CS%answer_date ) |
| 528 | 0 | 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 | 0 | answer_date=CS%answer_date ) |
| 531 | 0 | if ( CS%boundary_extrapolation ) then |
| 532 | call PQM_boundary_extrapolation_v1( n0, h0, u0, ppoly_r_E, ppoly_r_S, & | |
| 533 | 0 | ppoly_r_coefs, h_neglect ) |
| 534 | endif | |
| 535 | 0 | iMethod = INTEGRATION_PQM |
| 536 | case ( REMAPPING_VIA_CLASS ) | |
| 537 | call MOM_error( FATAL, 'MOM_remapping, build_reconstructions_1d: '//& | |
| 538 | 0 | '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 | 626996 | 'The selected remapping method is invalid' ) |
| 542 | end select | |
| 543 | ||
| 544 | 626996 | if (present(PCM_cell)) then |
| 545 | ! Change the coefficients to those for the piecewise constant method in indicated cells. | |
| 546 | 14861496 | do k=1,n0 ; if (PCM_cell(k)) then |
| 547 | 0 | ppoly_r_coefs(k,1) = u0(k) |
| 548 | 0 | ppoly_r_E(k,1:2) = u0(k) |
| 549 | 0 | ppoly_r_S(k,1:2) = 0.0 |
| 550 | 0 | do n=2,CS%degree+1 ; ppoly_r_coefs(k,n) = 0.0 ; enddo |
| 551 | endif ; enddo | |
| 552 | endif | |
| 553 | ||
| 554 | 626996 | end subroutine build_reconstructions_1d |
| 555 | ||
| 556 | !> Checks that edge values and reconstructions satisfy bounds | |
| 557 | 0 | subroutine check_reconstructions_1d(n0, h0, u0, deg, boundary_extrapolation, & |
| 558 | 0 | 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 | 0 | problem_detected = .false. |
| 573 | 0 | do i0 = 1, n0 |
| 574 | 0 | u_l = u0(max(1,i0-1)) |
| 575 | 0 | u_c = u0(i0) |
| 576 | 0 | u_r = u0(min(n0,i0+1)) |
| 577 | 0 | if (i0 > 1 .or. .not. boundary_extrapolation) then |
| 578 | 0 | u_min = min(u_l, u_c) |
| 579 | 0 | u_max = max(u_l, u_c) |
| 580 | 0 | if (ppoly_r_E(i0,1) < u_min) then |
| 581 | 0 | write(0,'(a,I0,5(1x,a,1pe24.16))') 'Left edge undershoot at ',i0,'u(i0-1)=',u_l,'u(i0)=',u_c, & |
| 582 | 0 | 'edge=',ppoly_r_E(i0,1),'err=',ppoly_r_E(i0,1)-u_min |
| 583 | 0 | problem_detected = .true. |
| 584 | endif | |
| 585 | 0 | if (ppoly_r_E(i0,1) > u_max) then |
| 586 | 0 | write(0,'(a,I0,5(1x,a,1pe24.16))') 'Left edge overshoot at ',i0,'u(i0-1)=',u_l,'u(i0)=',u_c, & |
| 587 | 0 | 'edge=',ppoly_r_E(i0,1),'err=',ppoly_r_E(i0,1)-u_max |
| 588 | 0 | problem_detected = .true. |
| 589 | endif | |
| 590 | endif | |
| 591 | 0 | if (i0 < n0 .or. .not. boundary_extrapolation) then |
| 592 | 0 | u_min = min(u_c, u_r) |
| 593 | 0 | u_max = max(u_c, u_r) |
| 594 | 0 | if (ppoly_r_E(i0,2) < u_min) then |
| 595 | 0 | write(0,'(a,I0,5(1x,a,1pe24.16))') 'Right edge undershoot at ',i0,'u(i0)=',u_c,'u(i0+1)=',u_r, & |
| 596 | 0 | 'edge=',ppoly_r_E(i0,2),'err=',ppoly_r_E(i0,2)-u_min |
| 597 | 0 | problem_detected = .true. |
| 598 | endif | |
| 599 | 0 | if (ppoly_r_E(i0,2) > u_max) then |
| 600 | 0 | write(0,'(a,I0,5(1x,a,1pe24.16))') 'Right edge overshoot at ',i0,'u(i0)=',u_c,'u(i0+1)=',u_r, & |
| 601 | 0 | 'edge=',ppoly_r_E(i0,2),'err=',ppoly_r_E(i0,2)-u_max |
| 602 | 0 | problem_detected = .true. |
| 603 | endif | |
| 604 | endif | |
| 605 | 0 | if (i0 > 1) then |
| 606 | 0 | if ( (u_c-u_l)*(ppoly_r_E(i0,1)-ppoly_r_E(i0-1,2)) < 0.) then |
| 607 | 0 | write(0,'(a,I0,5(1x,a,1pe24.16))') 'Non-monotonic edges at',i0,'u(i0-1)=',u_l,'u(i0)=',u_c, & |
| 608 | 0 | 'right edge=',ppoly_r_E(i0-1,2),'left edge=',ppoly_r_E(i0,1) |
| 609 | 0 | 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 | 0 | problem_detected = .true. |
| 611 | endif | |
| 612 | endif | |
| 613 | 0 | if (problem_detected) then |
| 614 | 0 | write(0,'(a,1p9e24.16)') 'Polynomial coeffs:',ppoly_r_coefs(i0,:) |
| 615 | 0 | write(0,'(3(a,1pe24.16,1x))') 'u_l=',u_l,'u_c=',u_c,'u_r=',u_r |
| 616 | 0 | write(0,'(a4,10a24)') 'i0','h0(i0)','u0(i0)','left edge','right edge','Polynomial coefficients' |
| 617 | 0 | do n = 1, n0 |
| 618 | 0 | 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 | 0 | 'Edge values or polynomial coefficients were inconsistent!') |
| 622 | endif | |
| 623 | enddo | |
| 624 | ||
| 625 | 0 | end 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. | |
| 646 | 626996 | subroutine intersect_src_tgt_grids( n0, h0, n1, h1, h_sub, h0_eff, & |
| 647 | 626996 | 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 | 626996 | h0_supply = h0(1) |
| 677 | 626996 | h1_supply = h1(1) |
| 678 | 626996 | src_has_volume = .true. |
| 679 | 626996 | tgt_has_volume = .true. |
| 680 | 626996 | i0 = 1 ; i1 = 1 |
| 681 | 626996 | i_start0 = 1 ; i_start1 = 1 |
| 682 | 626996 | i_max = 1 |
| 683 | 626996 | dh_max = 0. |
| 684 | 626996 | dh0_eff = 0. |
| 685 | ||
| 686 | ! First sub-cell is always vanished | |
| 687 | 626996 | h_sub(1) = 0. |
| 688 | 626996 | isrc_start(1) = 1 |
| 689 | 626996 | isrc_end(1) = 1 |
| 690 | 626996 | isrc_max(1) = 1 |
| 691 | 626996 | isub_src(1) = 1 |
| 692 | ||
| 693 | ! Loop over each sub-cell to calculate intersections with source and target grids | |
| 694 | 81929228 | 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 | 81302232 | 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 | 81302232 | 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 | 81302232 | isub_src(i_sub) = i0 |
| 709 | 81302232 | h_sub(i_sub) = dh |
| 710 | ||
| 711 | ! For recording the largest sub-cell within a source cell. | |
| 712 | 81302232 | if (dh >= dh_max) then |
| 713 | 69540560 | i_max = i_sub |
| 714 | 69540560 | 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 | 81929228 | 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 | 46550704 | 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 | 46550704 | isrc_start(i0) = i_start0 |
| 726 | 46550704 | isrc_end(i0) = i_sub |
| 727 | 46550704 | i_start0 = i_sub + 1 |
| 728 | ! Record the sub-cell that is the largest fraction of the source cell. | |
| 729 | 46550704 | isrc_max(i0) = i_max |
| 730 | 46550704 | i_max = i_sub + 1 |
| 731 | 46550704 | dh_max = 0. |
| 732 | ! Record the source cell thickness found by summing the sub-cell thicknesses. | |
| 733 | 46550704 | h0_eff(i0) = dh0_eff |
| 734 | ! Move the source index. | |
| 735 | 46550704 | if (i0 < n0) then |
| 736 | 46117255 | i0 = i0 + 1 |
| 737 | 46117255 | h0_supply = h0(i0) |
| 738 | 46117255 | dh0_eff = 0. |
| 739 | else | |
| 740 | 433449 | h0_supply = 0. |
| 741 | 433449 | src_has_volume = .false. |
| 742 | endif | |
| 743 | 34751528 | 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 | 33974163 | 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 | 33974163 | itgt_start(i1) = i_start1 |
| 750 | 33974163 | itgt_end(i1) = i_sub |
| 751 | 33974163 | i_start1 = i_sub + 1 |
| 752 | ! Move the target index. | |
| 753 | 33974163 | if (i1 < n1) then |
| 754 | 33457720 | i1 = i1 + 1 |
| 755 | 33457720 | h1_supply = h1(i1) |
| 756 | else | |
| 757 | 516443 | h1_supply = 0. |
| 758 | 516443 | tgt_has_volume = .false. |
| 759 | endif | |
| 760 | 777365 | elseif (src_has_volume) then |
| 761 | ! We ran out of target volume but still have source cells to consume | |
| 762 | 473996 | h_sub(i_sub) = h0_supply |
| 763 | ! Record the sub-cell start/end index that span the source cell i0. | |
| 764 | 473996 | isrc_start(i0) = i_start0 |
| 765 | 473996 | isrc_end(i0) = i_sub |
| 766 | 473996 | i_start0 = i_sub + 1 |
| 767 | ! Record the sub-cell that is the largest fraction of the source cell. | |
| 768 | 473996 | isrc_max(i0) = i_max |
| 769 | 473996 | i_max = i_sub + 1 |
| 770 | 473996 | dh_max = 0. |
| 771 | ! Record the source cell thickness found by summing the sub-cell thicknesses. | |
| 772 | 473996 | h0_eff(i0) = dh0_eff |
| 773 | 473996 | if (i0 < n0) then |
| 774 | 280449 | i0 = i0 + 1 |
| 775 | 280449 | h0_supply = h0(i0) |
| 776 | 280449 | dh0_eff = 0. |
| 777 | else | |
| 778 | 193547 | h0_supply = 0. |
| 779 | 193547 | src_has_volume = .false. |
| 780 | endif | |
| 781 | 303369 | elseif (tgt_has_volume) then |
| 782 | ! We ran out of source volume but still have target cells to consume | |
| 783 | 303369 | h_sub(i_sub) = h1_supply |
| 784 | ! Record the sub-cell start/end index that span the target cell i1. | |
| 785 | 303369 | itgt_start(i1) = i_start1 |
| 786 | 303369 | itgt_end(i1) = i_sub |
| 787 | 303369 | i_start1 = i_sub + 1 |
| 788 | ! Move the target index. | |
| 789 | 303369 | if (i1 < n1) then |
| 790 | 192816 | i1 = i1 + 1 |
| 791 | 192816 | h1_supply = h1(i1) |
| 792 | else | |
| 793 | 110553 | h1_supply = 0. |
| 794 | 110553 | tgt_has_volume = .false. |
| 795 | endif | |
| 796 | else | |
| 797 | 0 | stop 'intersect_src_tgt_grids: THIS SHOULD NEVER HAPPEN!' |
| 798 | endif | |
| 799 | ||
| 800 | enddo | |
| 801 | ||
| 802 | 626996 | end 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. | |
| 849 | 626996 | subroutine remap_src_to_sub_grid_om4(n0, h0, u0, ppoly0_E, ppoly0_coefs, n1, h_sub, & |
| 850 | 1253992 | h0_eff, isrc_start, isrc_end, isrc_max, isub_src, & |
| 851 | 626996 | 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 | 1253992 | 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 | 626996 | i0_last_thick_cell = 0 |
| 885 | 47651696 | do i0 = 1, n0 |
| 886 | 47024700 | u0_min(i0) = min(ppoly0_E(i0,1), ppoly0_E(i0,2)) |
| 887 | 47024700 | u0_max(i0) = max(ppoly0_E(i0,1), ppoly0_E(i0,2)) |
| 888 | 47651696 | 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 | 626996 | xa = 0. |
| 895 | 626996 | dh0_eff = 0. |
| 896 | 626996 | uh_sub(1) = 0. |
| 897 | 626996 | u_sub(1) = ppoly0_E(1,1) |
| 898 | 626996 | u02_err = 0. |
| 899 | 81302232 | do i_sub = 2, n0+n1 |
| 900 | ||
| 901 | ! Sub-cell thickness from loop above | |
| 902 | 80675236 | dh = h_sub(i_sub) |
| 903 | ||
| 904 | ! Source cell | |
| 905 | 80675236 | 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 | 80675236 | dh0_eff = dh0_eff + dh ! Cumulative thickness within the source cell |
| 911 | 80675236 | if (h0_eff(i0)>0.) then |
| 912 | 80404435 | xb = dh0_eff / h0_eff(i0) ! This expression yields xa <= xb <= 1.0 |
| 913 | 80404435 | xb = min(1., xb) ! This is only needed when the total target column is wider than the source column |
| 914 | 80404435 | u_sub(i_sub) = average_value_ppoly( n0, u0, ppoly0_E, ppoly0_coefs, method, i0, xa, xb) |
| 915 | else ! Vanished cell | |
| 916 | 270801 | xb = 1. |
| 917 | 270801 | u_sub(i_sub) = u0(i0) |
| 918 | endif | |
| 919 | 80675236 | 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 | 0 | u_orig = u_sub(i_sub) |
| 924 | 0 | u_sub(i_sub) = max( u_sub(i_sub), u0_min(i0) ) |
| 925 | 0 | u_sub(i_sub) = min( u_sub(i_sub), u0_max(i0) ) |
| 926 | 0 | u02_err = u02_err + dh*abs( u_sub(i_sub) - u_orig ) |
| 927 | endif | |
| 928 | 80675236 | uh_sub(i_sub) = dh * u_sub(i_sub) |
| 929 | ||
| 930 | 81302232 | 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 | 46397704 | dh0_eff = 0. |
| 933 | 46397704 | xa = 0. |
| 934 | else | |
| 935 | 34277532 | xa = xb ! Next integral will start at end of last |
| 936 | endif | |
| 937 | ||
| 938 | enddo | |
| 939 | 626996 | u_sub(n0+n1+1) = ppoly0_E(n0,2) ! This value is only needed when total target column |
| 940 | 626996 | 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 | 47651696 | do i0 = 1, i0_last_thick_cell |
| 949 | 47024700 | i_max = isrc_max(i0) |
| 950 | 47024700 | dh_max = h_sub(i_max) |
| 951 | 47651696 | 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 | 47024700 | duh = 0. |
| 954 | 127249902 | do i_sub = isrc_start(i0), isrc_end(i0) |
| 955 | 127249902 | if (i_sub /= i_max) duh = duh + uh_sub(i_sub) |
| 956 | enddo | |
| 957 | 47024700 | uh_sub(i_max) = u0(i0)*h0(i0) - duh |
| 958 | 47024700 | u02_err = u02_err + max( abs(uh_sub(i_max)), abs(u0(i0)*h0(i0)), abs(duh) ) |
| 959 | endif | |
| 960 | enddo | |
| 961 | endif | |
| 962 | ||
| 963 | 626996 | end subroutine remap_src_to_sub_grid_om4 |
| 964 | ||
| 965 | !> Remaps column of n0 values u0 on grid h0 to subgrid h_sub | |
| 966 | 0 | subroutine remap_src_to_sub_grid(n0, h0, u0, ppoly0_E, ppoly0_coefs, n1, h_sub, & |
| 967 | 0 | isrc_start, isrc_end, isrc_max, isub_src, & |
| 968 | 0 | 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 | 0 | 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 | 0 | i0_last_thick_cell = 0 |
| 1001 | 0 | do i0 = 1, n0 |
| 1002 | 0 | u0_min(i0) = min(ppoly0_E(i0,1), ppoly0_E(i0,2)) |
| 1003 | 0 | u0_max(i0) = max(ppoly0_E(i0,1), ppoly0_E(i0,2)) |
| 1004 | 0 | 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 | 0 | xa = 0. |
| 1011 | 0 | dh0_eff = 0. |
| 1012 | 0 | u02_err = 0. |
| 1013 | 0 | do i_sub = 1, n0+n1 |
| 1014 | ||
| 1015 | ! Sub-cell thickness from loop above | |
| 1016 | 0 | dh = h_sub(i_sub) |
| 1017 | ||
| 1018 | ! Source cell | |
| 1019 | 0 | 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 | 0 | dh0_eff = dh0_eff + dh ! Cumulative thickness within the source cell |
| 1025 | 0 | if (h0(i0)>0.) then |
| 1026 | 0 | xb = dh0_eff / h0(i0) ! This expression yields xa <= xb <= 1.0 |
| 1027 | 0 | xb = min(1., xb) ! This is only needed when the total target column is wider than the source column |
| 1028 | 0 | u_sub(i_sub) = average_value_ppoly( n0, u0, ppoly0_E, ppoly0_coefs, method, i0, xa, xb) |
| 1029 | else ! Vanished cell | |
| 1030 | 0 | xb = 1. |
| 1031 | 0 | u_sub(i_sub) = u0(i0) |
| 1032 | endif | |
| 1033 | 0 | 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 | 0 | u_orig = u_sub(i_sub) |
| 1038 | 0 | u_sub(i_sub) = max( u_sub(i_sub), u0_min(i0) ) |
| 1039 | 0 | u_sub(i_sub) = min( u_sub(i_sub), u0_max(i0) ) |
| 1040 | 0 | u02_err = u02_err + dh*abs( u_sub(i_sub) - u_orig ) |
| 1041 | endif | |
| 1042 | 0 | uh_sub(i_sub) = dh * u_sub(i_sub) |
| 1043 | ||
| 1044 | 0 | 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 | 0 | dh0_eff = 0. |
| 1047 | 0 | xa = 0. |
| 1048 | else | |
| 1049 | 0 | xa = xb ! Next integral will start at end of last |
| 1050 | endif | |
| 1051 | ||
| 1052 | enddo | |
| 1053 | 0 | i_sub = n0+n1+1 |
| 1054 | ! Sub-cell thickness from loop above | |
| 1055 | 0 | dh = h_sub(i_sub) |
| 1056 | ! Source cell | |
| 1057 | 0 | 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 | 0 | dh0_eff = dh0_eff + dh ! Cumulative thickness within the source cell |
| 1063 | 0 | if (h0(i0)>0.) then |
| 1064 | 0 | xb = dh0_eff / h0(i0) ! This expression yields xa <= xb <= 1.0 |
| 1065 | 0 | xb = min(1., xb) ! This is only needed when the total target column is wider than the source column |
| 1066 | 0 | u_sub(i_sub) = average_value_ppoly( n0, u0, ppoly0_E, ppoly0_coefs, method, i0, xa, xb) |
| 1067 | else ! Vanished cell | |
| 1068 | 0 | xb = 1. |
| 1069 | 0 | u_sub(i_sub) = u0(i0) |
| 1070 | endif | |
| 1071 | 0 | 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 | 0 | u_orig = u_sub(i_sub) |
| 1076 | 0 | u_sub(i_sub) = max( u_sub(i_sub), u0_min(i0) ) |
| 1077 | 0 | u_sub(i_sub) = min( u_sub(i_sub), u0_max(i0) ) |
| 1078 | 0 | u02_err = u02_err + dh*abs( u_sub(i_sub) - u_orig ) |
| 1079 | endif | |
| 1080 | 0 | 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 | 0 | do i0 = 1, i0_last_thick_cell |
| 1089 | 0 | i_max = isrc_max(i0) |
| 1090 | 0 | dh_max = h_sub(i_max) |
| 1091 | 0 | 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 | 0 | duh = 0. |
| 1094 | 0 | do i_sub = isrc_start(i0), isrc_end(i0) |
| 1095 | 0 | if (i_sub /= i_max) duh = duh + uh_sub(i_sub) |
| 1096 | enddo | |
| 1097 | 0 | uh_sub(i_max) = u0(i0)*h0(i0) - duh |
| 1098 | 0 | u02_err = u02_err + max( abs(uh_sub(i_max)), abs(u0(i0)*h0(i0)), abs(duh) ) |
| 1099 | endif | |
| 1100 | enddo | |
| 1101 | endif | |
| 1102 | ||
| 1103 | 0 | end 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 | |
| 1107 | 626996 | subroutine remap_sub_to_tgt_grid_om4(n0, n1, h1, h_sub, u_sub, uh_sub, & |
| 1108 | 626996 | 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 | 626996 | u1min = 0. ! Not necessary, but avoids an overzealous compiler ... |
| 1129 | 626996 | 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 | 626996 | uh_err = 0. |
| 1135 | 34904528 | do i1 = 1, n1 |
| 1136 | 34904528 | if (h1(i1) > 0.) then |
| 1137 | 32029152 | duh = 0. ; dh = 0. |
| 1138 | 32029152 | i_sub = itgt_start(i1) |
| 1139 | 32029152 | if (force_bounds_in_target) then |
| 1140 | 32029152 | u1min = u_sub(i_sub) |
| 1141 | 32029152 | u1max = u_sub(i_sub) |
| 1142 | endif | |
| 1143 | 111236004 | do i_sub = itgt_start(i1), itgt_end(i1) |
| 1144 | 79206852 | if (force_bounds_in_target) then |
| 1145 | 79206852 | u1min = min(u1min, u_sub(i_sub)) |
| 1146 | 79206852 | u1max = max(u1max, u_sub(i_sub)) |
| 1147 | endif | |
| 1148 | 79206852 | dh = dh + h_sub(i_sub) |
| 1149 | 79206852 | duh = duh + uh_sub(i_sub) |
| 1150 | ! This accumulates the contribution to the error bound for the sum of u*h | |
| 1151 | 111236004 | uh_err = uh_err + max(abs(duh),abs(uh_sub(i_sub)))*epsilon(duh) |
| 1152 | enddo | |
| 1153 | 32029152 | u1(i1) = duh / dh |
| 1154 | ! This is the contribution from the division to the error bound for the sum of u*h | |
| 1155 | 32029152 | uh_err = uh_err + abs(duh)*epsilon(duh) |
| 1156 | 32029152 | if (force_bounds_in_target) then |
| 1157 | 32029152 | u_orig = u1(i1) |
| 1158 | 32029152 | u1(i1) = max(u1min, min(u1max, u1(i1))) |
| 1159 | ! Adjusting to be bounded contributes to the error for the sum of u*h | |
| 1160 | 32029152 | uh_err = uh_err + dh*abs( u1(i1)-u_orig ) |
| 1161 | endif | |
| 1162 | else | |
| 1163 | 2248380 | u1(i1) = u_sub(itgt_start(i1)) |
| 1164 | endif | |
| 1165 | enddo | |
| 1166 | ||
| 1167 | 626996 | end 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 | |
| 1170 | 0 | subroutine remap_sub_to_tgt_grid(n0, n1, h1, h_sub, u_sub, uh_sub, & |
| 1171 | 0 | itgt_start, itgt_end, force_bounds_in_target, & |
| 1172 | 0 | 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 | 0 | u1min = 0. ! Not necessary, but avoids an overzealous compiler ... |
| 1197 | 0 | u1max = 0. ! ... warning about uninitialized variables |
| 1198 | 0 | u_ref = 0. ! An offset of 0. should do no harm |
| 1199 | 0 | 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 | 0 | uh_err = 0. |
| 1205 | 0 | do i1 = 1, n1 |
| 1206 | 0 | if (h1(i1) > 0.) then |
| 1207 | 0 | duh = 0. ; dh = 0. |
| 1208 | 0 | i_sub = itgt_start(i1) |
| 1209 | 0 | if (force_bounds_in_target) then |
| 1210 | 0 | u1min = u_sub(i_sub) |
| 1211 | 0 | u1max = u_sub(i_sub) |
| 1212 | endif | |
| 1213 | 0 | if (offset_summation) then |
| 1214 | 0 | u_ref = 0. ! An offset of 0. should do no harm |
| 1215 | 0 | h_max = 0. |
| 1216 | 0 | do i_sub = itgt_start(i1), itgt_end(i1) |
| 1217 | 0 | if (h_sub(i_sub) > h_max) then |
| 1218 | 0 | u_ref = u_sub(i_sub) |
| 1219 | 0 | h_max = h_sub(i_sub) |
| 1220 | endif | |
| 1221 | enddo | |
| 1222 | endif | |
| 1223 | 0 | do i_sub = itgt_start(i1), itgt_end(i1) |
| 1224 | 0 | if (force_bounds_in_target .or. better_force_bounds_in_target .and. h_sub(i_sub)>0.) then |
| 1225 | 0 | u1min = min(u1min, u_sub(i_sub)) |
| 1226 | 0 | u1max = max(u1max, u_sub(i_sub)) |
| 1227 | endif | |
| 1228 | 0 | dh = dh + h_sub(i_sub) |
| 1229 | ! Ideally u_ref would be already be substracted in uh_sub | |
| 1230 | 0 | 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 | 0 | uh_err = uh_err + max(abs(duh),abs(uh_sub(i_sub)))*epsilon(duh) |
| 1233 | enddo | |
| 1234 | 0 | 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 | 0 | uh_err = uh_err + abs(duh)*epsilon(duh) |
| 1237 | 0 | if (force_bounds_in_target) then |
| 1238 | 0 | u_orig = u1(i1) |
| 1239 | 0 | u1(i1) = max(u1min, min(u1max, u1(i1))) |
| 1240 | ! Adjusting to be bounded contributes to the error for the sum of u*h | |
| 1241 | 0 | uh_err = uh_err + dh*abs( u1(i1)-u_orig ) |
| 1242 | endif | |
| 1243 | else | |
| 1244 | 0 | u1(i1) = u_sub(itgt_start(i1)) |
| 1245 | endif | |
| 1246 | enddo | |
| 1247 | ||
| 1248 | 0 | end subroutine remap_sub_to_tgt_grid |
| 1249 | ||
| 1250 | !> Linearly interpolate interface data, u_src, from grid h_src to a grid h_dest | |
| 1251 | 61392 | subroutine 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 | 122784 | real :: frac_pos(ndest+1) ! Fractional position of the destination interface |
| 1265 | ! within the source layer [nondim], 0 <= frac_pos <= 1. | |
| 1266 | 122784 | 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 | 61392 | ks = 0 |
| 1271 | 61392 | dh = 0. |
| 1272 | 61392 | 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 | 2148720 | do k_dest=1,ndest+1 |
| 1277 | 6567816 | do while (dh<=x_dest .and. ks<nsrc) |
| 1278 | ! Move positions pointers forward until the interval 0 .. dh spans x_dest. | |
| 1279 | 4480488 | x_dest = x_dest - dh |
| 1280 | 4480488 | ks = ks + 1 |
| 1281 | 4480488 | dh = h_src(ks) ! Thickness of layer ks |
| 1282 | enddo | |
| 1283 | 2087328 | k_src(k_dest) = ks |
| 1284 | ||
| 1285 | 2087328 | if (dh>0.) then |
| 1286 | 2087328 | 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 | 0 | frac_pos(k_dest) = 0.5 |
| 1289 | endif | |
| 1290 | ||
| 1291 | 2148720 | if (k_dest <= ndest) then |
| 1292 | 2025936 | x_dest = x_dest + h_dest(k_dest) ! Position of interface k_dest+1 |
| 1293 | endif | |
| 1294 | enddo | |
| 1295 | ||
| 1296 | 2148720 | do k_dest=1,ndest+1 |
| 1297 | ! Linear interpolation between surrounding edge values. | |
| 1298 | 2087328 | ks = k_src(k_dest) |
| 1299 | 2148720 | 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 | 61392 | 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 | 61392 | do k_dest=1,ndest |
| 1307 | 61392 | if (h_dest(k_dest) > 0.) exit |
| 1308 | 0 | u_dest(k_dest) = 0.0 |
| 1309 | enddo | |
| 1310 | ||
| 1311 | ! Mask interfaces below vanished layers at the bottom | |
| 1312 | 509772 | do k_dest=ndest,1,-1 |
| 1313 | 509772 | if (h_dest(k_dest) > 0.) exit |
| 1314 | 448380 | u_dest(k_dest+1) = 0.0 |
| 1315 | enddo | |
| 1316 | endif | |
| 1317 | ||
| 1318 | 61392 | end 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 | |
| 1321 | 0 | subroutine 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 | 0 | uh_dest(:) = 0.0 |
| 1336 | ||
| 1337 | 0 | k_src = 0 |
| 1338 | 0 | k_dest = 0 |
| 1339 | 0 | h_dest_rem = 0. |
| 1340 | 0 | h_src_rem = 0. |
| 1341 | 0 | uh_src_rem = 0. |
| 1342 | 0 | src_ran_out = .false. |
| 1343 | ||
| 1344 | 0 | do while(.true.) |
| 1345 | 0 | if (h_src_rem==0. .and. k_src<nsrc) then |
| 1346 | ! Supply is empty so move to the next source cell | |
| 1347 | 0 | k_src = k_src + 1 |
| 1348 | 0 | h_src_rem = h_src(k_src) |
| 1349 | 0 | uh_src_rem = uh_src(k_src) |
| 1350 | 0 | if (h_src_rem==0.) cycle |
| 1351 | endif | |
| 1352 | 0 | if (h_dest_rem==0. .and. k_dest<ndest) then |
| 1353 | ! Sink has no capacity so move to the next destination cell | |
| 1354 | 0 | k_dest = k_dest + 1 |
| 1355 | 0 | h_dest_rem = h_dest(k_dest) |
| 1356 | 0 | uh_dest(k_dest) = 0. |
| 1357 | 0 | if (h_dest_rem==0.) cycle |
| 1358 | endif | |
| 1359 | 0 | if (k_src==nsrc .and. h_src_rem==0.) then |
| 1360 | 0 | if (src_ran_out) exit ! This is the second time implying there is no more src |
| 1361 | 0 | src_ran_out = .true. |
| 1362 | 0 | cycle |
| 1363 | endif | |
| 1364 | 0 | duh = 0. |
| 1365 | 0 | if (h_src_rem<h_dest_rem) then |
| 1366 | ! The source cell is fully within the destination cell | |
| 1367 | 0 | dh = h_src_rem |
| 1368 | 0 | if (dh>0.) duh = uh_src_rem |
| 1369 | 0 | h_src_rem = 0. |
| 1370 | 0 | uh_src_rem = 0. |
| 1371 | 0 | h_dest_rem = max(0., h_dest_rem - dh) |
| 1372 | 0 | elseif (h_src_rem>h_dest_rem) then |
| 1373 | ! Only part of the source cell can be used up | |
| 1374 | 0 | dh = h_dest_rem |
| 1375 | 0 | duh = (dh / h_src_rem) * uh_src_rem |
| 1376 | 0 | h_src_rem = max(0., h_src_rem - dh) |
| 1377 | 0 | uh_src_rem = uh_src_rem - duh |
| 1378 | 0 | h_dest_rem = 0. |
| 1379 | else ! h_src_rem==h_dest_rem | |
| 1380 | ! The source cell exactly fits the destination cell | |
| 1381 | 0 | duh = uh_src_rem |
| 1382 | 0 | h_src_rem = 0. |
| 1383 | 0 | uh_src_rem = 0. |
| 1384 | 0 | h_dest_rem = 0. |
| 1385 | endif | |
| 1386 | 0 | uh_dest(k_dest) = uh_dest(k_dest) + duh |
| 1387 | 0 | if (k_dest==ndest .and. (k_src==nsrc .or. h_dest_rem==0.)) exit |
| 1388 | enddo | |
| 1389 | ||
| 1390 | 0 | end 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. | |
| 1395 | 80404435 | real 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 | 80404435 | u_ave = 0. ! Avoids warnings about "potentially unset values"; u_ave is always calculated for legitimate schemes |
| 1414 | 80404435 | if (xb > xa) then |
| 1415 | 78107036 | select case ( method ) |
| 1416 | case ( INTEGRATION_PCM ) | |
| 1417 | 0 | u_ave = u0(i0) |
| 1418 | case ( INTEGRATION_PLM ) | |
| 1419 | u_ave = ( & | |
| 1420 | ppoly0_coefs(i0,1) & | |
| 1421 | 48160937 | + ppoly0_coefs(i0,2) * 0.5 * ( xb + xa ) ) |
| 1422 | case ( INTEGRATION_PPM ) | |
| 1423 | 29946099 | mx = 0.5 * ( xa + xb ) |
| 1424 | 29946099 | a_L = ppoly0_E(i0, 1) |
| 1425 | 29946099 | a_R = ppoly0_E(i0, 2) |
| 1426 | 29946099 | u_c = u0(i0) |
| 1427 | 29946099 | a_c = 0.5 * ( ( u_c - a_L ) + ( u_c - a_R ) ) ! a_6 / 6 |
| 1428 | 59892198 | if (mx<0.5) then |
| 1429 | ! This integration of the PPM reconstruction is expressed in distances from the left edge | |
| 1430 | 7463748 | xa2b2ab = (xa*xa+xb*xb)+xa*xb |
| 1431 | u_ave = a_L + ( ( a_R - a_L ) * mx & | |
| 1432 | 7463748 | + 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 | 22482351 | Ya = 1. - xa |
| 1436 | 22482351 | Yb = 1. - xb |
| 1437 | 22482351 | my = 0.5 * ( Ya + Yb ) |
| 1438 | 22482351 | Ya2b2ab = (Ya*Ya+Yb*Yb)+Ya*Yb |
| 1439 | u_ave = a_R + ( ( a_L - a_R ) * my & | |
| 1440 | 22482351 | + a_c * ( 3. * ( Yb + Ya ) - 2.*Ya2b2ab ) ) |
| 1441 | endif | |
| 1442 | case ( INTEGRATION_PQM ) | |
| 1443 | 0 | xa_2 = xa*xa |
| 1444 | 0 | xb_2 = xb*xb |
| 1445 | 0 | xa2pxb2 = xa_2 + xb_2 |
| 1446 | 0 | 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 | 0 | + ppoly0_coefs(i0,5) * 0.2 * ( ( xb*xb_2 + xa*xa_2 ) * xapxb + xa_2*xb_2 ) ) ) ) ) |
| 1453 | case default | |
| 1454 | 78107036 | call MOM_error( FATAL,'The selected integration method is invalid' ) |
| 1455 | end select | |
| 1456 | else ! dh == 0. | |
| 1457 | 2297399 | select case ( method ) |
| 1458 | case ( INTEGRATION_PCM ) | |
| 1459 | 0 | u_ave = ppoly0_coefs(i0,1) |
| 1460 | case ( INTEGRATION_PLM ) | |
| 1461 | !u_ave = ppoly0_coefs(i0,1) & | |
| 1462 | ! + xa * ppoly0_coefs(i0,2) | |
| 1463 | 39371 | a_L = ppoly0_E(i0, 1) |
| 1464 | 39371 | a_R = ppoly0_E(i0, 2) |
| 1465 | 39371 | Ya = 1. - xa |
| 1466 | 78742 | if (xa < 0.5) then |
| 1467 | 39371 | u_ave = a_L + xa * ( a_R - a_L ) |
| 1468 | else | |
| 1469 | 0 | 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 | 2258028 | a_L = ppoly0_E(i0, 1) |
| 1476 | 2258028 | a_R = ppoly0_E(i0, 2) |
| 1477 | 2258028 | u_c = u0(i0) |
| 1478 | 2258028 | a_c = 3. * ( ( u_c - a_L ) + ( u_c - a_R ) ) ! a_6 |
| 1479 | 2258028 | Ya = 1. - xa |
| 1480 | 4516056 | if (xa < 0.5) then |
| 1481 | 0 | u_ave = a_L + xa * ( ( a_R - a_L ) + a_c * Ya ) |
| 1482 | else | |
| 1483 | 2258028 | 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 | 0 | + xa * ppoly0_coefs(i0,5) ) ) ) |
| 1491 | case default | |
| 1492 | 0 | u_ave = 0. |
| 1493 | 2297399 | call MOM_error( FATAL,'The selected integration method is invalid' ) |
| 1494 | end select | |
| 1495 | endif | |
| 1496 | 80404435 | average_value_ppoly = u_ave |
| 1497 | ||
| 1498 | 80404435 | end 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. | |
| 1502 | 0 | subroutine check_remapped_values(n0, h0, u0, ppoly_r_E, deg, ppoly_r_coefs, & |
| 1503 | 0 | 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 | 0 | call measure_input_bounds( n0, h0, u0, ppoly_r_E, h0tot, h0err, u0tot, u0err, u0min, u0max ) |
| 1529 | 0 | call measure_output_bounds( n1, h1, u1, h1tot, h1err, u1tot, u1err, u1min, u1max ) |
| 1530 | ||
| 1531 | 0 | 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 | 0 | .or. (u1min<u0min .or. u1max>u0max) ) then |
| 1535 | 0 | write(0,*) 'iMethod = ',iMethod |
| 1536 | 0 | write(0,*) 'H: h0tot=',h0tot,'h1tot=',h1tot,'dh=',h1tot-h0tot,'h0err=',h0err,'h1err=',h1err |
| 1537 | 0 | if (abs(h1tot-h0tot)>h0err+h1err) & |
| 1538 | 0 | write(0,*) 'H non-conservation difference=',h1tot-h0tot,'allowed err=',h0err+h1err,' <-----!' |
| 1539 | 0 | write(0,*) 'UH: u0tot=',u0tot,'u1tot=',u1tot,'duh=',u1tot-u0tot,'u0err=',u0err,'u1err=',u1err,'uh_err=',uh_err |
| 1540 | 0 | if (abs(u1tot-u0tot)>(u0err+u1err)+uh_err) & |
| 1541 | 0 | write(0,*) 'U non-conservation difference=',u1tot-u0tot,'allowed err=',u0err+u1err+uh_err,' <-----!' |
| 1542 | 0 | write(0,*) 'U: u0min=',u0min,'u1min=',u1min |
| 1543 | 0 | if (u1min<u0min) write(0,*) 'U minimum overshoot=',u1min-u0min,' <-----!' |
| 1544 | 0 | write(0,*) 'U: u0max=',u0max,'u1max=',u1max |
| 1545 | 0 | if (u1max>u0max) write(0,*) 'U maximum overshoot=',u1max-u0max,' <-----!' |
| 1546 | 0 | write(0,'(a3,6a24)') 'k','h0','left edge','u0','right edge','h1','u1' |
| 1547 | 0 | do k = 1, max(n0,n1) |
| 1548 | 0 | if (k<=min(n0,n1)) then |
| 1549 | 0 | 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 | 0 | elseif (k>n0) then |
| 1551 | 0 | write(0,'(i3,96x,1p2e24.16)') k,h1(k),u1(k) |
| 1552 | else | |
| 1553 | 0 | 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 | 0 | write(0,'(a3,2a24)') 'k','u0','Polynomial coefficients' |
| 1557 | 0 | do k = 1, n0 |
| 1558 | 0 | write(0,'(i3,1p6e24.16)') k,u0(k),ppoly_r_coefs(k,:) |
| 1559 | enddo | |
| 1560 | call MOM_error( FATAL, 'MOM_remapping, '//trim(caller)//': '//& | |
| 1561 | 0 | 'Remapping result is inconsistent!' ) |
| 1562 | endif | |
| 1563 | ||
| 1564 | 0 | end subroutine check_remapped_values |
| 1565 | ||
| 1566 | !> Measure totals and bounds on source grid | |
| 1567 | 0 | subroutine 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 | 0 | eps = epsilon(h0(1)) |
| 1583 | 0 | h0tot = h0(1) |
| 1584 | 0 | h0err = 0. |
| 1585 | 0 | u0tot = h0(1) * u0(1) |
| 1586 | 0 | u0err = 0. |
| 1587 | 0 | u0min = min( edge_values(1,1), edge_values(1,2) ) |
| 1588 | 0 | u0max = max( edge_values(1,1), edge_values(1,2) ) |
| 1589 | 0 | do k = 2, n0 |
| 1590 | 0 | h0tot = h0tot + h0(k) |
| 1591 | 0 | h0err = h0err + eps * max(h0tot, h0(k)) |
| 1592 | 0 | u0tot = u0tot + h0(k) * u0(k) |
| 1593 | 0 | u0err = u0err + eps * max(abs(u0tot), abs(h0(k) * u0(k))) |
| 1594 | 0 | u0min = min( u0min, edge_values(k,1), edge_values(k,2) ) |
| 1595 | 0 | u0max = max( u0max, edge_values(k,1), edge_values(k,2) ) |
| 1596 | enddo | |
| 1597 | ||
| 1598 | 0 | end subroutine measure_input_bounds |
| 1599 | ||
| 1600 | !> Measure totals and bounds on destination grid | |
| 1601 | 0 | subroutine 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 | 0 | eps = epsilon(h1(1)) |
| 1616 | 0 | h1tot = h1(1) |
| 1617 | 0 | h1err = 0. |
| 1618 | 0 | u1tot = h1(1) * u1(1) |
| 1619 | 0 | u1err = 0. |
| 1620 | 0 | u1min = u1(1) |
| 1621 | 0 | u1max = u1(1) |
| 1622 | 0 | do k = 2, n1 |
| 1623 | 0 | h1tot = h1tot + h1(k) |
| 1624 | 0 | h1err = h1err + eps * max(h1tot, h1(k)) |
| 1625 | 0 | u1tot = u1tot + h1(k) * u1(k) |
| 1626 | 0 | u1err = u1err + eps * max(abs(u1tot), abs(h1(k) * u1(k))) |
| 1627 | 0 | u1min = min(u1min, u1(k)) |
| 1628 | 0 | u1max = max(u1max, u1(k)) |
| 1629 | enddo | |
| 1630 | ||
| 1631 | 0 | end subroutine measure_output_bounds |
| 1632 | ||
| 1633 | !> Calculates the change in interface positions based on h1 and h2 | |
| 1634 | 0 | subroutine 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 | 0 | x1 = 0. |
| 1645 | 0 | x2 = 0. |
| 1646 | 0 | dx(1) = 0. |
| 1647 | 0 | do K = 1, max(n1,n2) |
| 1648 | 0 | if (k <= n1) x1 = x1 + h1(k) ! Interface k+1, right of source cell k |
| 1649 | 0 | if (k <= n2) then |
| 1650 | 0 | x2 = x2 + h2(k) ! Interface k+1, right of target cell k |
| 1651 | 0 | dx(K+1) = x2 - x1 ! Change of interface k+1, target - source |
| 1652 | endif | |
| 1653 | enddo | |
| 1654 | ||
| 1655 | 0 | end subroutine dzFromH1H2 |
| 1656 | ||
| 1657 | !> Constructor for remapping control structure | |
| 1658 | 5 | subroutine 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 | 5 | h_neglect_edge=h_neglect_edge) |
| 1698 | ||
| 1699 | 5 | end 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. | |
| 1705 | 5 | subroutine 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 | 5 | degree = -99 |
| 1711 | 5 | if (associated(CS%reconstruction)) then |
| 1712 | ! We have a choice of being careless and allowing easy re-use (e.g. when testing) | |
| 1713 | 0 | CS%remapping_scheme = -911 |
| 1714 | 0 | call CS%reconstruction%destroy() |
| 1715 | 0 | 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 | 10 | select case ( uppercase(trim(string)) ) |
| 1721 | case ("PCM") | |
| 1722 | 0 | CS%remapping_scheme = REMAPPING_PCM |
| 1723 | 0 | degree = 0 |
| 1724 | case ("PLM") | |
| 1725 | 4 | CS%remapping_scheme = REMAPPING_PLM |
| 1726 | 4 | degree = 1 |
| 1727 | case ("PLM_HYBGEN") | |
| 1728 | 0 | CS%remapping_scheme = REMAPPING_PLM_HYBGEN |
| 1729 | 0 | degree = 1 |
| 1730 | case ("PPM_CW") | |
| 1731 | 0 | CS%remapping_scheme = REMAPPING_PPM_CW |
| 1732 | 0 | degree = 2 |
| 1733 | case ("PPM_H4") | |
| 1734 | 0 | CS%remapping_scheme = REMAPPING_PPM_H4 |
| 1735 | 0 | degree = 2 |
| 1736 | case ("PPM_IH4") | |
| 1737 | 1 | CS%remapping_scheme = REMAPPING_PPM_IH4 |
| 1738 | 1 | degree = 2 |
| 1739 | case ("PPM_HYBGEN") | |
| 1740 | 0 | CS%remapping_scheme = REMAPPING_PPM_HYBGEN |
| 1741 | 0 | degree = 2 |
| 1742 | case ("WENO_HYBGEN") | |
| 1743 | 0 | CS%remapping_scheme = REMAPPING_WENO_HYBGEN |
| 1744 | 0 | degree = 2 |
| 1745 | case ("PQM_IH4IH3") | |
| 1746 | 0 | CS%remapping_scheme = REMAPPING_PQM_IH4IH3 |
| 1747 | 0 | degree = 4 |
| 1748 | case ("PQM_IH6IH5") | |
| 1749 | 0 | CS%remapping_scheme = REMAPPING_PQM_IH6IH5 |
| 1750 | 0 | degree = 4 |
| 1751 | case ("C_PCM") | |
| 1752 | 0 | allocate( PCM :: CS%reconstruction ) |
| 1753 | 0 | CS%remapping_scheme = REMAPPING_VIA_CLASS |
| 1754 | case ("C_PLM_CW") | |
| 1755 | 0 | allocate( PLM_CW :: CS%reconstruction ) |
| 1756 | 0 | CS%remapping_scheme = REMAPPING_VIA_CLASS |
| 1757 | case ("C_PLM_HYBGEN") | |
| 1758 | 0 | allocate( PLM_hybgen :: CS%reconstruction ) |
| 1759 | 0 | CS%remapping_scheme = REMAPPING_VIA_CLASS |
| 1760 | case ("C_MPLM_WA") | |
| 1761 | 0 | allocate( MPLM_WA :: CS%reconstruction ) |
| 1762 | 0 | CS%remapping_scheme = REMAPPING_VIA_CLASS |
| 1763 | case ("C_EMPLM_WA") | |
| 1764 | 0 | allocate( EMPLM_WA :: CS%reconstruction ) |
| 1765 | 0 | CS%remapping_scheme = REMAPPING_VIA_CLASS |
| 1766 | case ("C_MPLM_WA_POLY") | |
| 1767 | 0 | allocate( MPLM_WA_poly :: CS%reconstruction ) |
| 1768 | 0 | CS%remapping_scheme = REMAPPING_VIA_CLASS |
| 1769 | case ("C_EMPLM_WA_POLY") | |
| 1770 | 0 | allocate( EMPLM_WA_poly :: CS%reconstruction ) |
| 1771 | 0 | CS%remapping_scheme = REMAPPING_VIA_CLASS |
| 1772 | case ("C_PLM_CWK") | |
| 1773 | 0 | allocate( PLM_CWK :: CS%reconstruction ) |
| 1774 | 0 | CS%remapping_scheme = REMAPPING_VIA_CLASS |
| 1775 | case ("C_MPLM_CWK") | |
| 1776 | 0 | allocate( MPLM_CWK :: CS%reconstruction ) |
| 1777 | 0 | CS%remapping_scheme = REMAPPING_VIA_CLASS |
| 1778 | case ("C_EMPLM_CWK") | |
| 1779 | 0 | allocate( EMPLM_CWK :: CS%reconstruction ) |
| 1780 | 0 | CS%remapping_scheme = REMAPPING_VIA_CLASS |
| 1781 | case ("C_PPM_CW") | |
| 1782 | 0 | allocate( PPM_CW :: CS%reconstruction ) |
| 1783 | 0 | CS%remapping_scheme = REMAPPING_VIA_CLASS |
| 1784 | case ("C_PPM_HYBGEN") | |
| 1785 | 0 | allocate( PPM_hybgen :: CS%reconstruction ) |
| 1786 | 0 | CS%remapping_scheme = REMAPPING_VIA_CLASS |
| 1787 | case ("C_PPM_CWK") | |
| 1788 | 0 | allocate( PPM_CWK :: CS%reconstruction ) |
| 1789 | 0 | CS%remapping_scheme = REMAPPING_VIA_CLASS |
| 1790 | case ("C_EPPM_CWK") | |
| 1791 | 0 | allocate( EPPM_CWK :: CS%reconstruction ) |
| 1792 | 0 | CS%remapping_scheme = REMAPPING_VIA_CLASS |
| 1793 | case ("C_PPM_H4_2019") | |
| 1794 | 0 | allocate( PPM_H4_2019 :: CS%reconstruction ) |
| 1795 | 0 | CS%remapping_scheme = REMAPPING_VIA_CLASS |
| 1796 | case ("C_PPM_H4_2018") | |
| 1797 | 0 | allocate( PPM_H4_2018 :: CS%reconstruction ) |
| 1798 | 0 | CS%remapping_scheme = REMAPPING_VIA_CLASS |
| 1799 | case ("C_PLM_WLS") | |
| 1800 | 0 | allocate( PLM_WLS :: CS%reconstruction ) |
| 1801 | 0 | CS%remapping_scheme = REMAPPING_VIA_CLASS |
| 1802 | case default | |
| 1803 | call MOM_error(FATAL, "setReconstructionType: "//& | |
| 1804 | 10 | "Unrecognized choice for REMAPPING_SCHEME ("//trim(string)//").") |
| 1805 | end select | |
| 1806 | ||
| 1807 | 5 | CS%degree = degree |
| 1808 | ||
| 1809 | 5 | end subroutine setReconstructionType |
| 1810 | ||
| 1811 | !> Destrcutor for remapping control structure | |
| 1812 | 1 | subroutine end_remapping(CS) |
| 1813 | type(remapping_CS), intent(inout) :: CS !< Remapping control structure | |
| 1814 | ||
| 1815 | 1 | CS%degree = 0 |
| 1816 | ||
| 1817 | 1 | end subroutine end_remapping |
| 1818 | ||
| 1819 | !> Test if interpolate_column() produces the wrong answer | |
| 1820 | 0 | subroutine 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 | 0 | real, dimension(ndest+1) :: u_dest ! Interpolated value at destination cell interfaces [A] |
| 1831 | ||
| 1832 | ! Interpolate from src to dest | |
| 1833 | 0 | call interpolate_column(nsrc, h_src, u_src, ndest, h_dest, u_dest, .true.) |
| 1834 | 0 | call test%real_arr(ndest, u_dest, u_true, msg) |
| 1835 | 0 | end subroutine test_interp |
| 1836 | ||
| 1837 | !> Test if reintegrate_column() produces the wrong answer | |
| 1838 | 0 | subroutine 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 | 0 | real, dimension(ndest) :: uh_dest ! Reintegrated value on destination cells [A H] |
| 1849 | ||
| 1850 | ! Interpolate from src to dest | |
| 1851 | 0 | call reintegrate_column(nsrc, h_src, uh_src, ndest, h_dest, uh_dest) |
| 1852 | 0 | call test%real_arr(ndest, uh_dest, uh_true, msg) |
| 1853 | ||
| 1854 | 0 | end subroutine test_reintegrate |
| 1855 | ||
| 1856 | !> Test class-based remapping for internal consistency on random data | |
| 1857 | 0 | subroutine 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 | 0 | real :: h0(n0) ! Source grid [H but really nondim] |
| 1866 | 0 | 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 | 0 | 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 | 0 | force_bounds_in_subcell=.false. ) |
| 1875 | ||
| 1876 | 0 | call random_seed(size=seed_size) |
| 1877 | 0 | allocate( seed(seed_Size) ) |
| 1878 | 0 | seed(:) = 102030405 ! Repeatable sequences |
| 1879 | 0 | call random_seed(put=seed) |
| 1880 | ||
| 1881 | 0 | error = .false. |
| 1882 | 0 | do iter = 1, niter |
| 1883 | 0 | call random_number( h0 ) ! In range 0-1 |
| 1884 | 0 | h0(:) = max(0., h0(:) - 0.05) ! Make 5% of values equal to zero |
| 1885 | 0 | call random_number( u0 ) ! In range 0-1 |
| 1886 | ||
| 1887 | 0 | call remapCS%reconstruction%reconstruct(h0, u0) |
| 1888 | 0 | if ( remapCS%reconstruction%check_reconstruction(h0, u0) ) then |
| 1889 | 0 | if ( .not. error ) then ! Only dump first error |
| 1890 | 0 | print *,'iter=',iter |
| 1891 | 0 | print *,'h0',h0 |
| 1892 | 0 | print *,'u0',u0 |
| 1893 | endif | |
| 1894 | 0 | error = .true. |
| 1895 | endif | |
| 1896 | ||
| 1897 | enddo | |
| 1898 | ||
| 1899 | 0 | write(label,'(I0)') niter |
| 1900 | 0 | call test%test( error, trim(label)//' consistency tests of '//scheme ) |
| 1901 | ||
| 1902 | 0 | call remapCS%reconstruction%destroy() |
| 1903 | ||
| 1904 | 0 | end subroutine test_recon_consistency |
| 1905 | ||
| 1906 | !> Test that remapping a uniform field remains uniform | |
| 1907 | 0 | subroutine 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 | 0 | real :: h0(n0), h1(n0) ! Source and target grids [H but really nondim] |
| 1916 | 0 | 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 | 0 | 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 | 0 | om4_remap_via_sub_cells=.false.) |
| 1929 | ||
| 1930 | 0 | call random_seed(size=seed_size) |
| 1931 | 0 | allocate( seed(seed_Size) ) |
| 1932 | 0 | seed(:) = 102030405 ! Repeatable sequences |
| 1933 | 0 | call random_seed(put=seed) |
| 1934 | ||
| 1935 | 0 | error = .false. |
| 1936 | 0 | do iter = 1, niter |
| 1937 | 0 | call random_number( h0 ) ! In range 0-1 |
| 1938 | 0 | h0(:) = max(0., h0(:) - 0.05) ! Make 5% of values equal to zero |
| 1939 | 0 | call random_number( h1 ) ! In range 0-1 |
| 1940 | 0 | h1(:) = max(0., h1(:) - 0.05) ! Make 5% of values equal to zero |
| 1941 | 0 | call random_number( u0(1) ) ! In range 0-1 |
| 1942 | 0 | u0(:) = u0(1) ! Make u0 uniform |
| 1943 | ||
| 1944 | 0 | call remapping_core_h( remapCS, n0, h0, u0, n0, h1, u1 ) |
| 1945 | 0 | if ( maxval( abs( u1(:) - u0(1) ) ) > 0. ) then |
| 1946 | 0 | if ( .not. error ) then ! Only dump first error |
| 1947 | 0 | print *,'iter=',iter |
| 1948 | 0 | print *,'u0(1)',u0(1) |
| 1949 | 0 | print *,'u1',u1 |
| 1950 | 0 | print *,'u1-u0(1)',u1 - u0(1) |
| 1951 | endif | |
| 1952 | 0 | error = .true. |
| 1953 | endif | |
| 1954 | ||
| 1955 | enddo | |
| 1956 | ||
| 1957 | 0 | write(label,'(I0)') niter |
| 1958 | 0 | call test%test( error, trim(label)//' uniformity tests of '//scheme ) |
| 1959 | ||
| 1960 | 0 | end 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 | |
| 1968 | 0 | subroutine 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 | 0 | real :: h0(n0), h1(n0) ! Source and target grids [H but really nondim] |
| 1977 | 0 | 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 | 0 | om4_remap_via_sub_cells=.false.) |
| 1988 | ||
| 1989 | 0 | error = .false. |
| 1990 | 0 | do iter = 1, niter |
| 1991 | 0 | call random_number( h0 ) ! In range 0-1 |
| 1992 | 0 | h0(:) = max(0., h0(:) - 0.00) ! Note we do NOT test with vanished layers |
| 1993 | 0 | h1(:) = h0(:) ! Exact copy |
| 1994 | 0 | call random_number( u0 ) ! In range 0-1 |
| 1995 | ||
| 1996 | 0 | call remapping_core_h( remapCS, n0, h0, u0, n0, h1, u1 ) |
| 1997 | 0 | if ( maxval( abs( u1(:) - u0(:) ) ) > epsilon(h0(1)) * maxval( abs( u0 ) ) ) then |
| 1998 | 0 | if ( .not. error ) then ! Only dump first error |
| 1999 | 0 | print *,'iter=',iter |
| 2000 | 0 | print *,'h0',h0 |
| 2001 | 0 | print *,'u0',u0 |
| 2002 | 0 | print *,'u1',u1 |
| 2003 | 0 | print *,'u1-u0',u1 - u0 |
| 2004 | endif | |
| 2005 | 0 | error = .true. |
| 2006 | endif | |
| 2007 | ||
| 2008 | enddo | |
| 2009 | ||
| 2010 | 0 | write(label,'(I0)') niter |
| 2011 | 0 | call test%test( error, trim(label)//' unchanged grid tests of '//scheme ) |
| 2012 | ||
| 2013 | 0 | call remapCS%reconstruction%destroy() |
| 2014 | ||
| 2015 | 0 | end subroutine test_unchanged_grid |
| 2016 | ||
| 2017 | !> Test class-based remapping bitwise reproduces original implementation | |
| 2018 | 0 | subroutine 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 | 0 | real :: h0(n0), h1(n1) ! Source and target grids [H but really nondim] |
| 2030 | 0 | 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 | 0 | integer, allocatable :: seed(:) ! Random number seed |
| 2035 | character(len=16) :: label ! Generated label | |
| 2036 | ||
| 2037 | 0 | call random_seed(size=seed_size) |
| 2038 | 0 | allocate( seed(seed_Size) ) |
| 2039 | 0 | seed(:) = 102030405 ! Repeatable sequences |
| 2040 | 0 | call random_seed(put=seed) |
| 2041 | ||
| 2042 | 0 | error = .false. |
| 2043 | 0 | do iter = 1, niter |
| 2044 | 0 | call random_number( h0 ) ! In range 0-1 |
| 2045 | 0 | h0(:) = max(0., h0(:) - 0.00) ! Make 5% of values equal to zero |
| 2046 | 0 | h0(:) = h0(:) / sum( h0 ) ! Approximately normalize to total depth of 1 |
| 2047 | 0 | call random_number(h1) ! In range 0-1 |
| 2048 | 0 | h1(:) = max(0., h1(:) - 0.00) ! Make 5% of values equal to zero |
| 2049 | 0 | h1(:) = h1(:) / sum( h1 ) ! Approximately normalize to total depth of 1 |
| 2050 | 0 | call random_number( u0 ) ! In range 0-1 |
| 2051 | ||
| 2052 | 0 | call remapping_core_h( CS1, n0, h0, u0, n1, h1, u1 ) |
| 2053 | 0 | call remapping_core_h( CS2, n0, h0, u0, n1, h1, u2 ) |
| 2054 | 0 | error = sum( abs( u2(:) - u1(:) ) ) > 0. |
| 2055 | 0 | if (error) then |
| 2056 | 0 | print *,'iter=',iter |
| 2057 | 0 | print *,'h1',h1 |
| 2058 | 0 | print *,'h0',h0 |
| 2059 | 0 | print *,'u0',u0 |
| 2060 | 0 | print *,'u1',u1 |
| 2061 | 0 | print *,'u2',u2 |
| 2062 | 0 | 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 | 0 | exit |
| 2068 | endif | |
| 2069 | enddo | |
| 2070 | ||
| 2071 | 0 | write(label,'(I0)') niter |
| 2072 | 0 | call test%test( error, trim(label)//' comparisons of '//msg ) |
| 2073 | ||
| 2074 | 0 | end 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 | |
| 2079 | 0 | logical 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 | 0 | real, allocatable :: h0(:), h1(:), h2(:) ! Thicknesses for test columns [H] |
| 2086 | 0 | real, allocatable :: u0(:), u1(:), u2(:) ! Values for test profiles [A] |
| 2087 | 0 | real, allocatable :: dx1(:) ! Change in interface position [H] |
| 2088 | type(remapping_CS) :: CS, CS2 !< Remapping control structures | |
| 2089 | 0 | real, allocatable, dimension(:,:) :: ppoly0_E ! Edge values of polynomials [A] |
| 2090 | 0 | real, allocatable, dimension(:,:) :: ppoly0_S ! Edge slopes of polynomials [A H-1] |
| 2091 | 0 | real, allocatable, dimension(:,:) :: ppoly0_coefs ! Coefficients of polynomials [A] |
| 2092 | 0 | real, allocatable, dimension(:) :: h_sub, h0_eff ! Subgrid and effective source thicknesses [H] |
| 2093 | 0 | real, allocatable, dimension(:) :: u_sub, uh_sub ! Subgrid values and totals [A, A H] |
| 2094 | real :: u02_err ! Error in remaping [A] | |
| 2095 | 0 | 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 | 0 | integer, allocatable :: seed(:) ! Random number seed |
| 2101 | 0 | 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 | 0 | type(PCM) :: PCM_instance |
| 2106 | 0 | type(PLM_CW) :: PLM_CW_instance |
| 2107 | 0 | type(PLM_hybgen) :: PLM_hybgen_instance |
| 2108 | 0 | type(MPLM_WA) :: MPLM_WA_instance |
| 2109 | 0 | type(EMPLM_WA) :: EMPLM_WA_instance |
| 2110 | 0 | type(MPLM_WA_poly) :: MPLM_WA_poly_instance |
| 2111 | 0 | type(EMPLM_WA_poly) :: EMPLM_WA_poly_instance |
| 2112 | 0 | type(PLM_CWK) :: PLM_CWK_instance |
| 2113 | 0 | type(MPLM_CWK) :: MPLM_CWK_instance |
| 2114 | 0 | type(EMPLM_CWK) :: EMPLM_CWK_instance |
| 2115 | 0 | type(PPM_H4_2019) :: PPM_H4_2019_instance |
| 2116 | 0 | type(PPM_H4_2018) :: PPM_H4_2018_instance |
| 2117 | 0 | type(PPM_CW) :: PPM_CW_instance |
| 2118 | 0 | type(PPM_hybgen) :: PPM_hybgen_instance |
| 2119 | 0 | type(PPM_CWK) :: PPM_CWK_instance |
| 2120 | 0 | type(EPPM_CWK) :: EPPM_CWK_instance |
| 2121 | 0 | type(PLM_WLS) :: PLM_WLS_instance |
| 2122 | ||
| 2123 | 0 | call test%set( verbose=verbose ) ! Sets the verbosity flag in test |
| 2124 | ! call test%set( stop_instantly=.true. ) ! While debugging | |
| 2125 | ||
| 2126 | 0 | answer_date = 20190101 ! 20181231 |
| 2127 | 0 | h_neglect = 1.0e-30 |
| 2128 | 0 | h_neglect_edge = h_neglect ; if (answer_date < 20190101) h_neglect_edge = 1.0e-10 |
| 2129 | ||
| 2130 | 0 | if (verbose) write(test%stdout,*) ' ===== MOM_remapping: remapping_unit_tests =================' |
| 2131 | ||
| 2132 | 0 | if (verbose) write(test%stdout,*) ' - - - - - 1st generation tests - - - - -' |
| 2133 | ||
| 2134 | call initialize_remapping(CS, 'PPM_H4', answer_date=answer_date, & | |
| 2135 | 0 | 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 | 0 | n0 = 4 |
| 2139 | 0 | allocate( h0(n0), u0(n0) ) |
| 2140 | 0 | h0 = (/0.75, 0.75, 0.75, 0.75/) |
| 2141 | 0 | u0 = (/9., 3., -3., -9./) |
| 2142 | ||
| 2143 | ! Profile 1: 3 layers of thickness 1.0 and total depth 3 | |
| 2144 | 0 | n1 = 3 |
| 2145 | 0 | allocate( h1(n1), u1(n1), dx1(n1+1) ) |
| 2146 | 0 | h1 = (/1.0, 1.0, 1.0/) |
| 2147 | ||
| 2148 | ! Profile 2: 6 layers of thickness 0.5 and total depth 3 | |
| 2149 | 0 | n2 = 6 |
| 2150 | 0 | allocate( h2(n2), u2(n2) ) |
| 2151 | 0 | h2 = (/0.5, 0.5, 0.5, 0.5, 0.5, 0.5/) |
| 2152 | ||
| 2153 | ! Mapping u1 from h1 to h2 | |
| 2154 | 0 | call dzFromH1H2( n0, h0, n1, h1, dx1 ) |
| 2155 | 0 | call remapping_core_w( CS, n0, h0, u0, n1, dx1, u1 ) |
| 2156 | 0 | call test%real_arr(3, u1, (/8.,0.,-8./), 'remapping_core_w() PPM_H4') |
| 2157 | ||
| 2158 | 0 | allocate(ppoly0_E(n0,2), ppoly0_S(n0,2), ppoly0_coefs(n0,CS%degree+1)) |
| 2159 | 0 | ppoly0_E(:,:) = 0.0 |
| 2160 | 0 | ppoly0_S(:,:) = 0.0 |
| 2161 | 0 | ppoly0_coefs(:,:) = 0.0 |
| 2162 | ||
| 2163 | 0 | call initialize_remapping(CS, 'PPM_H4', force_bounds_in_subcell=.false., answer_date=answer_date) |
| 2164 | ||
| 2165 | 0 | call remapping_core_h( CS, n0, h0, u0, n2, h2, u2, net_err=err ) |
| 2166 | 0 | call test%real_arr(6, u2, (/10.,6.,2.,-2.,-6.,-10./), 'remapping_core_h() 2') |
| 2167 | ||
| 2168 | 0 | call remapping_core_h( CS, n0, h0, u0, 6, (/.125,.125,.125,.125,.125,.125/), u2, net_err=err ) |
| 2169 | 0 | call test%real_arr(6, u2, (/11.5,10.5,9.5,8.5,7.5,6.5/), 'remapping_core_h() 3') |
| 2170 | ||
| 2171 | 0 | call remapping_core_h( CS, n0, h0, u0, 3, (/2.25,1.5,1./), u2, net_err=err ) |
| 2172 | 0 | call test%real_arr(3, u2, (/3.,-10.5,-12./), 'remapping_core_h() 4') |
| 2173 | ||
| 2174 | 0 | deallocate(h0, u0, h1, u1, h2, u2, ppoly0_E, ppoly0_S, ppoly0_coefs) |
| 2175 | 0 | call end_remapping(CS) |
| 2176 | ||
| 2177 | ! =============================================== | |
| 2178 | ! This section tests the reconstruction functions | |
| 2179 | ! =============================================== | |
| 2180 | 0 | if (verbose) write(test%stdout,*) ' - - - - - reconstruction tests - - - - -' |
| 2181 | ||
| 2182 | 0 | 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 | 0 | ppoly0_E(1:3,:), ppoly0_coefs(1:3,:) ) |
| 2186 | 0 | call test%real_arr(3, ppoly0_E(:,1), (/1.,2.,4./), 'PCM: left edges') |
| 2187 | 0 | call test%real_arr(3, ppoly0_E(:,2), (/1.,2.,4./), 'PCM: right edges') |
| 2188 | 0 | 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 | 0 | ppoly0_E(1:3,:), ppoly0_coefs(1:3,:), h_neglect ) |
| 2192 | 0 | call test%real_arr(3, ppoly0_E(:,1), (/1.,2.,5./), 'Unlim PLM: left edges') |
| 2193 | 0 | call test%real_arr(3, ppoly0_E(:,2), (/1.,4.,5./), 'Unlim PLM: right edges') |
| 2194 | 0 | call test%real_arr(3, ppoly0_coefs(:,1), (/1.,2.,5./), 'Unlim PLM: P0') |
| 2195 | 0 | 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 | 0 | ppoly0_E(1:3,:), ppoly0_coefs(1:3,:), h_neglect ) |
| 2199 | 0 | call test%real_arr(3, ppoly0_E(:,1), (/1.,1.,7./), 'Left lim PLM: left edges') |
| 2200 | 0 | call test%real_arr(3, ppoly0_E(:,2), (/1.,3.,7./), 'Left lim PLM: right edges') |
| 2201 | 0 | call test%real_arr(3, ppoly0_coefs(:,1), (/1.,1.,7./), 'Left lim PLM: P0') |
| 2202 | 0 | 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 | 0 | ppoly0_E(1:3,:), ppoly0_coefs(1:3,:), h_neglect ) |
| 2206 | 0 | call test%real_arr(3, ppoly0_E(:,1), (/1.,5.,7./), 'Right lim PLM: left edges') |
| 2207 | 0 | call test%real_arr(3, ppoly0_E(:,2), (/1.,7.,7./), 'Right lim PLM: right edges') |
| 2208 | 0 | call test%real_arr(3, ppoly0_coefs(:,1), (/1.,5.,7./), 'Right lim PLM: P0') |
| 2209 | 0 | 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 | 0 | ppoly0_E(1:3,:), ppoly0_coefs(1:3,:), h_neglect ) |
| 2213 | 0 | call test%real_arr(3, ppoly0_E(:,1), (/1.,2.,9./), 'Non-uniform line PLM: left edges') |
| 2214 | 0 | call test%real_arr(3, ppoly0_E(:,2), (/1.,6.,9./), 'Non-uniform line PLM: right edges') |
| 2215 | 0 | call test%real_arr(3, ppoly0_coefs(:,1), (/1.,2.,9./), 'Non-uniform line PLM: P0') |
| 2216 | 0 | 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 | 0 | 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 | 0 | call test%real_arr(5, ppoly0_E(:,1), (/0.,2.,4.,6.,8./), 'Line H4: left edges', tol=8.0e-15) |
| 2222 | 0 | call test%real_arr(5, ppoly0_E(:,2), (/2.,4.,6.,8.,10./), 'Line H4: right edges', tol=1.0e-14) |
| 2223 | ||
| 2224 | 0 | ppoly0_E(:,1) = (/0.,2.,4.,6.,8./) |
| 2225 | 0 | 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 | 0 | ppoly0_coefs(1:5,:), h_neglect, answer_date=answer_date ) |
| 2228 | 0 | call test%real_arr(5, ppoly0_coefs(:,1), (/1.,2.,4.,6.,9./), 'Line PPM: P0') |
| 2229 | 0 | call test%real_arr(5, ppoly0_coefs(:,2), (/0.,2.,2.,2.,0./), 'Line PPM: P1') |
| 2230 | 0 | 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 | 0 | 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 | 0 | call test%real_arr(5, ppoly0_E(:,1), (/3.,0.,3.,12.,27./), 'Parabola H4: left edges', tol=2.7e-14) |
| 2236 | 0 | call test%real_arr(5, ppoly0_E(:,2), (/0.,3.,12.,27.,48./), 'Parabola H4: right edges', tol=4.8e-14) |
| 2237 | 0 | ppoly0_E(:,1) = (/0.,0.,3.,12.,27./) |
| 2238 | 0 | 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 | 0 | ppoly0_coefs(1:5,:), h_neglect, answer_date=answer_date ) |
| 2241 | 0 | call test%real_arr(5, ppoly0_E(:,1), (/0.,0.,3.,12.,37./), 'Parabola PPM: left edges') |
| 2242 | 0 | call test%real_arr(5, ppoly0_E(:,2), (/0.,3.,12.,27.,37./), 'Parabola PPM: right edges') |
| 2243 | 0 | call test%real_arr(5, ppoly0_coefs(:,1), (/0.,0.,3.,12.,37./), 'Parabola PPM: P0') |
| 2244 | 0 | call test%real_arr(5, ppoly0_coefs(:,2), (/0.,0.,6.,12.,0./), 'Parabola PPM: P1') |
| 2245 | 0 | call test%real_arr(5, ppoly0_coefs(:,3), (/0.,3.,3.,3.,0./), 'Parabola PPM: P2') |
| 2246 | ||
| 2247 | 0 | ppoly0_E(:,1) = (/0.,0.,6.,10.,15./) |
| 2248 | 0 | 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 | 0 | ppoly0_coefs(1:5,:), h_neglect, answer_date=answer_date ) |
| 2251 | 0 | call test%real_arr(5, ppoly0_E(:,1), (/0.,3.,6.,16.,15./), 'Limits PPM: left edges') |
| 2252 | 0 | call test%real_arr(5, ppoly0_E(:,2), (/0.,6.,9.,16.,15./), 'Limits PPM: right edges') |
| 2253 | 0 | call test%real_arr(5, ppoly0_coefs(:,1), (/0.,3.,6.,16.,15./), 'Limits PPM: P0') |
| 2254 | 0 | call test%real_arr(5, ppoly0_coefs(:,2), (/0.,6.,0.,0.,0./), 'Limits PPM: P1') |
| 2255 | 0 | call test%real_arr(5, ppoly0_coefs(:,3), (/0.,-3.,3.,0.,0./), 'Limits PPM: P2') |
| 2256 | ||
| 2257 | 0 | deallocate(ppoly0_E, ppoly0_S, ppoly0_coefs, u2) |
| 2258 | ||
| 2259 | ! ============================================================== | |
| 2260 | ! This section tests the components of remapping_core_h() | |
| 2261 | ! ============================================================== | |
| 2262 | ||
| 2263 | 0 | 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 | 0 | 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 | 0 | isrc_start, isrc_end, isrc_max, itgt_start, itgt_end, isub_src ) |
| 2280 | 0 | if (verbose) write(test%stdout,*) "intersect_src_tgt_grids test 1: n0=2, n1=3" |
| 2281 | 0 | if (verbose) write(test%stdout,*) " h_src = | 3 | 3 |" |
| 2282 | 0 | if (verbose) write(test%stdout,*) " h_tgt = | 2 | 2 | 2 |" |
| 2283 | 0 | call test%real_arr(6, h_sub, (/0.,2.,1.,1.,2.,0./), 'h_sub') |
| 2284 | 0 | call test%real_arr(2, h0_eff, (/3.,3./), 'h0_eff') |
| 2285 | 0 | call test%int_arr(2, isrc_start, (/1,4/), 'isrc_start') |
| 2286 | 0 | call test%int_arr(2, isrc_end, (/3,5/), 'isrc_end') |
| 2287 | 0 | call test%int_arr(2, isrc_max, (/2,5/), 'isrc_max') |
| 2288 | 0 | call test%int_arr(3, itgt_start, (/1,3,5/), 'itgt_start') |
| 2289 | 0 | call test%int_arr(3, itgt_end, (/2,4,6/), 'itgt_end') |
| 2290 | 0 | call test%int_arr(6, isub_src, (/1,1,1,2,2,2/), 'isub_src') |
| 2291 | 0 | 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 | 0 | 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 | 0 | isrc_start, isrc_end, isrc_max, itgt_start, itgt_end, isub_src ) |
| 2308 | 0 | if (verbose) write(test%stdout,*) "intersect_src_tgt_grids test 2: n0=3, n1=2" |
| 2309 | 0 | if (verbose) write(test%stdout,*) " h_src = | 2 | 2 | 2 |" |
| 2310 | 0 | if (verbose) write(test%stdout,*) " h_tgt = | 3 | 3 |" |
| 2311 | 0 | call test%real_arr(6, h_sub, (/0.,2.,1.,1.,2.,0./), 'h_sub') |
| 2312 | 0 | call test%real_arr(3, h0_eff, (/2.,2.,2./), 'h0_eff') |
| 2313 | 0 | call test%int_arr(3, isrc_start, (/1,3,5/), 'isrc_start') |
| 2314 | 0 | call test%int_arr(3, isrc_end, (/2,4,5/), 'isrc_end') |
| 2315 | 0 | call test%int_arr(3, isrc_max, (/2,4,5/), 'isrc_max') |
| 2316 | 0 | call test%int_arr(2, itgt_start, (/1,4/), 'itgt_start') |
| 2317 | 0 | call test%int_arr(2, itgt_end, (/3,6/), 'itgt_end') |
| 2318 | 0 | call test%int_arr(6, isub_src, (/1,1,2,2,3,3/), 'isub_src') |
| 2319 | 0 | 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 | 0 | n0 = 2 ; n1 = 3 |
| 2323 | 0 | allocate( h0_eff(n0), isrc_start(n0), isrc_end(n0), isrc_max(n0), h0(n0), u0(n0) ) |
| 2324 | 0 | allocate( itgt_start(n1), itgt_end(n1), h1(n1), u1(n1) ) |
| 2325 | 0 | allocate( h_sub(n0+n1+1), isub_src(n0+n1+1) ) |
| 2326 | 0 | u0 = (/ 2. , 5. /) |
| 2327 | 0 | h0 = (/ 2. , 4. /) |
| 2328 | 0 | 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 | 0 | isrc_start, isrc_end, isrc_max, itgt_start, itgt_end, isub_src ) |
| 2340 | 0 | if (verbose) write(test%stdout,*) "intersect_src_tgt_grids test 3: n0=2, n1=3" |
| 2341 | 0 | if (verbose) write(test%stdout,*) " h_src = | 2 | 4 |" |
| 2342 | 0 | if (verbose) write(test%stdout,*) " h_tgt = | 2 | 2 | 2 |" |
| 2343 | 0 | call test%real_arr(6, h_sub, (/0.,2.,0.,2.,2.,0./), 'h_sub') |
| 2344 | 0 | call test%real_arr(2, h0_eff, (/2.,4./), 'h0_eff') |
| 2345 | 0 | call test%int_arr(2, isrc_start, (/1,3/), 'isrc_start') |
| 2346 | 0 | call test%int_arr(2, isrc_end, (/2,5/), 'isrc_end') |
| 2347 | 0 | call test%int_arr(2, isrc_max, (/2,5/), 'isrc_max') |
| 2348 | 0 | call test%int_arr(3, itgt_start, (/1,4,5/), 'itgt_start') |
| 2349 | 0 | call test%int_arr(3, itgt_end, (/3,4,6/), 'itgt_end') |
| 2350 | 0 | call test%int_arr(6, isub_src, (/1,1,2,2,2,2/), 'isub_src') |
| 2351 | 0 | 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 | 0 | call PLM_reconstruction(n0, h0, u0, ppoly0_E, ppoly0_coefs, h_neglect ) |
| 2358 | 0 | call PLM_boundary_extrapolation(n0, h0, u0, ppoly0_E, ppoly0_coefs, h_neglect) |
| 2359 | 0 | 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 | 0 | INTEGRATION_PLM, .false., u_sub, uh_sub, u02_err) |
| 2363 | 0 | 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 | 0 | INTEGRATION_PLM, .false., u_sub, uh_sub, u02_err) |
| 2367 | 0 | 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 | 0 | .false., .false., .false., u1, u02_err) |
| 2374 | 0 | 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 | 0 | .true., .false., .false., u1, u02_err) |
| 2377 | 0 | call test%real_arr(3, u1, (/2.,4.,6./), 'u1.b') |
| 2378 | 0 | deallocate( ppoly0_coefs, ppoly0_E, ppoly0_S, u_sub, uh_sub, h0, u0, h1, u1) |
| 2379 | 0 | 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 | 0 | n0 = 2 ; n1 = 3 |
| 2383 | 0 | allocate( h0_eff(n0), isrc_start(n0), isrc_end(n0), isrc_max(n0), h0(n0), u0(n0) ) |
| 2384 | 0 | allocate( itgt_start(n1), itgt_end(n1), h1(n1), u1(n1) ) |
| 2385 | 0 | allocate( h_sub(n0+n1+1), isub_src(n0+n1+1) ) |
| 2386 | 0 | u0 = (/ 2. , 5. /) |
| 2387 | 0 | h0 = (/ 2. , 4. /) |
| 2388 | 0 | 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 | 0 | isrc_start, isrc_end, isrc_max, itgt_start, itgt_end, isub_src ) |
| 2400 | 0 | if (verbose) write(test%stdout,*) "intersect_src_tgt_grids test 4: n0=2, n1=3" |
| 2401 | 0 | if (verbose) write(test%stdout,*) " h_src = | 2 | 4 |" |
| 2402 | 0 | if (verbose) write(test%stdout,*) " h_tgt = | 2 | 2 | 1 |" |
| 2403 | 0 | call test%real_arr(6, h_sub, (/0.,2.,0.,2.,1.,1./), 'h_sub') |
| 2404 | 0 | call test%real_arr(2, h0_eff, (/2.,3./), 'h0_eff') |
| 2405 | 0 | call test%int_arr(2, isrc_start, (/1,3/), 'isrc_start') |
| 2406 | 0 | call test%int_arr(2, isrc_end, (/2,6/), 'isrc_end') |
| 2407 | 0 | call test%int_arr(2, isrc_max, (/2,4/), 'isrc_max') |
| 2408 | 0 | call test%int_arr(3, itgt_start, (/1,4,5/), 'itgt_start') |
| 2409 | 0 | call test%int_arr(3, itgt_end, (/3,4,5/), 'itgt_end') |
| 2410 | 0 | call test%int_arr(6, isub_src, (/1,1,2,2,2,2/), 'isub_src') |
| 2411 | 0 | 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 | 0 | call PLM_reconstruction(n0, h0, u0, ppoly0_E, ppoly0_coefs, h_neglect ) |
| 2418 | 0 | call PLM_boundary_extrapolation(n0, h0, u0, ppoly0_E, ppoly0_coefs, h_neglect) |
| 2419 | 0 | 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 | 0 | INTEGRATION_PLM, .false., u_sub, uh_sub, u02_err) |
| 2423 | 0 | call test%real_arr(6, u_sub, (/1.,2.,3.,4.,5.5,6.5/), 'u_sub') |
| 2424 | 0 | deallocate( ppoly0_coefs, ppoly0_E, ppoly0_S, u_sub, uh_sub, h0, u0, h1, u1) |
| 2425 | 0 | 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 | 0 | n0 = 3 ; n1 = 2 |
| 2429 | 0 | allocate( h0_eff(n0), isrc_start(n0), isrc_end(n0), isrc_max(n0), h0(n0), u0(n0) ) |
| 2430 | 0 | allocate( itgt_start(n1), itgt_end(n1), h1(n1), u1(n1) ) |
| 2431 | 0 | allocate( h_sub(n0+n1+1), isub_src(n0+n1+1) ) |
| 2432 | 0 | u0 = (/ 2. , 4. , 5.5 /) |
| 2433 | 0 | h0 = (/ 2. , 2. , 1. /) |
| 2434 | 0 | 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 | 0 | isrc_start, isrc_end, isrc_max, itgt_start, itgt_end, isub_src ) |
| 2446 | 0 | if (verbose) write(test%stdout,*) "intersect_src_tgt_grids test 5: n0=3, n1=2" |
| 2447 | 0 | if (verbose) write(test%stdout,*) " h_src = | 2 | 2 | 1 |" |
| 2448 | 0 | if (verbose) write(test%stdout,*) " h_tgt = | 2 | 4 |" |
| 2449 | 0 | call test%real_arr(6, h_sub, (/0.,2.,0.,2.,1.,1./), 'h_sub') |
| 2450 | 0 | call test%real_arr(3, h0_eff, (/2.,2.,1./), 'h0_eff') |
| 2451 | 0 | call test%int_arr(3, isrc_start, (/1,3,5/), 'isrc_start') |
| 2452 | 0 | call test%int_arr(3, isrc_end, (/2,4,5/), 'isrc_end') |
| 2453 | 0 | call test%int_arr(3, isrc_max, (/2,4,5/), 'isrc_max') |
| 2454 | 0 | call test%int_arr(2, itgt_start, (/1,4/), 'itgt_start') |
| 2455 | 0 | call test%int_arr(2, itgt_end, (/3,6/), 'itgt_end') |
| 2456 | 0 | call test%int_arr(6, isub_src, (/1,1,2,2,3,3/), 'isub_src') |
| 2457 | 0 | 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 | 0 | call PLM_reconstruction(n0, h0, u0, ppoly0_E, ppoly0_coefs, h_neglect ) |
| 2464 | 0 | call PLM_boundary_extrapolation(n0, h0, u0, ppoly0_E, ppoly0_coefs, h_neglect) |
| 2465 | 0 | 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 | 0 | INTEGRATION_PLM, .false., u_sub, uh_sub, u02_err) |
| 2469 | 0 | 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 | 0 | INTEGRATION_PLM, .false., u_sub, uh_sub, u02_err) |
| 2473 | 0 | 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 | 0 | .false., .false., .false., u1, u02_err) |
| 2480 | 0 | 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 | 0 | .true., .false., .false., u1, u02_err) |
| 2483 | 0 | call test%real_arr(2, u1, (/2.,4.875/), 'u1.b') |
| 2484 | 0 | deallocate( ppoly0_coefs, ppoly0_E, ppoly0_S, u_sub, uh_sub, h0, u0, h1, u1) |
| 2485 | 0 | 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 | 0 | n0 = 3 ; n1 = 5 |
| 2489 | 0 | allocate( h0_eff(n0), isrc_start(n0), isrc_end(n0), isrc_max(n0), h0(n0), u0(n0) ) |
| 2490 | 0 | allocate( itgt_start(n1), itgt_end(n1), h1(n1), u1(n1) ) |
| 2491 | 0 | allocate( h_sub(n0+n1+1), isub_src(n0+n1+1) ) |
| 2492 | 0 | u0 = (/ 2. ,3., 4. /) |
| 2493 | 0 | h0 = (/ 2. ,0., 2. /) |
| 2494 | 0 | 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 | 0 | isrc_start, isrc_end, isrc_max, itgt_start, itgt_end, isub_src ) |
| 2506 | 0 | if (verbose) write(test%stdout,*) "intersect_src_tgt_grids test 6: n0=3, n1=5" |
| 2507 | 0 | if (verbose) write(test%stdout,*) " h_src = | 2 |0| 2 |" |
| 2508 | 0 | if (verbose) write(test%stdout,*) " h_tgt = | 1 |0| 1 |0| 2 |" |
| 2509 | 0 | call test%real_arr(9, h_sub, (/0.,1.,0.,1.,0.,0.,0.,2.,0./), 'h_sub') |
| 2510 | 0 | call test%real_arr(3, h0_eff, (/2.,0.,2./), 'h0_eff') |
| 2511 | 0 | call test%int_arr(3, isrc_start, (/1,5,6/), 'isrc_start') |
| 2512 | 0 | call test%int_arr(3, isrc_end, (/4,5,8/), 'isrc_end') |
| 2513 | 0 | call test%int_arr(3, isrc_max, (/4,5,8/), 'isrc_max') |
| 2514 | 0 | call test%int_arr(5, itgt_start, (/1,3,4,7,8/), 'itgt_start') |
| 2515 | 0 | call test%int_arr(5, itgt_end, (/2,3,6,7,9/), 'itgt_end') |
| 2516 | 0 | call test%int_arr(9, isub_src, (/1,1,1,1,2,3,3,3,3/), 'isub_src') |
| 2517 | 0 | 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 | 0 | call PLM_reconstruction(n0, h0, u0, ppoly0_E, ppoly0_coefs, h_neglect ) |
| 2524 | 0 | call PLM_boundary_extrapolation(n0, h0, u0, ppoly0_E, ppoly0_coefs, h_neglect) |
| 2525 | 0 | 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 | 0 | INTEGRATION_PLM, .false., u_sub, uh_sub, u02_err) |
| 2529 | 0 | 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 | 0 | INTEGRATION_PLM, .false., u_sub, uh_sub, u02_err) |
| 2533 | 0 | 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 | 0 | .false., .false., .false., u1, u02_err) |
| 2540 | 0 | 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 | 0 | .true., .false., .false., u1, u02_err) |
| 2543 | 0 | call test%real_arr(5, u1, (/1.5,2.,2.5,3.,4./), 'u1.b') |
| 2544 | 0 | deallocate( ppoly0_coefs, ppoly0_E, ppoly0_S, u_sub, uh_sub, h0, u0, h1, u1) |
| 2545 | 0 | 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 | 0 | if (verbose) write(test%stdout,*) '- - - - - - - - - - remapping_core_h() tests - - - - - - - - -' |
| 2551 | ||
| 2552 | 0 | allocate(u2(2)) |
| 2553 | ||
| 2554 | 0 | 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 | 0 | call remapping_core_h(CS, 4, (/0.,1.,1.,0./), (/5.,4.,2.,1./), 2, (/1.,1./), u2) |
| 2558 | 0 | call test%real_arr(2, u2, (/4.,2./), 'PLM: remapped h=0110->h=11 om4') |
| 2559 | 0 | call remapping_core_h(CS, 4, (/0.,1.,1.,0./), (/5.,4.,2.,1./), 2, (/1.,1./), u2) |
| 2560 | 0 | 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 | 0 | call remapping_core_h(CS, 4, (/0.,1.,1.,0./), (/5.,4.,2.,1./), 2, (/1.,4./), u2) |
| 2566 | 0 | call test%real_arr(2, u2, (/4.,1.25/), 'PLM: remapped h=0110->h=14 om4') |
| 2567 | 0 | call remapping_core_h(CS, 4, (/0.,1.,1.,0./), (/5.,4.,2.,1./), 2, (/1.,4./), u2) |
| 2568 | 0 | 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 | 0 | call remapping_set_param(CS, om4_remap_via_sub_cells=.true.) |
| 2574 | 0 | call remapping_core_h(CS, 4, (/0.,4.,4.,0./), (/5.,4.,2.,1./), 2, (/4.,2./), u2) |
| 2575 | 0 | call test%real_arr(2, u2, (/4.,2./), 'PLM: remapped h=0440->h=42 om4 (with known bug)') |
| 2576 | 0 | call remapping_set_param(CS, om4_remap_via_sub_cells=.false.) |
| 2577 | 0 | call remapping_core_h(CS, 4, (/0.,4.,4.,0./), (/5.,4.,2.,1./), 2, (/4.,2./), u2) |
| 2578 | 0 | 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 | 0 | call remapping_set_param(CS, om4_remap_via_sub_cells=.true.) |
| 2584 | 0 | call remapping_core_h(CS, 4, (/0.,5.,5.,0./), (/5.,4.,2.,1./), 2, (/2.,2./), u2) |
| 2585 | 0 | call test%real_arr(2, u2, (/4.5,3.5/), 'PLM: remapped h=0880->h=21 om4 (with known bug)') |
| 2586 | 0 | call remapping_set_param(CS, om4_remap_via_sub_cells=.false.) |
| 2587 | 0 | call remapping_core_h(CS, 4, (/0.,4.,4.,0./), (/5.,4.,2.,1./), 2, (/2.,1./), u2) |
| 2588 | 0 | call test%real_arr(2, u2, (/4.5,3.75/), 'PLM: remapped h=0440->h=21') |
| 2589 | ||
| 2590 | 0 | 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 | 0 | n0 = 8 |
| 2594 | 0 | allocate( h0(n0), u0(n0) ) |
| 2595 | 0 | h0 = (/0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0/) |
| 2596 | 0 | u0 = (/1.0, 1.5, 2.5, 3.5, 4.5, 5.5, 6.0, 6.0/) |
| 2597 | 0 | allocate( u1(8) ) |
| 2598 | ||
| 2599 | 0 | call initialize_remapping(CS, 'PLM', answer_date=99990101, h_neglect=1.e-17, h_neglect_edge=1.e-2) |
| 2600 | ||
| 2601 | 0 | do om4 = 0, 1 |
| 2602 | 0 | if ( om4 == 0 ) then |
| 2603 | 0 | CS%om4_remap_via_sub_cells = .false. |
| 2604 | 0 | om4_tag(:) = ' ' |
| 2605 | else | |
| 2606 | 0 | CS%om4_remap_via_sub_cells = .true. |
| 2607 | 0 | om4_tag(:) = ' om4' |
| 2608 | endif | |
| 2609 | ||
| 2610 | ! Unchanged grid | |
| 2611 | 0 | call remapping_core_h( CS, n0, h0, u0, 8, [0.,1.,1.,1.,1.,1.,0.,0.], u1) |
| 2612 | 0 | 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 | 0 | call remapping_core_h( CS, n0, h0, u0, 5, [1.,1.,1.,1.,1.], u1) |
| 2616 | 0 | 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 | 0 | call remapping_core_h( CS, n0, h0, u0, 4, [0.5,1.,1.5,2.], u1) |
| 2620 | 0 | 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 | 0 | call remapping_core_h( CS, n0, h0, u0, 6, [0.5,1.,0.,1.5,2.,0.], u1) |
| 2624 | 0 | 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 | 0 | call remapping_core_h( CS, n0, h0, u0, 5, [1.5,1.5,0.,4.,2.], u1) |
| 2628 | 0 | 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 | 0 | call remapping_core_h( CS, n0, h0, u0, 5, [1.,1.,1.,1.,0.5], u1) |
| 2632 | 0 | if ( om4 == 0 ) then |
| 2633 | 0 | 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 | 0 | 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 | 0 | call remapping_core_h( CS, n0, h0, u0, 3, [0.5,0.,1.], u1) |
| 2640 | 0 | if ( om4 == 0 ) then |
| 2641 | 0 | call test%real_arr(3, u1, (/1.25,1.5,2./), 'PLM: remapped h=01111100->h=h01') |
| 2642 | else | |
| 2643 | 0 | 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 | 0 | call end_remapping(CS) |
| 2649 | 0 | deallocate( h0, u0, u1 ) |
| 2650 | ||
| 2651 | ! ============================================================ | |
| 2652 | ! This section tests interpolation and reintegration functions | |
| 2653 | ! ============================================================ | |
| 2654 | 0 | 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 | 0 | 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 | 0 | 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 | 0 | 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 | 0 | 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 | 0 | 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 | 0 | 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 | 0 | 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 | 0 | 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 | 0 | 4, (/0.,2.,6.,0./), (/0.,1.,3.,8.,0./) ) |
| 2691 | ||
| 2692 | 0 | 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 | 0 | 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 | 0 | 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 | 0 | 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 | 0 | 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 | 0 | 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 | 0 | 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 | 0 | 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 | 0 | 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 | 0 | 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 | 0 | 3, (/0.,0.,0./), (/0.,0.,0./) ) |
| 2733 | ||
| 2734 | 0 | if (verbose) write(test%stdout,*) '- - - - - - - - - - Recon1d PCM tests - - - - - - - - -' |
| 2735 | 0 | call test%test( PCM_instance%unit_tests(verbose, test%stdout, test%stderr), 'PCM unit test') |
| 2736 | 0 | call test%test( MPLM_WA_instance%unit_tests(verbose, test%stdout, test%stderr), 'MPLM_WA unit test') |
| 2737 | 0 | call test%test( EMPLM_WA_instance%unit_tests(verbose, test%stdout, test%stderr), 'EMPLM_WA unit test') |
| 2738 | 0 | call test%test( MPLM_WA_poly_instance%unit_tests(verbose, test%stdout, test%stderr), 'MPLM_WA_poly unit test') |
| 2739 | 0 | call test%test( EMPLM_WA_poly_instance%unit_tests(verbose, test%stdout, test%stderr), 'EMPLM_WA_poly unit test') |
| 2740 | 0 | call test%test( PLM_hybgen_instance%unit_tests(verbose, test%stdout, test%stderr), 'PLM_hybgen unit test') |
| 2741 | 0 | call test%test( PLM_CW_instance%unit_tests(verbose, test%stdout, test%stderr), 'PLM_CW unit test') |
| 2742 | 0 | call test%test( PLM_CWK_instance%unit_tests(verbose, test%stdout, test%stderr), 'PLM_CWK unit test') |
| 2743 | 0 | call test%test( MPLM_CWK_instance%unit_tests(verbose, test%stdout, test%stderr), 'MPLM_CWK unit test') |
| 2744 | 0 | call test%test( EMPLM_CWK_instance%unit_tests(verbose, test%stdout, test%stderr), 'EMPLM_CWK unit test') |
| 2745 | 0 | call test%test( PPM_H4_2019_instance%unit_tests(verbose, test%stdout, test%stderr), 'PPM_H4_2019 unit test') |
| 2746 | 0 | call test%test( PPM_H4_2018_instance%unit_tests(verbose, test%stdout, test%stderr), 'PPM_H4_2018 unit test') |
| 2747 | 0 | call test%test( PPM_hybgen_instance%unit_tests(verbose, test%stdout, test%stderr), 'PPM_hybgen unit test') |
| 2748 | 0 | call test%test( PPM_CW_instance%unit_tests(verbose, test%stdout, test%stderr), 'PPM_CW unit test') |
| 2749 | 0 | call test%test( PPM_CWK_instance%unit_tests(verbose, test%stdout, test%stderr), 'PPM_CWK unit test') |
| 2750 | 0 | call test%test( EPPM_CWK_instance%unit_tests(verbose, test%stdout, test%stderr), 'EPPM_CWK unit test') |
| 2751 | 0 | call test%test( PLM_WLS_instance%unit_tests(verbose, test%stdout, test%stderr), 'PLM_WLS unit test') |
| 2752 | ||
| 2753 | ! Randomized, brute force tests | |
| 2754 | 0 | ntests = 3000 |
| 2755 | 0 | if (present(num_comp_samp)) ntests = num_comp_samp |
| 2756 | ||
| 2757 | 0 | call random_seed(size=seed_size) |
| 2758 | 0 | allocate( seed(seed_Size) ) |
| 2759 | 0 | seed(:) = 102030405 |
| 2760 | 0 | call random_seed(put=seed) |
| 2761 | ||
| 2762 | 0 | n0 = 9 |
| 2763 | ||
| 2764 | ! Internal consistency | |
| 2765 | 0 | call test_recon_consistency(test, 'C_PCM', n0, ntests, h_neglect) |
| 2766 | 0 | call test_recon_consistency(test, 'C_PLM_CW', n0, ntests, h_neglect) |
| 2767 | 0 | call test_recon_consistency(test, 'C_PLM_HYBGEN', n0, ntests, h_neglect) |
| 2768 | 0 | call test_recon_consistency(test, 'C_MPLM_WA', n0, ntests, h_neglect) |
| 2769 | 0 | call test_recon_consistency(test, 'C_EMPLM_WA', n0, ntests, h_neglect) |
| 2770 | 0 | call test_recon_consistency(test, 'C_MPLM_WA_poly', n0, ntests, h_neglect) |
| 2771 | 0 | call test_recon_consistency(test, 'C_EMPLM_WA_poly', n0, ntests, h_neglect) |
| 2772 | 0 | call test_recon_consistency(test, 'C_PLM_CWK', n0, ntests, h_neglect) |
| 2773 | 0 | call test_recon_consistency(test, 'C_MPLM_CWK', n0, ntests, h_neglect) |
| 2774 | 0 | call test_recon_consistency(test, 'C_EMPLM_CWK', n0, ntests, h_neglect) |
| 2775 | 0 | call test_recon_consistency(test, 'C_PPM_H4_2018', n0, ntests, h_neglect) |
| 2776 | 0 | call test_recon_consistency(test, 'C_PPM_H4_2019', n0, ntests, h_neglect) |
| 2777 | 0 | call test_recon_consistency(test, 'C_PPM_HYBGEN', n0, ntests, h_neglect) |
| 2778 | 0 | call test_recon_consistency(test, 'C_PPM_CW', n0, ntests, h_neglect) |
| 2779 | 0 | call test_recon_consistency(test, 'C_PPM_CWK', n0, ntests, h_neglect) |
| 2780 | 0 | call test_recon_consistency(test, 'C_EPPM_CWK', n0, ntests, h_neglect) |
| 2781 | 0 | call test_recon_consistency(test, 'C_PLM_WLS', n0, ntests, h_neglect) |
| 2782 | ||
| 2783 | 0 | call test_preserve_uniform(test, 'PCM', n0, ntests, h_neglect) |
| 2784 | 0 | 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 | 0 | call test_preserve_uniform(test, 'C_PLM_CW', n0, ntests, h_neglect) |
| 2794 | 0 | call test_preserve_uniform(test, 'C_PLM_HYBGEN', n0, ntests, h_neglect) |
| 2795 | 0 | call test_preserve_uniform(test, 'C_MPLM_WA', n0, ntests, h_neglect) |
| 2796 | 0 | call test_preserve_uniform(test, 'C_EMPLM_WA', n0, ntests, h_neglect) |
| 2797 | 0 | 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 | 0 | call test_preserve_uniform(test, 'C_PLM_CWK', n0, ntests, h_neglect) |
| 2800 | 0 | call test_preserve_uniform(test, 'C_MPLM_CWK', n0, ntests, h_neglect) |
| 2801 | 0 | call test_preserve_uniform(test, 'C_EMPLM_CWK', n0, ntests, h_neglect) |
| 2802 | 0 | call test_preserve_uniform(test, 'C_PPM_H4_2019', n0, ntests, h_neglect) |
| 2803 | 0 | call test_preserve_uniform(test, 'C_PPM_H4_2018', n0, ntests, h_neglect) |
| 2804 | 0 | call test_preserve_uniform(test, 'C_PPM_HYBGEN', n0, ntests, h_neglect) |
| 2805 | 0 | call test_preserve_uniform(test, 'C_PPM_CW', n0, ntests, h_neglect) |
| 2806 | 0 | call test_preserve_uniform(test, 'C_PPM_CWK', n0, ntests, h_neglect) |
| 2807 | 0 | call test_preserve_uniform(test, 'C_EPPM_CWK', n0, ntests, h_neglect) |
| 2808 | 0 | call test_preserve_uniform(test, 'C_PLM_WLS', n0, ntests, h_neglect) |
| 2809 | ||
| 2810 | 0 | call test_unchanged_grid(test, 'C_PCM', n0, ntests, h_neglect) |
| 2811 | 0 | call test_unchanged_grid(test, 'C_PLM_CW', n0, ntests, h_neglect) |
| 2812 | 0 | call test_unchanged_grid(test, 'C_PLM_HYBGEN', n0, ntests, h_neglect) |
| 2813 | 0 | call test_unchanged_grid(test, 'C_PLM_CWK', n0, ntests, h_neglect) |
| 2814 | 0 | call test_unchanged_grid(test, 'C_MPLM_CWK', n0, ntests, h_neglect) |
| 2815 | 0 | call test_unchanged_grid(test, 'C_EMPLM_CWK', n0, ntests, h_neglect) |
| 2816 | 0 | call test_unchanged_grid(test, 'C_PPM_HYBGEN', n0, ntests, h_neglect) |
| 2817 | 0 | call test_unchanged_grid(test, 'C_PPM_CW', n0, ntests, h_neglect) |
| 2818 | 0 | call test_unchanged_grid(test, 'C_PPM_CWK', n0, ntests, h_neglect) |
| 2819 | 0 | call test_unchanged_grid(test, 'C_EPPM_CWK', n0, ntests, h_neglect) |
| 2820 | 0 | 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 | 0 | allocate( h0(8), u0(8) ) |
| 2824 | 0 | h0 = (/0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0/) |
| 2825 | 0 | u0 = (/1.0, 1.5, 2.5, 3.5, 4.5, 5.5, 6.0, 6.0/) |
| 2826 | 0 | allocate( u1(8) ) |
| 2827 | 0 | call initialize_remapping(CS, 'C_PLM_CW', nk=8) |
| 2828 | 0 | call remapping_core_h( CS, 8, h0, u0, 8, [0.,1.,1.,1.,1.,1.,0.,0.], u1 ) |
| 2829 | 0 | call test%real_arr(8, u1, u0, 'remapping_core to unchanged grid with class') |
| 2830 | ||
| 2831 | 0 | call end_remapping(CS) |
| 2832 | 0 | deallocate( h0, u0, u1 ) |
| 2833 | ||
| 2834 | ! Brute force test that we have bitwise identical answers with the new classes | |
| 2835 | 0 | n0 = 7 |
| 2836 | 0 | 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 | 0 | 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 | 0 | h_neglect=h_neglect, h_neglect_edge=h_neglect_edge) |
| 2845 | 0 | 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 | 0 | 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 | 0 | h_neglect=h_neglect, h_neglect_edge=h_neglect_edge) |
| 2854 | 0 | 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 | 0 | 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 | 0 | h_neglect=h_neglect, h_neglect_edge=h_neglect_edge) |
| 2863 | 0 | 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 | 0 | 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 | 0 | h_neglect=h_neglect, h_neglect_edge=h_neglect_edge) |
| 2872 | 0 | 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 | 0 | force_bounds_in_subcell=.false.) |
| 2877 | 0 | call initialize_remapping(CS2, 'C_PCM', nk=n0, h_neglect=h_neglect) |
| 2878 | 0 | 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 | 0 | h_neglect=h_neglect, h_neglect_edge=h_neglect_edge) |
| 2884 | 0 | call initialize_remapping(CS2, 'C_MPLM_WA_POLY', nk=n0, h_neglect=h_neglect) |
| 2885 | 0 | 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 | 0 | h_neglect=h_neglect, h_neglect_edge=h_neglect_edge) |
| 2891 | 0 | call initialize_remapping(CS2, 'C_MPLM_WA_POLY', nk=n0, h_neglect=h_neglect) |
| 2892 | 0 | 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 | 0 | h_neglect=h_neglect, h_neglect_edge=h_neglect_edge) |
| 2898 | 0 | call initialize_remapping(CS2, 'C_EMPLM_WA_POLY', nk=n0, h_neglect=h_neglect) |
| 2899 | 0 | 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 | 0 | h_neglect=h_neglect, h_neglect_edge=h_neglect_edge) |
| 2905 | 0 | call initialize_remapping(CS2, 'C_EMPLM_WA_POLY', nk=n0, h_neglect=h_neglect) |
| 2906 | 0 | 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 | 0 | h_neglect=h_neglect, h_neglect_edge=h_neglect_edge) |
| 2912 | 0 | call initialize_remapping(CS2, 'C_PPM_H4_2018', nk=n0, h_neglect=h_neglect) |
| 2913 | 0 | 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 | 0 | h_neglect=h_neglect, h_neglect_edge=h_neglect_edge) |
| 2919 | 0 | call initialize_remapping(CS2, 'C_PPM_H4_2018', nk=n0, h_neglect=h_neglect) |
| 2920 | 0 | 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 | 0 | h_neglect=h_neglect, h_neglect_edge=h_neglect_edge) |
| 2926 | 0 | call initialize_remapping(CS2, 'C_PPM_H4_2019', nk=n0, h_neglect=h_neglect) |
| 2927 | 0 | 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 | 0 | h_neglect=h_neglect, h_neglect_edge=h_neglect_edge) |
| 2933 | 0 | call initialize_remapping(CS2, 'C_PPM_H4_2019', nk=n0, h_neglect=h_neglect) |
| 2934 | 0 | 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 | 0 | h_neglect=h_neglect, h_neglect_edge=h_neglect_edge) |
| 2940 | 0 | call initialize_remapping(CS2, 'C_PLM_hybgen', nk=n0, h_neglect=h_neglect) |
| 2941 | 0 | 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 | 0 | h_neglect=h_neglect, h_neglect_edge=h_neglect_edge) |
| 2947 | 0 | call initialize_remapping(CS2, 'C_PPM_HYBGEN', nk=n0, h_neglect=h_neglect) |
| 2948 | 0 | call compare_two_schemes(test, CS, CS2, n0, n1, ntests, 'PPM_HYBGEN bounded <-> C_PPM_HYGEN') |
| 2949 | ||
| 2950 | 0 | call end_remapping(CS) |
| 2951 | 0 | call end_remapping(CS2) |
| 2952 | ||
| 2953 | 0 | remapping_unit_tests = test%summarize('remapping_unit_tests') |
| 2954 | ||
| 2955 | 0 | end function remapping_unit_tests |
| 2956 | ||
| 2957 | 0 | end module MOM_remapping |