From 3c0b0a139eadb5a03dc97053f740a4e705e5c0eb Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Tue, 12 Mar 2024 14:52:16 -0400 Subject: [PATCH] Refactor/split remap_via_sub_cells() adding intersect_src_tgt_grids() - The first block of remap_via_sub_cells() computes the intersection between two columns (grids). I've split out this block into a stand alone function so that we can re-use it in a to-be-written analog of remap_via_sub_cells() that will remap integrated quantities. - Added some in-line documentation of algorithm used by remap_via_sub_cells() - Added new column-wise function intersect_src_tgt_grids() - Added tests of intersect_src_tgt_grids() that check against known results - Had to add a new helper function to MOM_remapping to report state of integer arrays. --- src/ALE/MOM_remapping.F90 | 588 +++++++++++++++++++++++++++++--------- 1 file changed, 452 insertions(+), 136 deletions(-) diff --git a/src/ALE/MOM_remapping.F90 b/src/ALE/MOM_remapping.F90 index 95c0b2b081..5b4b936874 100644 --- a/src/ALE/MOM_remapping.F90 +++ b/src/ALE/MOM_remapping.F90 @@ -477,8 +477,6 @@ subroutine remap_via_sub_cells( n0, h0, u0, ppoly0_E, ppoly0_coefs, n1, h1, meth integer :: i_sub ! Index of sub-cell integer :: i0 ! Index into h0(1:n0), source column integer :: i1 ! Index into h1(1:n1), target column - integer :: i_start0 ! Used to record which sub-cells map to source cells - integer :: i_start1 ! Used to record which sub-cells map to target cells integer :: i_max ! Used to record which sub-cell is the largest contribution of a source cell real :: dh_max ! Used to record which sub-cell is the largest contribution of a source cell [H] real, dimension(n0+n1+1) :: h_sub ! Width of each each sub-cell [H] @@ -494,7 +492,6 @@ subroutine remap_via_sub_cells( n0, h0, u0, ppoly0_E, ppoly0_coefs, n1, h1, meth integer, dimension(n1) :: itgt_start ! Index of first sub-cell within each target cell integer, dimension(n1) :: itgt_end ! Index of last sub-cell within each target cell real :: xa, xb ! Non-dimensional position within a source cell (0..1) [nondim] - real :: h0_supply, h1_supply ! The amount of width available for constructing sub-cells [H] real :: dh ! The width of the sub-cell [H] real :: duh ! The total amount of accumulated stuff (u*h) [A H] real :: dh0_eff ! Running sum of source cell thickness [H] @@ -509,8 +506,6 @@ subroutine remap_via_sub_cells( n0, h0, u0, ppoly0_E, ppoly0_coefs, n1, h1, meth real :: u0tot, u1tot, u2tot ! Integrated reconstruction values [H A] real :: u_orig ! The original value of the reconstruction in a cell [A] real :: u0min, u0max, u1min, u1max, u2min, u2max ! Minimum and maximum values of reconstructions [A] - logical :: src_has_volume !< True if h0 has not been consumed - logical :: tgt_has_volume !< True if h1 has not been consumed i0_last_thick_cell = 0 do i0 = 1, n0 @@ -519,134 +514,13 @@ subroutine remap_via_sub_cells( n0, h0, u0, ppoly0_E, ppoly0_coefs, n1, h1, meth if (h0(i0)>0.) i0_last_thick_cell = i0 enddo - ! Initialize algorithm - h0_supply = h0(1) - h1_supply = h1(1) - src_has_volume = .true. - tgt_has_volume = .true. - i0 = 1 ; i1 = 1 - i_start0 = 1 ; i_start1 = 1 - i_max = 1 - dh_max = 0. - dh0_eff = 0. - - ! First sub-cell is always vanished - h_sub(1) = 0. - isrc_start(1) = 1 - isrc_end(1) = 1 - isrc_max(1) = 1 - isub_src(1) = 1 - - ! Loop over each sub-cell to calculate intersections with source and target grids - do i_sub = 2, n0+n1+1 - - ! This is the width of the sub-cell, determined by which ever column has the least - ! supply available to consume. - dh = min(h0_supply, h1_supply) - - ! This is the running sum of the source cell thickness. After summing over each - ! sub-cell, the sum of sub-cell thickness might differ from the original source - ! cell thickness due to round off. - dh0_eff = dh0_eff + min(dh, h0_supply) - - ! Record the source index (i0) that this sub-cell integral belongs to. This - ! is needed to index the reconstruction coefficients for the source cell - ! used in the integrals of the sub-cell width. - isub_src(i_sub) = i0 - h_sub(i_sub) = dh - - ! For recording the largest sub-cell within a source cell. - if (dh >= dh_max) then - i_max = i_sub - dh_max = dh - endif - - ! Which ever column (source or target) has the least width left to consume determined - ! the width, dh, of sub-cell i_sub in the expression for dh above. - if (h0_supply <= h1_supply .and. src_has_volume) then - ! h0_supply is smaller than h1_supply) so we consume h0_supply and increment the - ! source cell index. - h1_supply = h1_supply - dh ! Although this is a difference the result will - ! be non-negative because of the conditional. - ! Record the sub-cell start/end index that span the source cell i0. - isrc_start(i0) = i_start0 - isrc_end(i0) = i_sub - i_start0 = i_sub + 1 - ! Record the sub-cell that is the largest fraction of the source cell. - isrc_max(i0) = i_max - i_max = i_sub + 1 - dh_max = 0. - ! Record the source cell thickness found by summing the sub-cell thicknesses. - h0_eff(i0) = dh0_eff - ! Move the source index. - if (i0 < n0) then - i0 = i0 + 1 - h0_supply = h0(i0) - dh0_eff = 0. - else - h0_supply = 0. - src_has_volume = .false. - endif - elseif (h0_supply >= h1_supply .and. tgt_has_volume) then - ! h1_supply is smaller than h0_supply) so we consume h1_supply and increment the - ! target cell index. - h0_supply = h0_supply - dh ! Although this is a difference the result will - ! be non-negative because of the conditional. - ! Record the sub-cell start/end index that span the target cell i1. - itgt_start(i1) = i_start1 - itgt_end(i1) = i_sub - i_start1 = i_sub + 1 - ! Move the target index. - if (i1 < n1) then - i1 = i1 + 1 - h1_supply = h1(i1) - else - h1_supply = 0. - tgt_has_volume = .false. - endif - elseif (src_has_volume) then - ! We ran out of target volume but still have source cells to consume - h_sub(i_sub) = h0_supply - ! Record the sub-cell start/end index that span the source cell i0. - isrc_start(i0) = i_start0 - isrc_end(i0) = i_sub - i_start0 = i_sub + 1 - ! Record the sub-cell that is the largest fraction of the source cell. - isrc_max(i0) = i_max - i_max = i_sub + 1 - dh_max = 0. - ! Record the source cell thickness found by summing the sub-cell thicknesses. - h0_eff(i0) = dh0_eff - if (i0 < n0) then - i0 = i0 + 1 - h0_supply = h0(i0) - dh0_eff = 0. - else - h0_supply = 0. - src_has_volume = .false. - endif - elseif (tgt_has_volume) then - ! We ran out of source volume but still have target cells to consume - h_sub(i_sub) = h1_supply - ! Record the sub-cell start/end index that span the target cell i1. - itgt_start(i1) = i_start1 - itgt_end(i1) = i_sub - i_start1 = i_sub + 1 - ! Move the target index. - if (i1 < n1) then - i1 = i1 + 1 - h1_supply = h1(i1) - else - h1_supply = 0. - tgt_has_volume = .false. - endif - else - stop 'remap_via_sub_cells: THIS SHOULD NEVER HAPPEN!' - endif - - enddo + ! Calculate sub-layer thicknesses and indices connecting sub-layers to source and target grids + call intersect_src_tgt_grids( n0, h0, n1, h1, h_sub, h0_eff, & + isrc_start, isrc_end, isrc_max, itgt_start, itgt_end, isub_src ) ! Loop over each sub-cell to calculate average/integral values within each sub-cell. + ! Uses: h_sub, isub_src, h0_eff + ! Sets: u_sub, uh_sub xa = 0. dh0_eff = 0. uh_sub(1) = 0. @@ -712,6 +586,8 @@ subroutine remap_via_sub_cells( n0, h0, u0, ppoly0_E, ppoly0_coefs, n1, h1, meth ! Loop over each source cell substituting the integral/average for the thickest sub-cell (within ! the source cell) with the residual of the source cell integral minus the other sub-cell integrals ! aka a genius algorithm for accurate conservation when remapping from Robert Hallberg (@Hallberg-NOAA). + ! Uses: i0_last_thick_cell, isrc_max, h_sub, isrc_start, isrc_end, uh_sub, u0, h0 + ! Updates: uh_sub do i0 = 1, i0_last_thick_cell i_max = isrc_max(i0) dh_max = h_sub(i_max) @@ -728,6 +604,8 @@ subroutine remap_via_sub_cells( n0, h0, u0, ppoly0_E, ppoly0_coefs, n1, h1, meth endif ! Loop over each target cell summing the integrals from sub-cells within the target cell. + ! Uses: itgt_start, itgt_end, h_sub, uh_sub, u_sub + ! Sets: u1 uh_err = 0. do i1 = 1, n1 if (h1(i1) > 0.) then @@ -847,6 +725,190 @@ subroutine remap_via_sub_cells( n0, h0, u0, ppoly0_E, ppoly0_coefs, n1, h1, meth end subroutine remap_via_sub_cells +!> Returns the intersection of source and targets grids along with and auxiliary lists or indices. +!! +!! For source grid with thicknesses h0(1:n0) and target grid with thicknesses h1(1:n1) the intersection +!! or "subgrid" has thicknesses h_sub(1:n0+n1+1). +!! h0 and h1 must have the same units. h_sub will return with the same units as h0 and h1. +!! +!! Notes on the algorithm: +!! Internally, grids are defined by the interfaces (although we describe grids via thicknesses for accuracy). +!! The intersection or union of two grids is thus defined by the super set of both lists of interfaces. +!! Because both source and target grids can contain vanished cells, we do not eliminate repeated +!! interfaces from the union. +!! That is, the total number of interfaces of the sub-cells is equal to the total numer of interfaces of +!! the source grid (n0+1) plus the total number of interfaces of the target grid (n1+1), i.e. n0+n1+2. +!! Whenever target and source interfaces align, then the retention of identical interfaces leads to a +!! vanished subcell. +!! The remapping uses a common point of reference to the left (top) so there is always a vanished subcell +!! at the left (top). +!! If the total column thicknesses are the same, then the right (bottom) interfaces are also aligned and +!! so the last subcell will also be vanished. +subroutine intersect_src_tgt_grids( n0, h0, n1, h1, h_sub, h0_eff, & + isrc_start, isrc_end, isrc_max, itgt_start, itgt_end, isub_src ) + integer, intent(in) :: n0 !< Number of cells in source grid + real, intent(in) :: h0(n0) !< Source grid widths (size n0) [H] + integer, intent(in) :: n1 !< Number of cells in target grid + real, intent(in) :: h1(n1) !< Target grid widths (size n1) [H] + real, intent(out) :: h_sub(n0+n1+1) !< Overlapping sub-cell thicknesses, h_sub [H] + real, intent(out) :: h0_eff(n0) !< Effective thickness of source cells [H] + integer, intent(out) :: isrc_start(n0) !< Index of first sub-cell within each source cell + integer, intent(out) :: isrc_end(n0) !< Index of last sub-cell within each source cell + integer, intent(out) :: isrc_max(n0) !< Index of thickest sub-cell within each source cell + integer, intent(out) :: itgt_start(n1) !< Index of first sub-cell within each target cell + integer, intent(out) :: itgt_end(n1) !< Index of last sub-cell within each target cell + integer, intent(out) :: isub_src(n0+n1+1) !< Index of source cell for each sub-cell + ! Local variables + integer :: i_sub ! Index of sub-cell + integer :: i0 ! Index into h0(1:n0), source column + integer :: i1 ! Index into h1(1:n1), target column + integer :: i_start0 ! Used to record which sub-cells map to source cells + integer :: i_start1 ! Used to record which sub-cells map to target cells + integer :: i_max ! Used to record which sub-cell is the largest contribution of a source cell + real :: dh_max ! Used to record which sub-cell is the largest contribution of a source cell [H] + real :: h0_supply, h1_supply ! The amount of width available for constructing sub-cells [H] + real :: dh ! The width of the sub-cell [H] + real :: dh0_eff ! Running sum of source cell thickness [H] + ! For error checking/debugging + logical, parameter :: adjust_thickest_subcell = .true. ! To fix round-off conservation issues + logical, parameter :: debug_bounds = .false. ! For debugging overshoots etc. + integer :: i0_last_thick_cell + logical :: src_has_volume !< True if h0 has not been consumed + logical :: tgt_has_volume !< True if h1 has not been consumed + + i0_last_thick_cell = 0 + do i0 = 1, n0 + if (h0(i0)>0.) i0_last_thick_cell = i0 + enddo + + ! Initialize algorithm + h0_supply = h0(1) + h1_supply = h1(1) + src_has_volume = .true. + tgt_has_volume = .true. + i0 = 1 ; i1 = 1 + i_start0 = 1 ; i_start1 = 1 + i_max = 1 + dh_max = 0. + dh0_eff = 0. + + ! First sub-cell is always vanished + h_sub(1) = 0. + isrc_start(1) = 1 + isrc_end(1) = 1 + isrc_max(1) = 1 + isub_src(1) = 1 + + ! Loop over each sub-cell to calculate intersections with source and target grids + do i_sub = 2, n0+n1+1 + + ! This is the width of the sub-cell, determined by which ever column has the least + ! supply available to consume. + dh = min(h0_supply, h1_supply) + + ! This is the running sum of the source cell thickness. After summing over each + ! sub-cell, the sum of sub-cell thickness might differ from the original source + ! cell thickness due to round off. + dh0_eff = dh0_eff + min(dh, h0_supply) + + ! Record the source index (i0) that this sub-cell integral belongs to. This + ! is needed to index the reconstruction coefficients for the source cell + ! used in the integrals of the sub-cell width. + isub_src(i_sub) = i0 + h_sub(i_sub) = dh + + ! For recording the largest sub-cell within a source cell. + if (dh >= dh_max) then + i_max = i_sub + dh_max = dh + endif + + ! Which ever column (source or target) has the least width left to consume determined + ! the width, dh, of sub-cell i_sub in the expression for dh above. + if (h0_supply <= h1_supply .and. src_has_volume) then + ! h0_supply is smaller than h1_supply) so we consume h0_supply and increment the + ! source cell index. + h1_supply = h1_supply - dh ! Although this is a difference the result will + ! be non-negative because of the conditional. + ! Record the sub-cell start/end index that span the source cell i0. + isrc_start(i0) = i_start0 + isrc_end(i0) = i_sub + i_start0 = i_sub + 1 + ! Record the sub-cell that is the largest fraction of the source cell. + isrc_max(i0) = i_max + i_max = i_sub + 1 + dh_max = 0. + ! Record the source cell thickness found by summing the sub-cell thicknesses. + h0_eff(i0) = dh0_eff + ! Move the source index. + if (i0 < n0) then + i0 = i0 + 1 + h0_supply = h0(i0) + dh0_eff = 0. + else + h0_supply = 0. + src_has_volume = .false. + endif + elseif (h0_supply >= h1_supply .and. tgt_has_volume) then + ! h1_supply is smaller than h0_supply) so we consume h1_supply and increment the + ! target cell index. + h0_supply = h0_supply - dh ! Although this is a difference the result will + ! be non-negative because of the conditional. + ! Record the sub-cell start/end index that span the target cell i1. + itgt_start(i1) = i_start1 + itgt_end(i1) = i_sub + i_start1 = i_sub + 1 + ! Move the target index. + if (i1 < n1) then + i1 = i1 + 1 + h1_supply = h1(i1) + else + h1_supply = 0. + tgt_has_volume = .false. + endif + elseif (src_has_volume) then + ! We ran out of target volume but still have source cells to consume + h_sub(i_sub) = h0_supply + ! Record the sub-cell start/end index that span the source cell i0. + isrc_start(i0) = i_start0 + isrc_end(i0) = i_sub + i_start0 = i_sub + 1 + ! Record the sub-cell that is the largest fraction of the source cell. + isrc_max(i0) = i_max + i_max = i_sub + 1 + dh_max = 0. + ! Record the source cell thickness found by summing the sub-cell thicknesses. + h0_eff(i0) = dh0_eff + if (i0 < n0) then + i0 = i0 + 1 + h0_supply = h0(i0) + dh0_eff = 0. + else + h0_supply = 0. + src_has_volume = .false. + endif + elseif (tgt_has_volume) then + ! We ran out of source volume but still have target cells to consume + h_sub(i_sub) = h1_supply + ! Record the sub-cell start/end index that span the target cell i1. + itgt_start(i1) = i_start1 + itgt_end(i1) = i_sub + i_start1 = i_sub + 1 + ! Move the target index. + if (i1 < n1) then + i1 = i1 + 1 + h1_supply = h1(i1) + else + h1_supply = 0. + tgt_has_volume = .false. + endif + else + stop 'remap_via_sub_cells: THIS SHOULD NEVER HAPPEN!' + endif + + enddo +end subroutine intersect_src_tgt_grids + !> Linearly interpolate interface data, u_src, from grid h_src to a grid h_dest subroutine interpolate_column(nsrc, h_src, u_src, ndest, h_dest, u_dest, mask_edges) integer, intent(in) :: nsrc !< Number of source cells @@ -1348,6 +1410,8 @@ logical function remapping_unit_tests(verbose) real, allocatable, dimension(:,:) :: ppoly0_E ! Edge values of polynomials [A] real, allocatable, dimension(:,:) :: ppoly0_S ! Edge slopes of polynomials [A H-1] real, allocatable, dimension(:,:) :: ppoly0_coefs ! Coefficients of polynomials [A] + real, allocatable, dimension(:) :: h_sub, h0_eff ! Subgrid and effective source thicknesses [H] + integer, allocatable, dimension(:) :: isrc_start, isrc_end, isrc_max, itgt_start, itgt_end, isub_src ! Indices integer :: answer_date ! The vintage of the expressions to test integer :: i real, parameter :: hNeglect_dflt = 1.0e-30 ! A thickness [H ~> m or kg m-2] that can be @@ -1435,11 +1499,14 @@ logical function remapping_unit_tests(verbose) 3, (/2.25,1.5,1./), INTEGRATION_PPM, .false., u2, err ) if (verbose) call dumpGrid(3,h2,x2,u2) - if (.not. remapping_unit_tests) write(stdout,*) 'Pass' + deallocate(ppoly0_E, ppoly0_S, ppoly0_coefs) + + ! =============================================== + ! This section tests the reconstruction functions + ! =============================================== write(stdout,*) '===== MOM_remapping: new remapping_unit_tests ==================' - deallocate(ppoly0_E, ppoly0_S, ppoly0_coefs) allocate(ppoly0_coefs(5,6)) allocate(ppoly0_E(5,2)) allocate(ppoly0_S(5,2)) @@ -1552,6 +1619,10 @@ logical function remapping_unit_tests(verbose) remapping_unit_tests = remapping_unit_tests .or. & test_answer(v, 5, ppoly0_coefs(:,3), (/0.,-3.,3.,0.,0./), 'Limits PPM: P2') + ! ============================================================== + ! This section tests the components of remapping_via_sub_cells() + ! ============================================================== + call PLM_reconstruction(4, (/0.,1.,1.,0./), (/5.,4.,2.,1./), ppoly0_E(1:4,:), & ppoly0_coefs(1:4,:), h_neglect ) remapping_unit_tests = remapping_unit_tests .or. & @@ -1566,10 +1637,225 @@ logical function remapping_unit_tests(verbose) deallocate(ppoly0_E, ppoly0_S, ppoly0_coefs) - ! This line carries out tests on some older remapping schemes. - remapping_unit_tests = remapping_unit_tests .or. remapping_attic_unit_tests(verbose) + ! Test 1: n0=2, n1=3 Maps uniform grids with one extra target layer and no implicitly-vanished interior sub-layers + ! h_src = | 3 | 3 | + ! h_tgt = | 2 | 2 | 2 | + ! h_sub = |0| 2 | 1 | 1 | 2 |0| + ! isrc_start |1 | 4 | + ! isrc_end | 3 | 5 | + ! isrc_max | 2 | 5 | + ! itgt_start |1 | 3 | 5 | + ! itgt_end | 2 | 4 | 6| + ! isub_src |1| 1 | 1 | 2 | 2 |2| + 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) ) + call intersect_src_tgt_grids( 2, (/3., 3./), & ! n0, h0 + 3, (/2., 2., 2./), & ! n1, h1 + h_sub, h0_eff, & + isrc_start, isrc_end, isrc_max, itgt_start, itgt_end, isub_src ) + write(stdout,*) "intersect_src_tgt_grids test 1: n0=2, n1=3" + write(stdout,*) " h_src = | 3 | 3 |" + write(stdout,*) " h_tgt = | 2 | 2 | 2 |" + remapping_unit_tests = remapping_unit_tests .or. & + test_answer(v, 6, h_sub, (/0.,2.,1.,1.,2.,0./), 'h_sub') + remapping_unit_tests = remapping_unit_tests .or. & + test_answer(v, 2, h0_eff, (/3.,3./), 'h0_eff') + remapping_unit_tests = remapping_unit_tests .or. & + test_ianswer(v, 2, isrc_start, (/1,4/), 'isrc_start') + remapping_unit_tests = remapping_unit_tests .or. & + test_ianswer(v, 2, isrc_end, (/3,5/), 'isrc_end') + remapping_unit_tests = remapping_unit_tests .or. & + test_ianswer(v, 2, isrc_max, (/2,5/), 'isrc_max') + remapping_unit_tests = remapping_unit_tests .or. & + test_ianswer(v, 3, itgt_start, (/1,3,5/), 'itgt_start') + remapping_unit_tests = remapping_unit_tests .or. & + test_ianswer(v, 3, itgt_end, (/2,4,6/), 'itgt_end') + remapping_unit_tests = remapping_unit_tests .or. & + test_ianswer(v, 6, isub_src, (/1,1,1,2,2,2/), 'isub_src') + deallocate( h_sub, h0_eff, isrc_start, isrc_end, isrc_max, itgt_start, itgt_end, isub_src ) + + ! Test 2: n0=3, n1=2 Reverses "test 1" with more source than target layers + ! h_src = | 2 | 2 | 2 | + ! h_tgt = | 3 | 3 | + ! h_sub = |0| 2 | 1 | 1 | 2 |0| + ! isrc_start |1 | 3 | 5 | + ! isrc_end | 2 | 4 | 5 | + ! isrc_max | 2 | 4 | 5 | + ! itgt_start |1 | 4 | + ! itgt_end | 3 | 6| + ! isub_src |1| 1 | 2 | 2 | 3 |3| + 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) ) + call intersect_src_tgt_grids( 3, (/2., 2., 2./), & ! n0, h0 + 2, (/3., 3./), & ! n1, h1 + h_sub, h0_eff, & + isrc_start, isrc_end, isrc_max, itgt_start, itgt_end, isub_src ) + write(stdout,*) "intersect_src_tgt_grids test 2: n0=3, n1=2" + write(stdout,*) " h_src = | 2 | 2 | 2 |" + write(stdout,*) " h_tgt = | 3 | 3 |" + remapping_unit_tests = remapping_unit_tests .or. & + test_answer(v, 6, h_sub, (/0.,2.,1.,1.,2.,0./), 'h_sub') + remapping_unit_tests = remapping_unit_tests .or. & + test_answer(v, 3, h0_eff, (/2.,2.,2./), 'h0_eff') + remapping_unit_tests = remapping_unit_tests .or. & + test_ianswer(v, 3, isrc_start, (/1,3,5/), 'isrc_start') + remapping_unit_tests = remapping_unit_tests .or. & + test_ianswer(v, 3, isrc_end, (/2,4,5/), 'isrc_end') + remapping_unit_tests = remapping_unit_tests .or. & + test_ianswer(v, 3, isrc_max, (/2,4,5/), 'isrc_max') + remapping_unit_tests = remapping_unit_tests .or. & + test_ianswer(v, 2, itgt_start, (/1,4/), 'itgt_start') + remapping_unit_tests = remapping_unit_tests .or. & + test_ianswer(v, 2, itgt_end, (/3,6/), 'itgt_end') + remapping_unit_tests = remapping_unit_tests .or. & + test_ianswer(v, 6, isub_src, (/1,1,2,2,3,3/), 'isub_src') + deallocate( h_sub, h0_eff, isrc_start, isrc_end, isrc_max, itgt_start, itgt_end, isub_src ) + + ! Test 3: n0=2, n1=3 With aligned interfaces that lead to implicitly-vanished interior sub-layers + ! h_src = | 2 | 4 | + ! h_tgt = | 2 | 2 | 2 | + ! h_sub = |0| 2 |0| 2 | 2 |0| + ! isrc_start |1 |3 | + ! isrc_end | 2 | 5 | + ! isrc_max | 2 | 5 | + ! itgt_start |1 | 4 | 5 | + ! itgt_end | 3| 4 | 6| + ! isub_src |1| 1 |2| 2 | 2 |2| + 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) ) + call intersect_src_tgt_grids( 2, (/2., 4./), & ! n0, h0 + 3, (/2., 2., 2./), & ! n1, h1 + h_sub, h0_eff, & + isrc_start, isrc_end, isrc_max, itgt_start, itgt_end, isub_src ) + write(stdout,*) "intersect_src_tgt_grids test 3: n0=2, n1=3" + write(stdout,*) " h_src = | 2 | 4 |" + write(stdout,*) " h_tgt = | 2 | 2 | 2 |" + remapping_unit_tests = remapping_unit_tests .or. & + test_answer(v, 6, h_sub, (/0.,2.,0.,2.,2.,0./), 'h_sub') + remapping_unit_tests = remapping_unit_tests .or. & + test_answer(v, 2, h0_eff, (/2.,4./), 'h0_eff') + remapping_unit_tests = remapping_unit_tests .or. & + test_ianswer(v, 2, isrc_start, (/1,3/), 'isrc_start') + remapping_unit_tests = remapping_unit_tests .or. & + test_ianswer(v, 2, isrc_end, (/2,5/), 'isrc_end') + remapping_unit_tests = remapping_unit_tests .or. & + test_ianswer(v, 2, isrc_max, (/2,5/), 'isrc_max') + remapping_unit_tests = remapping_unit_tests .or. & + test_ianswer(v, 3, itgt_start, (/1,4,5/), 'itgt_start') + remapping_unit_tests = remapping_unit_tests .or. & + test_ianswer(v, 3, itgt_end, (/3,4,6/), 'itgt_end') + remapping_unit_tests = remapping_unit_tests .or. & + test_ianswer(v, 6, isub_src, (/1,1,2,2,2,2/), 'isub_src') + deallocate( h_sub, h0_eff, isrc_start, isrc_end, isrc_max, itgt_start, itgt_end, isub_src ) + + ! Test 4: n0=2, n1=3 Incomplete target column, sum(h_tgt)sum(h_src), useful for diagnostics + ! h_src = | 2 | 2 | 1 | + ! h_tgt = | 2 | 4 | + ! h_sub = |0| 2 |0| 2 | 1 | 1 | + ! isrc_start |1 |3 | 5 | + ! isrc_end | 2 | 4 | 5 | + ! isrc_max | 2 | 4 | 5 | | + ! itgt_start |1 | 4 | + ! itgt_end | 3| 6 | + ! isub_src |1| 1 |2| 2 | 3 | 3 | + 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) ) + call intersect_src_tgt_grids( 3, (/2., 2., 1./), & ! n0, h0 + 2, (/2., 4./), & ! n1, h1 + h_sub, h0_eff, & + isrc_start, isrc_end, isrc_max, itgt_start, itgt_end, isub_src ) + write(stdout,*) "intersect_src_tgt_grids test 5: n0=3, n1=2" + write(stdout,*) " h_src = | 2 | 2 | 1 |" + write(stdout,*) " h_tgt = | 2 | 4 |" + remapping_unit_tests = remapping_unit_tests .or. & + test_answer(v, 6, h_sub, (/0.,2.,0.,2.,1.,1./), 'h_sub') + remapping_unit_tests = remapping_unit_tests .or. & + test_answer(v, 3, h0_eff, (/2.,2.,1./), 'h0_eff') + remapping_unit_tests = remapping_unit_tests .or. & + test_ianswer(v, 3, isrc_start, (/1,3,5/), 'isrc_start') + remapping_unit_tests = remapping_unit_tests .or. & + test_ianswer(v, 3, isrc_end, (/2,4,5/), 'isrc_end') + remapping_unit_tests = remapping_unit_tests .or. & + test_ianswer(v, 3, isrc_max, (/2,4,5/), 'isrc_max') + remapping_unit_tests = remapping_unit_tests .or. & + test_ianswer(v, 2, itgt_start, (/1,4/), 'itgt_start') + remapping_unit_tests = remapping_unit_tests .or. & + test_ianswer(v, 2, itgt_end, (/3,6/), 'itgt_end') + remapping_unit_tests = remapping_unit_tests .or. & + test_ianswer(v, 6, isub_src, (/1,1,2,2,3,3/), 'isub_src') + deallocate( h_sub, h0_eff, isrc_start, isrc_end, isrc_max, itgt_start, itgt_end, isub_src ) + + ! Test 6: n0=3, n1=5 Source and targets with vanished layers + ! h_src = | 2 |0| 2 | + ! h_tgt = | 1 |0| 1 |0| 2 | + ! h_sub = |0| 1 |0| 1 |0|0|0| 2 |0| + ! isrc_start |1 |5|6 | + ! isrc_end | 4 |5| 8 | + ! isrc_max | 4 |5| 8 | + ! itgt_start |1 |3| 4 |7| 8 | + ! itgt_end | 2 |3| 6|7| 9| + ! isub_src |1| 1 |1| 1 |2|3|3| 3 |3| + allocate( h_sub(9), h0_eff(3), isrc_start(3), isrc_end(3), isrc_max(3), itgt_start(5), itgt_end(5), isub_src(9) ) + call intersect_src_tgt_grids( 3, (/2., 0., 2./), & ! n0, h0 + 5, (/1., 0., 1., 0., 2./), & ! n1, h1 + h_sub, h0_eff, & + isrc_start, isrc_end, isrc_max, itgt_start, itgt_end, isub_src ) + write(stdout,*) "intersect_src_tgt_grids test 6: n0=3, n1=5" + write(stdout,*) " h_src = | 2 |0| 2 |" + write(stdout,*) " h_tgt = | 1 |0| 1 |0| 2 |" + remapping_unit_tests = remapping_unit_tests .or. & + test_answer(v, 9, h_sub, (/0.,1.,0.,1.,0.,0.,0.,2.,0./), 'h_sub') + remapping_unit_tests = remapping_unit_tests .or. & + test_answer(v, 3, h0_eff, (/2.,0.,2./), 'h0_eff') + remapping_unit_tests = remapping_unit_tests .or. & + test_ianswer(v, 3, isrc_start, (/1,5,6/), 'isrc_start') + remapping_unit_tests = remapping_unit_tests .or. & + test_ianswer(v, 3, isrc_end, (/4,5,8/), 'isrc_end') + remapping_unit_tests = remapping_unit_tests .or. & + test_ianswer(v, 3, isrc_max, (/4,5,8/), 'isrc_max') + remapping_unit_tests = remapping_unit_tests .or. & + test_ianswer(v, 5, itgt_start, (/1,3,4,7,8/), 'itgt_start') + remapping_unit_tests = remapping_unit_tests .or. & + test_ianswer(v, 5, itgt_end, (/2,3,6,7,9/), 'itgt_end') + remapping_unit_tests = remapping_unit_tests .or. & + test_ianswer(v, 9, isub_src, (/1,1,1,1,2,3,3,3,3/), 'isub_src') + deallocate( h_sub, h0_eff, isrc_start, isrc_end, isrc_max, itgt_start, itgt_end, isub_src ) - if (.not. remapping_unit_tests) write(stdout,*) 'Pass' + ! ============================================================ + ! This section tests interpolation and reintegration functions + ! ============================================================ write(stdout,*) '=== MOM_remapping: interpolation and reintegration unit tests ===' if (verbose) write(stdout,*) '- - - - - - - - - - interpolation tests - - - - - - - - -' @@ -1671,6 +1957,9 @@ logical function remapping_unit_tests(verbose) 3, (/0.,0.,0./), (/0., 0., 0./) ) remapping_unit_tests = remapping_unit_tests .or. fail + ! This line carries out tests on some older remapping schemes. + remapping_unit_tests = remapping_unit_tests .or. remapping_attic_unit_tests(verbose) + if (.not. remapping_unit_tests) write(stdout,*) 'Pass' end function remapping_unit_tests @@ -1706,6 +1995,33 @@ logical function test_answer(verbose, n, u, u_true, label, tol) end function test_answer +!> Returns true if any cell of i and i_true are not identical. Returns false otherwise. +logical function test_ianswer(verbose, n, i_calc, i_true, label) + logical, intent(in) :: verbose !< If true, write results to stdout + integer, intent(in) :: n !< Number of cells in u + integer, dimension(n), intent(in) :: i_calc !< Values to test [A] + integer, dimension(n), intent(in) :: i_true !< Values to test against (correct answer) [A] + character(len=*), intent(in) :: label !< Message + ! Local variables + integer :: k + + test_ianswer = .false. + do k = 1, n + if (i_calc(k) /= i_true(k)) test_ianswer = .true. + enddo + if (verbose) then + write(stdout,'(a12," : calculated =",30i3)') label,i_calc + write(stdout,'(12x," correct =",30i3)') i_true + if (test_ianswer) write(stdout,'("FAIL ->",5x," error =",30i3)') i_calc(:) - i_true(:) + endif + if (test_ianswer) then + write(stderr,'(a12," : calculated =",30i3)') label,i_calc + write(stderr,'(12x," correct =",30i3)') i_true + write(stderr,'("FAIL ->",5x," error =",30i3)') i_calc(:) - i_true(:) + endif + +end function test_ianswer + !> Returns true if a test of interpolate_column() produces the wrong answer logical function test_interp(verbose, msg, nsrc, h_src, u_src, ndest, h_dest, u_true) logical, intent(in) :: verbose !< If true, write results to stdout