mom_remapping Module Reference

Detailed Description

Provides column-wise vertical remapping functions.

Data Types
type	remapping_cs
	Container for remapping parameters. More...

Functions/Subroutines
subroutine, public	remapping_set_param (CS, remapping_scheme, boundary_extrapolation, check_reconstruction, check_remapping, force_bounds_in_subcell, answers_2018)
	Set parameters within remapping object. More...
subroutine, public	extract_member_remapping_cs (CS, remapping_scheme, degree, boundary_extrapolation, check_reconstruction, check_remapping, force_bounds_in_subcell)
subroutine	buildgridfromh (nz, h, x)
	Calculate edge coordinate x from cell width h. More...
logical function	ispossumerrsignificant (n1, sum1, n2, sum2)
	Compare two summation estimates of positive data and judge if due to more than round-off. When two sums are calculated from different vectors that should add up to the same value, the results can differ by round off. The round off error can be bounded to be proportional to the number of operations. This function returns true if the difference between sum1 and sum2 is larger than than the estimated round off bound. More...
subroutine, public	remapping_core_h (CS, n0, h0, u0, n1, h1, u1, h_neglect, h_neglect_edge)
	Remaps column of values u0 on grid h0 to grid h1 assuming the top edge is aligned. More...
subroutine, public	remapping_core_w (CS, n0, h0, u0, n1, dx, u1, h_neglect, h_neglect_edge)
	Remaps column of values u0 on grid h0 to implied grid h1 where the interfaces of h1 differ from those of h0 by dx. More...
subroutine, public	build_reconstructions_1d (CS, n0, h0, u0, ppoly_r_coefs, ppoly_r_E, ppoly_r_S, iMethod, h_neglect, h_neglect_edge)
	Creates polynomial reconstructions of u0 on the source grid h0. More...
subroutine	check_reconstructions_1d (n0, h0, u0, deg, boundary_extrapolation, ppoly_r_coefs, ppoly_r_E, ppoly_r_S)
	Checks that edge values and reconstructions satisfy bounds. More...
subroutine	remap_via_sub_cells (n0, h0, u0, ppoly0_E, ppoly0_coefs, n1, h1, method, force_bounds_in_subcell, u1, uh_err, ah_sub, aisub_src, aiss, aise)
	Remaps column of n0 values u0 on grid h0 to grid h1 with n1 cells by calculating the n0+n1+1 sub-integrals of the intersection of h0 and h1, and the summing the appropriate integrals into the h1*u1 values. h0 and h1 must have the same units. More...
real function, public	average_value_ppoly (n0, u0, ppoly0_E, ppoly0_coefs, method, i0, xa, xb)
	Returns the average value of a reconstruction within a single source cell, i0, between the non-dimensional positions xa and xb (xa<=xb) with dimensional separation dh. More...
subroutine	measure_input_bounds (n0, h0, u0, edge_values, h0tot, h0err, u0tot, u0err, u0min, u0max)
	Measure totals and bounds on source grid. More...
subroutine	measure_output_bounds (n1, h1, u1, h1tot, h1err, u1tot, u1err, u1min, u1max)
	Measure totals and bounds on destination grid. More...
subroutine	remapbyprojection (n0, h0, u0, ppoly0_E, ppoly0_coefs, n1, h1, method, u1, h_neglect)
	Remaps column of values u0 on grid h0 to grid h1 by integrating over the projection of each h1 cell onto the h0 grid. More...
subroutine	remapbydeltaz (n0, h0, u0, ppoly0_E, ppoly0_coefs, n1, dx1, method, u1, h1, h_neglect)
	Remaps column of values u0 on grid h0 to implied grid h1 where the interfaces of h1 differ from those of h0 by dx. The new grid is defined relative to the original grid by change dx1(:) = xNew(:) - xOld(:) and the remapping calculated so that hNew(k) qNew(k) = hOld(k) qOld(k) + F(k+1) - F(k) where F(k) = dx1(k) qAverage and where qAverage is the average qOld in the region zOld(k) to zNew(k). More...
subroutine	integraterecononinterval (n0, h0, u0, ppoly0_E, ppoly0_coefs, method, xL, xR, hC, uAve, jStart, xStart, h_neglect)
	Integrate the reconstructed column profile over a single cell. More...
subroutine, public	dzfromh1h2 (n1, h1, n2, h2, dx)
	Calculates the change in interface positions based on h1 and h2. More...
subroutine, public	initialize_remapping (CS, remapping_scheme, boundary_extrapolation, check_reconstruction, check_remapping, force_bounds_in_subcell, answers_2018)
	Constructor for remapping control structure. More...
subroutine	setreconstructiontype (string, CS)
	Changes the method of reconstruction Use this routine to parse a string parameter specifying the reconstruction and re-allocates work arrays appropriately. It is called from initialize_remapping but can be called from an external module too. More...
subroutine, public	end_remapping (CS)
	Destrcutor for remapping control structure. More...
logical function, public	remapping_unit_tests (verbose)
	Runs unit tests on remapping functions. Should only be called from a single/root thread Returns True if a test fails, otherwise False. More...
logical function	test_answer (verbose, n, u, u_true, label, tol)
	Returns true if any cell of u and u_true are not identical. Returns false otherwise. More...
subroutine	dumpgrid (n, h, x, u)
	Convenience function for printing grid to screen. More...

Variables
integer, parameter	remapping_pcm = 0
	O(h^1) remapping scheme.
integer, parameter	remapping_plm = 1
	O(h^2) remapping scheme.
integer, parameter	remapping_ppm_h4 = 2
	O(h^3) remapping scheme.
integer, parameter	remapping_ppm_ih4 = 3
	O(h^3) remapping scheme.
integer, parameter	remapping_pqm_ih4ih3 = 4
	O(h^4) remapping scheme.
integer, parameter	remapping_pqm_ih6ih5 = 5
	O(h^5) remapping scheme.
integer, parameter	integration_pcm = 0
	Piecewise Constant Method.
integer, parameter	integration_plm = 1
	Piecewise Linear Method.
integer, parameter	integration_ppm = 3
	Piecewise Parabolic Method.
integer, parameter	integration_pqm = 5
	Piecewise Quartic Method.
character(len=40)	mdl = "MOM_remapping"
	This module's name.
character(len=256), public	remappingschemesdoc = "PCM (1st-order accurate)\n"// "PLM (2nd-order accurate)\n"// "PPM_H4 (3rd-order accurate)\n"// "PPM_IH4 (3rd-order accurate)\n"// "PQM_IH4IH3 (4th-order accurate)\n"// "PQM_IH6IH5 (5th-order accurate)\n"
	Documentation for external callers.
character(len=3), public	remappingdefaultscheme = "PLM"
	Default remapping method.
real, parameter	hneglect_dflt = 1.E-30
	A thickness [H ~> m or kg m-2] that can be added to thicknesses in a denominator without changing the numerical result, except where a division by zero would otherwise occur.
logical, parameter	old_algorithm = .false.
	Use the old "broken" algorithm. This is a temporary measure to assist debugging until we delete the old algorithm.

Function/Subroutine Documentation

average_value_ppoly()

real function, public mom_remapping::average_value_ppoly	(	integer, intent(in)	n0,
		real, dimension(:), intent(in)	u0,
		real, dimension(:,:), intent(in)	ppoly0_E,
		real, dimension(:,:), intent(in)	ppoly0_coefs,
		integer, intent(in)	method,
		integer, intent(in)	i0,
		real, intent(in)	xa,
		real, intent(in)	xb
	)

Returns the average value of a reconstruction within a single source cell, i0, between the non-dimensional positions xa and xb (xa<=xb) with dimensional separation dh.

Parameters

[in]	n0	Number of cells in source grid
[in]	u0	Cell means
[in]	ppoly0_e	Edge value of polynomial
[in]	ppoly0_coefs	Coefficients of polynomial
[in]	method	Remapping scheme to use
[in]	i0	Source cell index
[in]	xa	Non-dimensional start position within source cell
[in]	xb	Non-dimensional end position within source cell

Definition at line 929 of file MOM_remapping.F90.

  integer,       intent(in)    :: n0     !< Number of cells in source grid
  real,          intent(in)    :: u0(:)  !< Cell means
  real,          intent(in)    :: ppoly0_E(:,:)     !< Edge value of polynomial
  real,          intent(in)    :: ppoly0_coefs(:,:) !< Coefficients of polynomial
  integer,       intent(in)    :: method !< Remapping scheme to use
  integer,       intent(in)    :: i0     !< Source cell index
  real,          intent(in)    :: xa     !< Non-dimensional start position within source cell
  real,          intent(in)    :: xb     !< Non-dimensional end position within source cell
  ! Local variables
  real :: u_ave, xa_2, xb_2, xa2pxb2, xapxb
  real, parameter :: r_3 = 1.0/3.0 ! Used in evaluation of integrated polynomials
 
  real :: mx, a_L, a_R, u_c, Ya, Yb, my, xa2b2ab, Ya2b2ab, a_c
 
  if (xb > xa) then
    select case ( method )
      case ( integration_pcm )
        u_ave = u0(i0)
      case ( integration_plm )
        u_ave = (                                           &
            ppoly0_coefs(i0,1)                       &
          + ppoly0_coefs(i0,2) * 0.5 * ( xb + xa ) )
      case ( integration_ppm )
        mx = 0.5 * ( xa + xb )
        a_l = ppoly0_e(i0, 1)
        a_r = ppoly0_e(i0, 2)
        u_c = u0(i0)
        a_c = 0.5 * ( ( u_c - a_l ) + ( u_c - a_r ) ) ! a_6 / 6
        if (mx<0.5) then
          ! This integration of the PPM reconstruction is expressed in distances from the left edge
          xa2b2ab = (xa*xa+xb*xb)+xa*xb
          u_ave = a_l + ( ( a_r - a_l ) * mx &
                          + a_c * ( 3. * ( xb + xa ) - 2.*xa2b2ab ) )
        else
          ! This integration of the PPM reconstruction is expressed in distances from the right edge
          ya = 1. - xa
          yb = 1. - xb
          my = 0.5 * ( ya + yb )
          ya2b2ab = (ya*ya+yb*yb)+ya*yb
          u_ave = a_r  + ( ( a_l - a_r ) * my &
                           + a_c * ( 3. * ( yb + ya ) - 2.*ya2b2ab ) )
        endif
      case ( integration_pqm )
        xa_2 = xa*xa
        xb_2 = xb*xb
        xa2pxb2 = xa_2 + xb_2
        xapxb = xa + xb
        u_ave = (                                                                               &
              ppoly0_coefs(i0,1)                                                         &
          + ( ppoly0_coefs(i0,2) * 0.5 * ( xapxb )                                       &
          + ( ppoly0_coefs(i0,3) * r_3 * ( xa2pxb2 + xa*xb )                             &
          + ( ppoly0_coefs(i0,4) * 0.25* ( xa2pxb2 * xapxb )                             &
          +   ppoly0_coefs(i0,5) * 0.2 * ( ( xb*xb_2 + xa*xa_2 ) * xapxb + xa_2*xb_2 ) ) ) ) )
      case default
        call mom_error( fatal,'The selected integration method is invalid' )
    end select
  else ! dh == 0.
    select case ( method )
      case ( integration_pcm )
        u_ave =        ppoly0_coefs(i0,1)
      case ( integration_plm )
       !u_ave =        ppoly0_coefs(i0,1)   &
       !      + xa *   ppoly0_coefs(i0,2)
        a_l = ppoly0_e(i0, 1)
        a_r = ppoly0_e(i0, 2)
        ya = 1. - xa
        if (xa < 0.5) then
          u_ave = a_l + xa * ( a_r - a_l )
        else
          u_ave = a_r + ya * ( a_l - a_r )
        endif
      case ( integration_ppm )
       !u_ave =        ppoly0_coefs(i0,1)   &
       !      + xa * ( ppoly0_coefs(i0,2)   &
       !      + xa *   ppoly0_coefs(i0,3) )
        a_l = ppoly0_e(i0, 1)
        a_r = ppoly0_e(i0, 2)
        u_c = u0(i0)
        a_c = 3. * ( ( u_c - a_l ) + ( u_c - a_r ) ) ! a_6
        ya = 1. - xa
        if (xa < 0.5) then
          u_ave = a_l + xa * ( ( a_r - a_l ) + a_c * ya )
        else
          u_ave = a_r + ya * ( ( a_l - a_r ) + a_c * xa )
        endif
      case ( integration_pqm )
        u_ave =        ppoly0_coefs(i0,1)   &
              + xa * ( ppoly0_coefs(i0,2)   &
              + xa * ( ppoly0_coefs(i0,3)   &
              + xa * ( ppoly0_coefs(i0,4)   &
              + xa *   ppoly0_coefs(i0,5) ) ) )
      case default
        call mom_error( fatal,'The selected integration method is invalid' )
    end select
  endif
  average_value_ppoly = u_ave
 

build_reconstructions_1d()

subroutine, public mom_remapping::build_reconstructions_1d	(	type(remapping_cs), intent(in)	CS,
		integer, intent(in)	n0,
		real, dimension(n0), intent(in)	h0,
		real, dimension(n0), intent(in)	u0,
		real, dimension(n0,cs%degree+1), intent(out)	ppoly_r_coefs,
		real, dimension(n0,2), intent(out)	ppoly_r_E,
		real, dimension(n0,2), intent(out)	ppoly_r_S,
		integer, intent(out)	iMethod,
		real, intent(in), optional	h_neglect,
		real, intent(in), optional	h_neglect_edge
	)

Creates polynomial reconstructions of u0 on the source grid h0.

Parameters

[in]	cs	Remapping control structure
[in]	n0	Number of cells on source grid
[in]	h0	Cell widths on source grid
[in]	u0	Cell averages on source grid
[out]	ppoly_r_coefs	Coefficients of polynomial
[out]	ppoly_r_e	Edge value of polynomial
[out]	ppoly_r_s	Edge slope of polynomial
[out]	imethod	Integration method
[in]	h_neglect	A negligibly small width for the purpose of cell reconstructions in the same units as h0.
[in]	h_neglect_edge	A negligibly small width for the purpose of edge value calculations in the same units as h0.

Definition at line 355 of file MOM_remapping.F90.

  type(remapping_CS),    intent(in)  :: CS !< Remapping control structure
  integer,               intent(in)  :: n0 !< Number of cells on source grid
  real, dimension(n0),   intent(in)  :: h0 !< Cell widths on source grid
  real, dimension(n0),   intent(in)  :: u0 !< Cell averages on source grid
  real, dimension(n0,CS%degree+1), &
                         intent(out) :: ppoly_r_coefs !< Coefficients of polynomial
  real, dimension(n0,2), intent(out) :: ppoly_r_E !< Edge value of polynomial
  real, dimension(n0,2), intent(out) :: ppoly_r_S !< Edge slope of polynomial
  integer,               intent(out) :: iMethod !< Integration method
  real, optional,        intent(in)  :: h_neglect !< A negligibly small width for the
                                         !! purpose of cell reconstructions
                                         !! in the same units as h0.
  real, optional,        intent(in)  :: h_neglect_edge !< A negligibly small width
                                         !! for the purpose of edge value
                                         !! calculations in the same units as h0.
  ! Local variables
  integer :: local_remapping_scheme
  integer :: remapping_scheme !< Remapping scheme
  logical :: boundary_extrapolation !< Extrapolate at boundaries if true
 
  ! Reset polynomial
  ppoly_r_e(:,:) = 0.0
  ppoly_r_s(:,:) = 0.0
  ppoly_r_coefs(:,:) = 0.0
  imethod = -999
 
  local_remapping_scheme = cs%remapping_scheme
  if (n0<=1) then
    local_remapping_scheme = remapping_pcm
  elseif (n0<=3) then
    local_remapping_scheme = min( local_remapping_scheme, remapping_plm )
  elseif (n0<=4) then
    local_remapping_scheme = min( local_remapping_scheme, remapping_ppm_h4 )
  endif
  select case ( local_remapping_scheme )
    case ( remapping_pcm )
      call pcm_reconstruction( n0, u0, ppoly_r_e, ppoly_r_coefs)
      imethod = integration_pcm
    case ( remapping_plm )
      call plm_reconstruction( n0, h0, u0, ppoly_r_e, ppoly_r_coefs, h_neglect )
      if ( cs%boundary_extrapolation ) then
        call plm_boundary_extrapolation( n0, h0, u0, ppoly_r_e, ppoly_r_coefs, h_neglect)
      endif
      imethod = integration_plm
    case ( remapping_ppm_h4 )
      call edge_values_explicit_h4( n0, h0, u0, ppoly_r_e, h_neglect_edge, answers_2018=cs%answers_2018 )
      call ppm_reconstruction( n0, h0, u0, ppoly_r_e, ppoly_r_coefs, h_neglect, answers_2018=cs%answers_2018 )
      if ( cs%boundary_extrapolation ) then
        call ppm_boundary_extrapolation( n0, h0, u0, ppoly_r_e, ppoly_r_coefs, h_neglect )
      endif
      imethod = integration_ppm
    case ( remapping_ppm_ih4 )
      call edge_values_implicit_h4( n0, h0, u0, ppoly_r_e, h_neglect_edge, answers_2018=cs%answers_2018 )
      call ppm_reconstruction( n0, h0, u0, ppoly_r_e, ppoly_r_coefs, h_neglect, answers_2018=cs%answers_2018 )
      if ( cs%boundary_extrapolation ) then
        call ppm_boundary_extrapolation( n0, h0, u0, ppoly_r_e, ppoly_r_coefs, h_neglect )
      endif
      imethod = integration_ppm
    case ( remapping_pqm_ih4ih3 )
      call edge_values_implicit_h4( n0, h0, u0, ppoly_r_e, h_neglect_edge, answers_2018=cs%answers_2018 )
      call edge_slopes_implicit_h3( n0, h0, u0, ppoly_r_s, h_neglect, answers_2018=cs%answers_2018 )
      call pqm_reconstruction( n0, h0, u0, ppoly_r_e, ppoly_r_s, ppoly_r_coefs, h_neglect, &
                               answers_2018=cs%answers_2018 )
      if ( cs%boundary_extrapolation ) then
        call pqm_boundary_extrapolation_v1( n0, h0, u0, ppoly_r_e, ppoly_r_s, &
                                            ppoly_r_coefs, h_neglect )
      endif
      imethod = integration_pqm
    case ( remapping_pqm_ih6ih5 )
      call edge_values_implicit_h6( n0, h0, u0, ppoly_r_e, h_neglect_edge, answers_2018=cs%answers_2018 )
      call edge_slopes_implicit_h5( n0, h0, u0, ppoly_r_s, h_neglect, answers_2018=cs%answers_2018 )
      call pqm_reconstruction( n0, h0, u0, ppoly_r_e, ppoly_r_s, ppoly_r_coefs, h_neglect, &
                               answers_2018=cs%answers_2018 )
      if ( cs%boundary_extrapolation ) then
        call pqm_boundary_extrapolation_v1( n0, h0, u0, ppoly_r_e, ppoly_r_s, &
                                            ppoly_r_coefs, h_neglect )
      endif
      imethod = integration_pqm
    case default
      call mom_error( fatal, 'MOM_remapping, build_reconstructions_1d: '//&
           'The selected remapping method is invalid' )
  end select
 

buildgridfromh()

subroutine mom_remapping::buildgridfromh ( integer, intent(in)  nz, real, dimension(nz), intent(in)  h, real, dimension(nz+1), intent(inout)  x  )

private

Calculate edge coordinate x from cell width h.

Parameters

[in]	nz	Number of cells
[in]	h	Cell widths
[in,out]	x	Edge coordiantes starting at x(1)=0

Definition at line 141 of file MOM_remapping.F90.

  integer,               intent(in)    :: nz !< Number of cells
  real, dimension(nz),   intent(in)    :: h  !< Cell widths
  real, dimension(nz+1), intent(inout) :: x  !< Edge coordiantes starting at x(1)=0
  ! Local variables
  integer :: k
 
  x(1) = 0.0
  do k = 1,nz
    x(k+1) = x(k) + h(k)
  enddo
 

check_reconstructions_1d()

subroutine mom_remapping::check_reconstructions_1d ( integer, intent(in)  n0, real, dimension(n0), intent(in)  h0, real, dimension(n0), intent(in)  u0, integer, intent(in)  deg, logical, intent(in)  boundary_extrapolation, real, dimension(n0,deg+1), intent(out)  ppoly_r_coefs, real, dimension(n0,2), intent(out)  ppoly_r_E, real, dimension(n0,2), intent(out)  ppoly_r_S  )

private

Checks that edge values and reconstructions satisfy bounds.

Parameters

[in]	n0	Number of cells on source grid
[in]	h0	Cell widths on source grid
[in]	u0	Cell averages on source grid
[in]	deg	Degree of polynomial reconstruction
[in]	boundary_extrapolation	Extrapolate at boundaries if true
[out]	ppoly_r_coefs	Coefficients of polynomial
[out]	ppoly_r_e	Edge value of polynomial
[out]	ppoly_r_s	Edge slope of polynomial

Definition at line 444 of file MOM_remapping.F90.

  integer,                  intent(in)  :: n0 !< Number of cells on source grid
  real, dimension(n0),      intent(in)  :: h0 !< Cell widths on source grid
  real, dimension(n0),      intent(in)  :: u0 !< Cell averages on source grid
  integer,                  intent(in)  :: deg !< Degree of polynomial reconstruction
  logical,                  intent(in)  :: boundary_extrapolation !< Extrapolate at boundaries if true
  real, dimension(n0,deg+1),intent(out) :: ppoly_r_coefs !< Coefficients of polynomial
  real, dimension(n0,2),    intent(out) :: ppoly_r_E !< Edge value of polynomial
  real, dimension(n0,2),    intent(out) :: ppoly_r_S !< Edge slope of polynomial
  ! Local variables
  integer :: i0, n
  real :: u_l, u_c, u_r ! Cell averages
  real :: u_min, u_max
  logical :: problem_detected
 
  problem_detected = .false.
  do i0 = 1, n0
    u_l = u0(max(1,i0-1))
    u_c = u0(i0)
    u_r = u0(min(n0,i0+1))
    if (i0 > 1 .or. .not. boundary_extrapolation) then
      u_min = min(u_l, u_c)
      u_max = max(u_l, u_c)
      if (ppoly_r_e(i0,1) < u_min) then
        write(0,'(a,i4,5(x,a,1pe24.16))') 'Left edge undershoot at',i0,'u(i0-1)=',u_l,'u(i0)=',u_c, &
                                          'edge=',ppoly_r_e(i0,1),'err=',ppoly_r_e(i0,1)-u_min
        problem_detected = .true.
      endif
      if (ppoly_r_e(i0,1) > u_max) then
        write(0,'(a,i4,5(x,a,1pe24.16))') 'Left edge overshoot at',i0,'u(i0-1)=',u_l,'u(i0)=',u_c, &
                                          'edge=',ppoly_r_e(i0,1),'err=',ppoly_r_e(i0,1)-u_max
        problem_detected = .true.
      endif
    endif
    if (i0 < n0 .or. .not. boundary_extrapolation) then
      u_min = min(u_c, u_r)
      u_max = max(u_c, u_r)
      if (ppoly_r_e(i0,2) < u_min) then
        write(0,'(a,i4,5(x,a,1pe24.16))') 'Right edge undershoot at',i0,'u(i0)=',u_c,'u(i0+1)=',u_r, &
                                          'edge=',ppoly_r_e(i0,2),'err=',ppoly_r_e(i0,2)-u_min
        problem_detected = .true.
      endif
      if (ppoly_r_e(i0,2) > u_max) then
        write(0,'(a,i4,5(x,a,1pe24.16))') 'Right edge overshoot at',i0,'u(i0)=',u_c,'u(i0+1)=',u_r, &
                                          'edge=',ppoly_r_e(i0,2),'err=',ppoly_r_e(i0,2)-u_max
        problem_detected = .true.
      endif
    endif
    if (i0 > 1) then
      if ( (u_c-u_l)*(ppoly_r_e(i0,1)-ppoly_r_e(i0-1,2)) < 0.) then
        write(0,'(a,i4,5(x,a,1pe24.16))') 'Non-monotonic edges at',i0,'u(i0-1)=',u_l,'u(i0)=',u_c, &
                                          'right edge=',ppoly_r_e(i0-1,2),'left edge=',ppoly_r_e(i0,1)
        write(0,'(5(a,1pe24.16,x))') 'u(i0)-u(i0-1)',u_c-u_l,'edge diff=',ppoly_r_e(i0,1)-ppoly_r_e(i0-1,2)
        problem_detected = .true.
      endif
    endif
    if (problem_detected) then
      write(0,'(a,1p9e24.16)') 'Polynomial coeffs:',ppoly_r_coefs(i0,:)
      write(0,'(3(a,1pe24.16,x))') 'u_l=',u_l,'u_c=',u_c,'u_r=',u_r
      write(0,'(a4,10a24)') 'i0','h0(i0)','u0(i0)','left edge','right edge','Polynomial coefficients'
      do n = 1, n0
        write(0,'(i4,1p10e24.16)') n,h0(n),u0(n),ppoly_r_e(n,1),ppoly_r_e(n,2),ppoly_r_coefs(n,:)
      enddo
      call mom_error(fatal, 'MOM_remapping, check_reconstructions_1d: '// &
                   'Edge values or polynomial coefficients were inconsistent!')
    endif
  enddo
 

dumpgrid()

subroutine mom_remapping::dumpgrid ( integer, intent(in)  n, real, dimension(:), intent(in)  h, real, dimension(:), intent(in)  x, real, dimension(:), intent(in)  u  )

private

Convenience function for printing grid to screen.

Parameters

[in]	n	Number of cells
[in]	h	Cell thickness
[in]	x	Interface delta
[in]	u	Cell average values

Definition at line 1918 of file MOM_remapping.F90.

  integer, intent(in) :: n !< Number of cells
  real, dimension(:), intent(in) :: h !< Cell thickness
  real, dimension(:), intent(in) :: x !< Interface delta
  real, dimension(:), intent(in) :: u !< Cell average values
  integer :: i
  write(stdout,'("i=",20i10)') (i,i=1,n+1)
  write(stdout,'("x=",20es10.2)') (x(i),i=1,n+1)
  write(stdout,'("i=",5x,20i10)') (i,i=1,n)
  write(stdout,'("h=",5x,20es10.2)') (h(i),i=1,n)
  write(stdout,'("u=",5x,20es10.2)') (u(i),i=1,n)

dzfromh1h2()

subroutine, public mom_remapping::dzfromh1h2	(	integer, intent(in)	n1,
		real, dimension(:), intent(in)	h1,
		integer, intent(in)	n2,
		real, dimension(:), intent(in)	h2,
		real, dimension(:), intent(out)	dx
	)

Calculates the change in interface positions based on h1 and h2.

Parameters

[in]	n1	Number of cells on source grid
[in]	h1	Cell widths of source grid (size n1)
[in]	n2	Number of cells on target grid
[in]	h2	Cell widths of target grid (size n2)
[out]	dx	Change in interface position (size n2+1)

Definition at line 1525 of file MOM_remapping.F90.

  integer,            intent(in)  :: n1 !< Number of cells on source grid
  real, dimension(:), intent(in)  :: h1 !< Cell widths of source grid (size n1)
  integer,            intent(in)  :: n2 !< Number of cells on target grid
  real, dimension(:), intent(in)  :: h2 !< Cell widths of target grid (size n2)
  real, dimension(:), intent(out) :: dx !< Change in interface position (size n2+1)
  ! Local variables
  integer :: k
  real :: x1, x2
 
  x1 = 0.
  x2 = 0.
  dx(1) = 0.
  do k = 1, max(n1,n2)
    if (k <= n1) x1 = x1 + h1(k) ! Interface k+1, right of source cell k
    if (k <= n2) then
      x2 = x2 + h2(k) ! Interface k+1, right of target cell k
      dx(k+1) = x2 - x1 ! Change of interface k+1, target - source
    endif
  enddo
 

end_remapping()

subroutine, public mom_remapping::end_remapping

(

type(remapping_cs), intent(inout)

CS

)

Destrcutor for remapping control structure.

Parameters

[in,out]

cs

Remapping control structure

Definition at line 1606 of file MOM_remapping.F90.

  type(remapping_CS), intent(inout) :: CS !< Remapping control structure
 
  cs%degree = 0
 

extract_member_remapping_cs()

subroutine, public mom_remapping::extract_member_remapping_cs	(	type(remapping_cs), intent(in)	CS,
		integer, intent(out), optional	remapping_scheme,
		integer, intent(out), optional	degree,
		logical, intent(out), optional	boundary_extrapolation,
		logical, intent(out), optional	check_reconstruction,
		logical, intent(out), optional	check_remapping,
		logical, intent(out), optional	force_bounds_in_subcell
	)

Parameters

[in]	cs	Control structure for remapping module
[out]	remapping_scheme	Determines which reconstruction scheme to use
[out]	degree	Degree of polynomial reconstruction
[out]	boundary_extrapolation	If true, extrapolate boundaries
[out]	check_reconstruction	If true, reconstructions are checked for consistency.
[out]	check_remapping	If true, the result of remapping are checked for conservation and bounds.
[out]	force_bounds_in_subcell	If true, the intermediate values used in remapping are forced to be bounded.

Definition at line 120 of file MOM_remapping.F90.

  type(remapping_CS), intent(in) :: CS !< Control structure for remapping module
  integer, optional, intent(out) :: remapping_scheme        !< Determines which reconstruction scheme to use
  integer, optional, intent(out) :: degree                  !< Degree of polynomial reconstruction
  logical, optional, intent(out) :: boundary_extrapolation  !< If true, extrapolate boundaries
  logical, optional, intent(out) :: check_reconstruction    !< If true, reconstructions are checked for consistency.
  logical, optional, intent(out) :: check_remapping         !< If true, the result of remapping are checked
                                                            !!  for conservation and bounds.
  logical, optional, intent(out) :: force_bounds_in_subcell !< If true, the intermediate values used in
                                                            !! remapping are forced to be bounded.
 
  if (present(remapping_scheme)) remapping_scheme = cs%remapping_scheme
  if (present(degree)) degree = cs%degree
  if (present(boundary_extrapolation)) boundary_extrapolation = cs%boundary_extrapolation
  if (present(check_reconstruction)) check_reconstruction = cs%check_reconstruction
  if (present(check_remapping)) check_remapping = cs%check_remapping
  if (present(force_bounds_in_subcell)) force_bounds_in_subcell = cs%force_bounds_in_subcell
 

initialize_remapping()

subroutine, public mom_remapping::initialize_remapping	(	type(remapping_cs), intent(inout)	CS,
		character(len=*), intent(in)	remapping_scheme,
		logical, intent(in), optional	boundary_extrapolation,
		logical, intent(in), optional	check_reconstruction,
		logical, intent(in), optional	check_remapping,
		logical, intent(in), optional	force_bounds_in_subcell,
		logical, intent(in), optional	answers_2018
	)

Constructor for remapping control structure.

Parameters

[in,out]	cs	Remapping control structure
[in]	remapping_scheme	Remapping scheme to use
[in]	boundary_extrapolation	Indicate to extrapolate in boundary cells
[in]	check_reconstruction	Indicate to check reconstructions
[in]	check_remapping	Indicate to check results of remapping
[in]	force_bounds_in_subcell	Force subcells values to be bounded
[in]	answers_2018	If true use older, less acccurate expressions.

Definition at line 1549 of file MOM_remapping.F90.

  ! Arguments
  type(remapping_CS), intent(inout) :: CS !< Remapping control structure
  character(len=*),   intent(in)    :: remapping_scheme !< Remapping scheme to use
  logical, optional,  intent(in)    :: boundary_extrapolation !< Indicate to extrapolate in boundary cells
  logical, optional,  intent(in)    :: check_reconstruction !< Indicate to check reconstructions
  logical, optional,  intent(in)    :: check_remapping !< Indicate to check results of remapping
  logical, optional,  intent(in)    :: force_bounds_in_subcell !< Force subcells values to be bounded
  logical, optional,  intent(in)    :: answers_2018 !< If true use older, less acccurate expressions.
 
  ! Note that remapping_scheme is mandatory for initialize_remapping()
  call remapping_set_param(cs, remapping_scheme=remapping_scheme, boundary_extrapolation=boundary_extrapolation,  &
               check_reconstruction=check_reconstruction, check_remapping=check_remapping, &
               force_bounds_in_subcell=force_bounds_in_subcell, answers_2018=answers_2018)
 

integraterecononinterval()

subroutine mom_remapping::integraterecononinterval ( integer, intent(in)  n0, real, dimension(:), intent(in)  h0, real, dimension(:), intent(in)  u0, real, dimension(:,:), intent(in)  ppoly0_E, real, dimension(:,:), intent(in)  ppoly0_coefs, integer, intent(in)  method, real, intent(in)  xL, real, intent(in)  xR, real, intent(in)  hC, real, intent(out)  uAve, integer, intent(inout)  jStart, real, intent(inout)  xStart, real, intent(in), optional  h_neglect  )

private

Integrate the reconstructed column profile over a single cell.

Parameters

[in]	n0	Number of cells in source grid
[in]	h0	Source grid sizes (size n0)
[in]	u0	Source cell averages
[in]	ppoly0_e	Edge value of polynomial
[in]	ppoly0_coefs	Coefficients of polynomial
[in]	method	Remapping scheme to use
[in]	xl	Left edges of target cell
[in]	xr	Right edges of target cell
[in]	hc	Cell width hC = xR - xL
[out]	uave	Average value on target cell
[in,out]	jstart	The index of the cell to start searching from On exit, contains index of last cell used
[in,out]	xstart	The left edge position of cell jStart On first entry should be 0.
[in]	h_neglect	A negligibly small width for the purpose of cell reconstructions in the same units as h.

Definition at line 1225 of file MOM_remapping.F90.

  integer,              intent(in)    :: n0     !< Number of cells in source grid
  real, dimension(:),   intent(in)    :: h0     !< Source grid sizes (size n0)
  real, dimension(:),   intent(in)    :: u0     !< Source cell averages
  real, dimension(:,:), intent(in)    :: ppoly0_E !< Edge value of polynomial
  real, dimension(:,:), intent(in)    :: ppoly0_coefs !< Coefficients of polynomial
  integer,              intent(in)    :: method !< Remapping scheme to use
  real,                 intent(in)    :: xL     !< Left edges of target cell
  real,                 intent(in)    :: xR     !< Right edges of target cell
  real,                 intent(in)    :: hC     !< Cell width hC = xR - xL
  real,                 intent(out)   :: uAve   !< Average value on target cell
  integer,              intent(inout) :: jStart !< The index of the cell to start searching from
                                   !< On exit, contains index of last cell used
  real,                 intent(inout) :: xStart !< The left edge position of cell jStart
                                   !< On first entry should be 0.
  real,       optional, intent(in)    :: h_neglect !< A negligibly small width for the
                                          !! purpose of cell reconstructions
                                          !! in the same units as h.
  ! Local variables
  integer :: j, k
  integer :: jL, jR       ! indexes of source cells containing target
                          ! cell edges
  real    :: q            ! complete integration
  real    :: xi0, xi1     ! interval of integration (local -- normalized
                          ! -- coordinates)
  real    :: x0jLl, x0jLr ! Left/right position of cell jL
  real    :: x0jRl, x0jRr ! Left/right position of cell jR
  real    :: hAct         ! The distance actually used in the integration
                          ! (notionally xR - xL) which differs due to roundoff.
  real    :: x0_2, x1_2, x02px12, x0px1 ! Used in evaluation of integrated polynomials
  real     :: hNeglect    ! A negligible thicness in the same units as h.
  real, parameter :: r_3 = 1.0/3.0 ! Used in evaluation of integrated polynomials
 
  hneglect = hneglect_dflt ; if (present(h_neglect)) hneglect = h_neglect
 
  q = -1.e30
  x0jll = -1.e30
  x0jrl = -1.e30
 
  ! Find the left most cell in source grid spanned by the target cell
  jl = -1
  x0jlr = xstart
  do j = jstart, n0
    x0jll = x0jlr
    x0jlr = x0jll + h0(j)
    ! Left edge is found in cell j
    if ( ( xl >= x0jll ) .AND. ( xl <= x0jlr ) ) then
      jl = j
      exit ! once target grid cell is found, exit loop
    endif
  enddo
  jstart = jl
  xstart = x0jll
 
! ! HACK to handle round-off problems. Need only at  j=n0.
! ! This moves the effective cell boundary outwards a smidgen.
! if (xL>x0jLr) x0jLr = xL
 
  ! If, at this point, jL is equal to -1, it means the vanished
  ! cell lies outside the source grid. In other words, it means that
  ! the source and target grids do not cover the same physical domain
  ! and there is something very wrong !
  if ( jl == -1 ) call mom_error(fatal, &
          'MOM_remapping, integrateReconOnInterval: '//&
          'The location of the left-most cell could not be found')
 
 
  ! ============================================================
  ! Check whether target cell is vanished. If it is, the cell
  ! average is simply the interpolated value at the location
  ! of the vanished cell. If it isn't, we need to integrate the
  ! quantity within the cell and divide by the cell width to
  ! determine the cell average.
  ! ============================================================
  ! 1. Cell is vanished
 !if ( abs(xR - xL) <= epsilon(xR)*max(abs(xR),abs(xL)) ) then
  if ( abs(xr - xl) == 0.0 ) then
 
    ! We check whether the source cell (i.e. the cell in which the
    ! vanished target cell lies) is vanished. If it is, the interpolated
    ! value is set to be mean of the edge values (which should be the same).
    ! If it isn't, we simply interpolate.
    if ( h0(jl) == 0.0 ) then
      uave = 0.5 * ( ppoly0_e(jl,1) + ppoly0_e(jl,2) )
    else
      ! WHY IS THIS NOT WRITTEN AS xi0 = ( xL - x0jLl ) / h0(jL) ---AJA
      xi0 = xl / ( h0(jl) + hneglect ) - x0jll / ( h0(jl) + hneglect )
 
      select case ( method )
        case ( integration_pcm )
          uave =         ppoly0_coefs(jl,1)
        case ( integration_plm )
          uave =         ppoly0_coefs(jl,1)   &
               + xi0 *   ppoly0_coefs(jl,2)
        case ( integration_ppm )
          uave =         ppoly0_coefs(jl,1)   &
               + xi0 * ( ppoly0_coefs(jl,2)   &
               + xi0 *   ppoly0_coefs(jl,3) )
        case ( integration_pqm )
          uave =         ppoly0_coefs(jl,1)   &
               + xi0 * ( ppoly0_coefs(jl,2)   &
               + xi0 * ( ppoly0_coefs(jl,3)   &
               + xi0 * ( ppoly0_coefs(jl,4)   &
               + xi0 *   ppoly0_coefs(jl,5) ) ) )
        case default
          call mom_error( fatal,'The selected integration method is invalid' )
      end select
 
    endif ! end checking whether source cell is vanished
 
  ! 2. Cell is not vanished
  else
 
    ! Find the right most cell in source grid spanned by the target cell
    jr = -1
    x0jrr = xstart
    do j = jstart,n0
      x0jrl = x0jrr
      x0jrr = x0jrl + h0(j)
      ! Right edge is found in cell j
      if ( ( xr >= x0jrl ) .AND. ( xr <= x0jrr ) ) then
        jr = j
        exit  ! once target grid cell is found, exit loop
      endif
    enddo ! end loop on source grid cells
 
    ! If xR>x0jRr then the previous loop reached j=n0 and the target
    ! position, xR, was beyond the right edge of the source grid (h0).
    ! This can happen due to roundoff, in which case we set jR=n0.
    if (xr>x0jrr) jr = n0
 
    ! To integrate, two cases must be considered: (1) the target cell is
    ! entirely contained within a cell of the source grid and (2) the target
    ! cell spans at least two cells of the source grid.
 
    if ( jl == jr ) then
      ! The target cell is entirely contained within a cell of the source
      ! grid. This situation is represented by the following schematic, where
      ! the cell in which xL and xR are located has index jL=jR :
      !
      ! ----|-----o--------o----------|-------------
      !           xL       xR
      !
      ! Determine normalized coordinates
#ifdef __USE_ROUNDOFF_SAFE_ADJUSTMENTS__
      xi0 = max( 0., min( 1., ( xl - x0jll ) / ( h0(jl) + hneglect ) ) )
      xi1 = max( 0., min( 1., ( xr - x0jll ) / ( h0(jl) + hneglect ) ) )
#else
      xi0 = xl / h0(jl) - x0jll / ( h0(jl) + hneglect )
      xi1 = xr / h0(jl) - x0jll / ( h0(jl) + hneglect )
#endif
 
      hact = h0(jl) * ( xi1 - xi0 )
 
      ! Depending on which polynomial is used, integrate quantity
      ! between xi0 and xi1. Integration is carried out in normalized
      ! coordinates, hence: \int_xL^xR p(x) dx = h \int_xi0^xi1 p(xi) dxi
      select case ( method )
        case ( integration_pcm )
          q = ( xr - xl ) * ppoly0_coefs(jl,1)
        case ( integration_plm )
          q = ( xr - xl ) * (                            &
              ppoly0_coefs(jl,1)                         &
            + ppoly0_coefs(jl,2) * 0.5 * ( xi1 + xi0 ) )
        case ( integration_ppm )
          q = ( xr - xl ) * (                            &
                ppoly0_coefs(jl,1)                       &
            + ( ppoly0_coefs(jl,2) * 0.5 * ( xi1 + xi0 ) &
            +   ppoly0_coefs(jl,3) * r_3 * ( ( xi1*xi1 + xi0*xi0 ) + xi0*xi1 ) ) )
        case ( integration_pqm )
          x0_2 = xi0*xi0
          x1_2 = xi1*xi1
          x02px12 = x0_2 + x1_2
          x0px1 = xi1 + xi0
          q = ( xr - xl ) * (                                    &
                ppoly0_coefs(jl,1)                               &
            + ( ppoly0_coefs(jl,2) * 0.5 * ( xi1 + xi0 )         &
            + ( ppoly0_coefs(jl,3) * r_3 * ( x02px12 + xi0*xi1 ) &
            +   ppoly0_coefs(jl,4) * 0.25* ( x02px12 * x0px1 )   &
            +   ppoly0_coefs(jl,5) * 0.2 * ( ( xi1*x1_2 + xi0*x0_2 ) * x0px1 + x0_2*x1_2 ) ) ) )
        case default
          call mom_error( fatal,'The selected integration method is invalid' )
      end select
 
    else
    ! The target cell spans at least two cells of the source grid.
    ! This situation is represented by the following schematic, where
    ! the cells in which xL and xR are located have indexes jL and jR,
    ! respectively :
    !
    ! ----|-----o---|--- ... --|---o----------|-------------
    !           xL                 xR
    !
    ! We first integrate from xL up to the right boundary of cell jL, then
    ! add the integrated amounts of cells located between jL and jR and then
    ! integrate from the left boundary of cell jR up to xR
 
      q = 0.0
 
      ! Integrate from xL up to right boundary of cell jL
#ifdef __USE_ROUNDOFF_SAFE_ADJUSTMENTS__
      xi0 = max( 0., min( 1., ( xl - x0jll ) / ( h0(jl) + hneglect ) ) )
#else
      xi0 = (xl - x0jll) / ( h0(jl) + hneglect )
#endif
      xi1 = 1.0
 
      hact = h0(jl) * ( xi1 - xi0 )
 
      select case ( method )
        case ( integration_pcm )
          q = q + ( x0jlr - xl ) * ppoly0_coefs(jl,1)
        case ( integration_plm )
          q = q + ( x0jlr - xl ) * (                     &
              ppoly0_coefs(jl,1)                         &
            + ppoly0_coefs(jl,2) * 0.5 * ( xi1 + xi0 ) )
        case ( integration_ppm )
          q = q + ( x0jlr - xl ) * (                     &
                ppoly0_coefs(jl,1)                       &
            + ( ppoly0_coefs(jl,2) * 0.5 * ( xi1 + xi0 ) &
            +   ppoly0_coefs(jl,3) * r_3 * ( ( xi1*xi1 + xi0*xi0 ) + xi0*xi1 ) ) )
        case ( integration_pqm )
          x0_2 = xi0*xi0
          x1_2 = xi1*xi1
          x02px12 = x0_2 + x1_2
          x0px1 = xi1 + xi0
          q = q + ( x0jlr - xl ) * (                             &
                ppoly0_coefs(jl,1)                               &
            + ( ppoly0_coefs(jl,2) * 0.5 * ( xi1 + xi0 )         &
            + ( ppoly0_coefs(jl,3) * r_3 * ( x02px12 + xi0*xi1 ) &
            +   ppoly0_coefs(jl,4) * 0.25* ( x02px12 * x0px1 )   &
            +   ppoly0_coefs(jl,5) * 0.2 * ( ( xi1*x1_2 + xi0*x0_2 ) * x0px1 + x0_2*x1_2 ) ) ) )
        case default
          call mom_error( fatal, 'The selected integration method is invalid' )
      end select
 
      ! Integrate contents within cells strictly comprised between jL and jR
      if ( jr > (jl+1) ) then
        do k = jl+1,jr-1
          q = q + h0(k) * u0(k)
          hact = hact + h0(k)
        enddo
      endif
 
      ! Integrate from left boundary of cell jR up to xR
      xi0 = 0.0
#ifdef __USE_ROUNDOFF_SAFE_ADJUSTMENTS__
      xi1 = max( 0., min( 1., ( xr - x0jrl ) / ( h0(jr) + hneglect ) ) )
#else
      xi1 = (xr - x0jrl) / ( h0(jr) + hneglect )
#endif
 
      hact = hact + h0(jr) * ( xi1 - xi0 )
 
      select case ( method )
        case ( integration_pcm )
          q = q + ( xr - x0jrl ) * ppoly0_coefs(jr,1)
        case ( integration_plm )
          q = q + ( xr - x0jrl ) * (                     &
              ppoly0_coefs(jr,1)                         &
            + ppoly0_coefs(jr,2) * 0.5 * ( xi1 + xi0 ) )
        case ( integration_ppm )
          q = q + ( xr - x0jrl ) * (                     &
                ppoly0_coefs(jr,1)                       &
            + ( ppoly0_coefs(jr,2) * 0.5 * ( xi1 + xi0 ) &
            +   ppoly0_coefs(jr,3) * r_3 * ( ( xi1*xi1 + xi0*xi0 ) + xi0*xi1 ) ) )
        case ( integration_pqm )
          x0_2 = xi0*xi0
          x1_2 = xi1*xi1
          x02px12 = x0_2 + x1_2
          x0px1 = xi1 + xi0
          q = q + ( xr - x0jrl ) * (                             &
                ppoly0_coefs(jr,1)                               &
            + ( ppoly0_coefs(jr,2) * 0.5 * ( xi1 + xi0 )         &
            + ( ppoly0_coefs(jr,3) * r_3 * ( x02px12 + xi0*xi1 ) &
            +   ppoly0_coefs(jr,4) * 0.25* ( x02px12 * x0px1 )   &
            +   ppoly0_coefs(jr,5) * 0.2 * ( ( xi1*x1_2 + xi0*x0_2 ) * x0px1 + x0_2*x1_2 ) ) ) )
        case default
          call mom_error( fatal,'The selected integration method is invalid' )
      end select
 
    endif ! end integration for non-vanished cells
 
    ! The cell average is the integrated value divided by the cell width
#ifdef __USE_ROUNDOFF_SAFE_ADJUSTMENTS__
if (hact==0.) then
    uave = ppoly0_coefs(jl,1)
else
    uave = q / hact
endif
#else
    uave = q / hc
#endif
 
  endif ! endif clause to check if cell is vanished
 

ispossumerrsignificant()

logical function mom_remapping::ispossumerrsignificant ( integer, intent(in)  n1, real, intent(in)  sum1, integer, intent(in)  n2, real, intent(in)  sum2  )

private

Compare two summation estimates of positive data and judge if due to more than round-off. When two sums are calculated from different vectors that should add up to the same value, the results can differ by round off. The round off error can be bounded to be proportional to the number of operations. This function returns true if the difference between sum1 and sum2 is larger than than the estimated round off bound.

Note

This estimate/function is only valid for summation of positive data.

Parameters

[in]	n1	Number of values in sum1
[in]	n2	Number of values in sum2
[in]	sum1	Sum of n1 values
[in]	sum2	Sum of n2 values

Returns

True if difference in sums is large

Definition at line 163 of file MOM_remapping.F90.

  integer, intent(in) :: n1   !< Number of values in sum1
  integer, intent(in) :: n2   !< Number of values in sum2
  real,    intent(in) :: sum1 !< Sum of n1 values
  real,    intent(in) :: sum2 !< Sum of n2 values
  logical             :: isPosSumErrSignificant !< True if difference in sums is large
  ! Local variables
  real :: sumErr, allowedErr, eps
 
  if (sum1<0.) call mom_error(fatal,'isPosSumErrSignificant: sum1<0 is not allowed!')
  if (sum2<0.) call mom_error(fatal,'isPosSumErrSignificant: sum2<0 is not allowed!')
  sumerr = abs(sum1-sum2)
  eps = epsilon(sum1)
  allowederr = eps*0.5*(real(n1-1)*sum1+real(n2-1)*sum2)
  if (sumerr>allowederr) then
    write(0,*) 'isPosSumErrSignificant: sum1,sum2=',sum1,sum2
    write(0,*) 'isPosSumErrSignificant: eps=',eps
    write(0,*) 'isPosSumErrSignificant: err,n*eps=',sumerr,allowederr
    write(0,*) 'isPosSumErrSignificant: err/eps,n1,n2,n1+n2=',sumerr/eps,n1,n2,n1+n2
    ispossumerrsignificant = .true.
  else
    ispossumerrsignificant = .false.
  endif

measure_input_bounds()

subroutine mom_remapping::measure_input_bounds ( integer, intent(in)  n0, real, dimension(n0), intent(in)  h0, real, dimension(n0), intent(in)  u0, real, dimension(n0,2), intent(in)  edge_values, real, intent(out)  h0tot, real, intent(out)  h0err, real, intent(out)  u0tot, real, intent(out)  u0err, real, intent(out)  u0min, real, intent(out)  u0max  )

private

Measure totals and bounds on source grid.

Parameters

[in]	n0	Number of cells on source grid
[in]	h0	Cell widths on source grid
[in]	u0	Cell averages on source grid
[in]	edge_values	Cell edge values on source grid
[out]	h0tot	Sum of cell widths
[out]	h0err	Magnitude of round-off error in h0tot
[out]	u0tot	Sum of cell widths times values
[out]	u0err	Magnitude of round-off error in u0tot
[out]	u0min	Minimum value in reconstructions of u0
[out]	u0max	Maximum value in reconstructions of u0

Definition at line 1030 of file MOM_remapping.F90.

  integer,               intent(in)  :: n0 !< Number of cells on source grid
  real, dimension(n0),   intent(in)  :: h0 !< Cell widths on source grid
  real, dimension(n0),   intent(in)  :: u0 !< Cell averages on source grid
  real, dimension(n0,2), intent(in)  :: edge_values !< Cell edge values on source grid
  real,                  intent(out) :: h0tot !< Sum of cell widths
  real,                  intent(out) :: h0err !< Magnitude of round-off error in h0tot
  real,                  intent(out) :: u0tot !< Sum of cell widths times values
  real,                  intent(out) :: u0err !< Magnitude of round-off error in u0tot
  real,                  intent(out) :: u0min !< Minimum value in reconstructions of u0
  real,                  intent(out) :: u0max !< Maximum value in reconstructions of u0
  ! Local variables
  integer :: k
  real :: eps
 
  eps = epsilon(h0(1))
  h0tot = h0(1)
  h0err = 0.
  u0tot = h0(1) * u0(1)
  u0err = 0.
  u0min = min( edge_values(1,1), edge_values(1,2) )
  u0max = max( edge_values(1,1), edge_values(1,2) )
  do k = 2, n0
    h0tot = h0tot + h0(k)
    h0err = h0err + eps * max(h0tot, h0(k))
    u0tot = u0tot + h0(k) * u0(k)
    u0err = u0err + eps * max(abs(u0tot), abs(h0(k) * u0(k)))
    u0min = min( u0min, edge_values(k,1), edge_values(k,2) )
    u0max = max( u0max, edge_values(k,1), edge_values(k,2) )
  enddo
 

measure_output_bounds()

subroutine mom_remapping::measure_output_bounds ( integer, intent(in)  n1, real, dimension(n1), intent(in)  h1, real, dimension(n1), intent(in)  u1, real, intent(out)  h1tot, real, intent(out)  h1err, real, intent(out)  u1tot, real, intent(out)  u1err, real, intent(out)  u1min, real, intent(out)  u1max  )

private

Measure totals and bounds on destination grid.

Parameters

[in]	n1	Number of cells on destination grid
[in]	h1	Cell widths on destination grid
[in]	u1	Cell averages on destination grid
[out]	h1tot	Sum of cell widths
[out]	h1err	Magnitude of round-off error in h1tot
[out]	u1tot	Sum of cell widths times values
[out]	u1err	Magnitude of round-off error in u1tot
[out]	u1min	Minimum value in reconstructions of u1
[out]	u1max	Maximum value in reconstructions of u1

Definition at line 1064 of file MOM_remapping.F90.

  integer,               intent(in)  :: n1 !< Number of cells on destination grid
  real, dimension(n1),   intent(in)  :: h1 !< Cell widths on destination grid
  real, dimension(n1),   intent(in)  :: u1 !< Cell averages on destination grid
  real,                  intent(out) :: h1tot !< Sum of cell widths
  real,                  intent(out) :: h1err !< Magnitude of round-off error in h1tot
  real,                  intent(out) :: u1tot !< Sum of cell widths times values
  real,                  intent(out) :: u1err !< Magnitude of round-off error in u1tot
  real,                  intent(out) :: u1min !< Minimum value in reconstructions of u1
  real,                  intent(out) :: u1max !< Maximum value in reconstructions of u1
  ! Local variables
  integer :: k
  real :: eps
 
  eps = epsilon(h1(1))
  h1tot = h1(1)
  h1err = 0.
  u1tot = h1(1) * u1(1)
  u1err = 0.
  u1min = u1(1)
  u1max = u1(1)
  do k = 2, n1
    h1tot = h1tot + h1(k)
    h1err = h1err + eps * max(h1tot, h1(k))
    u1tot = u1tot + h1(k) * u1(k)
    u1err = u1err + eps * max(abs(u1tot), abs(h1(k) * u1(k)))
    u1min = min(u1min, u1(k))
    u1max = max(u1max, u1(k))
  enddo
 

remap_via_sub_cells()

subroutine mom_remapping::remap_via_sub_cells ( integer, intent(in)  n0, real, dimension(n0), intent(in)  h0, real, dimension(n0), intent(in)  u0, real, dimension(n0,2), intent(in)  ppoly0_E, real, dimension(:,:), intent(in)  ppoly0_coefs, integer, intent(in)  n1, real, dimension(n1), intent(in)  h1, integer, intent(in)  method, logical, intent(in)  force_bounds_in_subcell, real, dimension(n1), intent(out)  u1, real, intent(out)  uh_err, real, dimension(n0+n1+1), intent(out), optional  ah_sub, integer, dimension(n0+n1+1), intent(out), optional  aisub_src, integer, dimension(n0), intent(out), optional  aiss, integer, dimension(n0), intent(out), optional  aise  )

private

Remaps column of n0 values u0 on grid h0 to grid h1 with n1 cells by calculating the n0+n1+1 sub-integrals of the intersection of h0 and h1, and the summing the appropriate integrals into the h1*u1 values. h0 and h1 must have the same units.

Parameters

[in]	n0	Number of cells in source grid
[in]	h0	Source grid widths (size n0)
[in]	u0	Source cell averages (size n0)
[in]	ppoly0_e	Edge value of polynomial
[in]	ppoly0_coefs	Coefficients of polynomial
[in]	n1	Number of cells in target grid
[in]	h1	Target grid widths (size n1)
[in]	method	Remapping scheme to use
[in]	force_bounds_in_subcell	Force sub-cell values to be bounded
[out]	u1	Target cell averages (size n1)
[out]	uh_err	Estimate of bound on error in sum of u*h
[out]	ah_sub	h_sub
[out]	aisub_src	i_sub_src
[out]	aiss	isrc_start
[out]	aise	isrc_ens

Definition at line 518 of file MOM_remapping.F90.

  integer,           intent(in)    :: n0     !< Number of cells in source grid
  real,              intent(in)    :: h0(n0)  !< Source grid widths (size n0)
  real,              intent(in)    :: u0(n0)  !< Source cell averages (size n0)
  real,              intent(in)    :: ppoly0_E(n0,2)            !< Edge value of polynomial
  real,              intent(in)    :: ppoly0_coefs(:,:) !< Coefficients of polynomial
  integer,           intent(in)    :: n1     !< Number of cells in target grid
  real,              intent(in)    :: h1(n1)  !< Target grid widths (size n1)
  integer,           intent(in)    :: method !< Remapping scheme to use
  logical,           intent(in)    :: force_bounds_in_subcell !< Force sub-cell values to be bounded
  real,              intent(out)   :: u1(n1)  !< Target cell averages (size n1)
  real,              intent(out)   :: uh_err !< Estimate of bound on error in sum of u*h
  real, optional,    intent(out)   :: ah_sub(n0+n1+1) !< h_sub
  integer, optional, intent(out)   :: aisub_src(n0+n1+1) !< i_sub_src
  integer, optional, intent(out)   :: aiss(n0) !< isrc_start
  integer, optional, intent(out)   :: aise(n0) !< isrc_ens
  ! 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
  real, dimension(n0+n1+1) :: h_sub ! Width of each each sub-cell
  real, dimension(n0+n1+1) :: uh_sub ! Integral of u*h over each sub-cell
  real, dimension(n0+n1+1) :: u_sub ! Average of u over each sub-cell
  integer, dimension(n0+n1+1) :: isub_src ! Index of source cell for each sub-cell
  integer, dimension(n0) :: isrc_start ! Index of first sub-cell within each source cell
  integer, dimension(n0) :: isrc_end ! Index of last sub-cell within each source cell
  integer, dimension(n0) :: isrc_max ! Index of thickest sub-cell within each source cell
  real, dimension(n0) :: h0_eff ! Effective thickness of source cells
  real, dimension(n0) :: u0_min ! Minimum value of reconstructions in source cell
  real, dimension(n0) :: u0_max ! Minimum value of reconstructions in source cell
  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)
  real :: h0_supply, h1_supply ! The amount of width available for constructing sub-cells
  real :: dh ! The width of the sub-cell
  real :: duh ! The total amount of accumulated stuff (u*h)
  real :: dh0_eff ! Running sum of source cell thickness
  ! For error checking/debugging
  logical, parameter :: force_bounds_in_target = .true. ! To fix round-off issues
  logical, parameter :: adjust_thickest_subcell = .true. ! To fix round-off conservation issues
  logical, parameter :: debug_bounds = .false. ! For debugging overshoots etc.
  integer :: k, i0_last_thick_cell
  real :: h0tot, h0err, h1tot, h1err, h2tot, h2err, u02_err
  real :: u0tot, u0err, u0min, u0max, u1tot, u1err, u1min, u1max, u2tot, u2err, u2min, u2max, u_orig
  logical :: src_has_volume !< True if h0 has not been consumed
  logical :: tgt_has_volume !< True if h1 has not been consumed
 
  if (old_algorithm) isrc_max(:)=1
 
  i0_last_thick_cell = 0
  do i0 = 1, n0
    u0_min(i0) = min(ppoly0_e(i0,1), ppoly0_e(i0,2))
    u0_max(i0) = max(ppoly0_e(i0,1), ppoly0_e(i0,2))
    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 (old_algorithm) then
        if (i0 < i0_last_thick_cell) then
          i0 = i0 + 1
          h0_supply = h0(i0)
          dh0_eff = 0.
          do while (h0_supply==0. .and. i0<i0_last_thick_cell)
            ! This loop skips over vanished source cells
            i0 = i0 + 1
            h0_supply = h0(i0)
          enddo
        else
          h0_supply = 1.e30
        endif
      else
        if (i0 < n0) then
          i0 = i0 + 1
          h0_supply = h0(i0)
          dh0_eff = 0.
        else
          h0_supply = 0.
          src_has_volume = .false.
        endif
      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
        if (old_algorithm) then
          h1_supply = 1.e30
        else
          h1_supply = 0.
          tgt_has_volume = .false.
        endif
      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
 
  ! Loop over each sub-cell to calculate average/integral values within each sub-cell.
  xa = 0.
  dh0_eff = 0.
  uh_sub(1) = 0.
  u_sub(1) = ppoly0_e(1,1)
  u02_err = 0.
  do i_sub = 2, n0+n1
 
    ! Sub-cell thickness from loop above
    dh = h_sub(i_sub)
 
    ! Source cell
    i0 = isub_src(i_sub)
 
    ! Evaluate average and integral for sub-cell i_sub.
    ! Integral is over distance dh but expressed in terms of non-dimensional
    ! positions with source cell from xa to xb  (0 <= xa <= xb <= 1).
    dh0_eff = dh0_eff + dh ! Cumulative thickness within the source cell
    if (h0_eff(i0)>0.) then
      xb = dh0_eff / h0_eff(i0) ! This expression yields xa <= xb <= 1.0
      xb = min(1., xb) ! This is only needed when the total target column is wider than the source column
      u_sub(i_sub) = average_value_ppoly( n0, u0, ppoly0_e, ppoly0_coefs, method, i0, xa, xb)
    else ! Vanished cell
      xb = 1.
      u_sub(i_sub) = u0(i0)
    endif
    if (debug_bounds) then
      if (method<5 .and.(u_sub(i_sub)<u0_min(i0) .or. u_sub(i_sub)>u0_max(i0))) then
        write(0,*) 'Sub cell average is out of bounds',i_sub,'method=',method
        write(0,*) 'xa,xb: ',xa,xb
        write(0,*) 'Edge values: ',ppoly0_e(i0,:),'mean',u0(i0)
        write(0,*) 'a_c: ',(u0(i0)-ppoly0_e(i0,1))+(u0(i0)-ppoly0_e(i0,2))
        write(0,*) 'Polynomial coeffs: ',ppoly0_coefs(i0,:)
        write(0,*) 'Bounds min=',u0_min(i0),'max=',u0_max(i0)
        write(0,*) 'Average: ',u_sub(i_sub),'rel to min=',u_sub(i_sub)-u0_min(i0),'rel to max=',u_sub(i_sub)-u0_max(i0)
        call mom_error( fatal, 'MOM_remapping, remap_via_sub_cells: '//&
             'Sub-cell average is out of bounds!' )
      endif
    endif
    if (force_bounds_in_subcell) then
      ! These next two lines should not be needed but when using PQM we found roundoff
      ! can lead to overshoots. These lines sweep issues under the rug which need to be
      ! properly .. later. -AJA
      u_orig = u_sub(i_sub)
      u_sub(i_sub) = max( u_sub(i_sub), u0_min(i0) )
      u_sub(i_sub) = min( u_sub(i_sub), u0_max(i0) )
      u02_err = u02_err + dh*abs( u_sub(i_sub) - u_orig )
    endif
    uh_sub(i_sub) = dh * u_sub(i_sub)
 
    if (isub_src(i_sub+1) /= i0) then
      ! If the next sub-cell is in a different source cell, reset the position counters
      dh0_eff = 0.
      xa = 0.
    else
      xa = xb ! Next integral will start at end of last
    endif
 
  enddo
  u_sub(n0+n1+1) = ppoly0_e(n0,2)                   ! This value is only needed when total target column
  uh_sub(n0+n1+1) = ppoly0_e(n0,2) * h_sub(n0+n1+1) ! is wider than the source column
 
  if (adjust_thickest_subcell) then
    ! 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).
    do i0 = 1, i0_last_thick_cell
      i_max = isrc_max(i0)
      dh_max = h_sub(i_max)
      if (dh_max > 0.) then
        ! duh will be the sum of sub-cell integrals within the source cell except for the thickest sub-cell.
        duh = 0.
        do i_sub = isrc_start(i0), isrc_end(i0)
          if (i_sub /= i_max) duh = duh + uh_sub(i_sub)
        enddo
        uh_sub(i_max) = u0(i0)*h0(i0) - duh
        u02_err = u02_err + max( abs(uh_sub(i_max)), abs(u0(i0)*h0(i0)), abs(duh) )
      endif
    enddo
  endif
 
  ! Loop over each target cell summing the integrals from sub-cells within the target cell.
  uh_err = 0.
  do i1 = 1, n1
    if (h1(i1) > 0.) then
      duh = 0. ; dh = 0.
      i_sub = itgt_start(i1)
      if (force_bounds_in_target) then
        u1min = u_sub(i_sub)
        u1max = u_sub(i_sub)
      endif
      do i_sub = itgt_start(i1), itgt_end(i1)
        if (force_bounds_in_target) then
          u1min = min(u1min, u_sub(i_sub))
          u1max = max(u1max, u_sub(i_sub))
        endif
        dh = dh + h_sub(i_sub)
        duh = duh + uh_sub(i_sub)
        ! This accumulates the contribution to the error bound for the sum of u*h
        uh_err = uh_err + max(abs(duh),abs(uh_sub(i_sub)))*epsilon(duh)
      enddo
      u1(i1) = duh / dh
      ! This is the contribution from the division to the error bound for the sum of u*h
      uh_err = uh_err + abs(duh)*epsilon(duh)
      if (force_bounds_in_target) then
        u_orig = u1(i1)
        u1(i1) = max(u1min, min(u1max, u1(i1)))
        ! Adjusting to be bounded contributes to the error for the sum of u*h
        uh_err = uh_err + dh*abs( u1(i1)-u_orig )
      endif
    else
      u1(i1) = u_sub(itgt_start(i1))
    endif
  enddo
 
  ! Check errors and bounds
  if (debug_bounds) then
    call measure_input_bounds( n0, h0, u0, ppoly0_e, h0tot, h0err, u0tot, u0err, u0min, u0max )
    call measure_output_bounds( n1, h1, u1, h1tot, h1err, u1tot, u1err, u1min, u1max )
    call measure_output_bounds( n0+n1+1, h_sub, u_sub, h2tot, h2err, u2tot, u2err, u2min, u2max )
    if (method<5) then ! We except PQM until we've debugged it
    if (     (abs(u1tot-u0tot)>(u0err+u1err)+uh_err+u02_err .and. abs(h1tot-h0tot)<h0err+h1err) &
        .or. (abs(u2tot-u0tot)>u0err+u2err+u02_err .and. abs(h2tot-h0tot)<h0err+h2err) &
        .or. (u1min<u0min .or. u1max>u0max) ) then
      write(0,*) 'method = ',method
      write(0,*) 'Source to sub-cells:'
      write(0,*) 'H: h0tot=',h0tot,'h2tot=',h2tot,'dh=',h2tot-h0tot,'h0err=',h0err,'h2err=',h2err
      if (abs(h2tot-h0tot)>h0err+h2err) &
        write(0,*) 'H non-conservation difference=',h2tot-h0tot,'allowed err=',h0err+h2err,' <-----!'
      write(0,*) 'UH: u0tot=',u0tot,'u2tot=',u2tot,'duh=',u2tot-u0tot,'u0err=',u0err,'u2err=',u2err,&
                 'adjustment err=',u02_err
      if (abs(u2tot-u0tot)>u0err+u2err) &
        write(0,*) 'U non-conservation difference=',u2tot-u0tot,'allowed err=',u0err+u2err,' <-----!'
      write(0,*) 'Sub-cells to target:'
      write(0,*) 'H: h2tot=',h2tot,'h1tot=',h1tot,'dh=',h1tot-h2tot,'h2err=',h2err,'h1err=',h1err
      if (abs(h1tot-h2tot)>h2err+h1err) &
        write(0,*) 'H non-conservation difference=',h1tot-h2tot,'allowed err=',h2err+h1err,' <-----!'
      write(0,*) 'UH: u2tot=',u2tot,'u1tot=',u1tot,'duh=',u1tot-u2tot,'u2err=',u2err,'u1err=',u1err,'uh_err=',uh_err
      if (abs(u1tot-u2tot)>u2err+u1err) &
        write(0,*) 'U non-conservation difference=',u1tot-u2tot,'allowed err=',u2err+u1err,' <-----!'
      write(0,*) 'Source to target:'
      write(0,*) 'H: h0tot=',h0tot,'h1tot=',h1tot,'dh=',h1tot-h0tot,'h0err=',h0err,'h1err=',h1err
      if (abs(h1tot-h0tot)>h0err+h1err) &
        write(0,*) 'H non-conservation difference=',h1tot-h0tot,'allowed err=',h0err+h1err,' <-----!'
      write(0,*) 'UH: u0tot=',u0tot,'u1tot=',u1tot,'duh=',u1tot-u0tot,'u0err=',u0err,'u1err=',u1err,'uh_err=',uh_err
      if (abs(u1tot-u0tot)>(u0err+u1err)+uh_err) &
        write(0,*) 'U non-conservation difference=',u1tot-u0tot,'allowed err=',u0err+u1err+uh_err,' <-----!'
      write(0,*) 'U: u0min=',u0min,'u1min=',u1min,'u2min=',u2min
      if (u1min<u0min) write(0,*) 'U minimum overshoot=',u1min-u0min,' <-----!'
      if (u2min<u0min) write(0,*) 'U2 minimum overshoot=',u2min-u0min,' <-----!'
      write(0,*) 'U: u0max=',u0max,'u1max=',u1max,'u2max=',u2max
      if (u1max>u0max) write(0,*) 'U maximum overshoot=',u1max-u0max,' <-----!'
      if (u2max>u0max) write(0,*) 'U2 maximum overshoot=',u2max-u0max,' <-----!'
      write(0,'(a3,6a24,2a3)') 'k','h0','left edge','u0','right edge','h1','u1','is','ie'
      do k = 1, max(n0,n1)
        if (k<=min(n0,n1)) then
          write(0,'(i3,1p6e24.16,2i3)') k,h0(k),ppoly0_e(k,1),u0(k),ppoly0_e(k,2),h1(k),u1(k),itgt_start(k),itgt_end(k)
        elseif (k>n0) then
          write(0,'(i3,96x,1p2e24.16,2i3)') k,h1(k),u1(k),itgt_start(k),itgt_end(k)
        else
          write(0,'(i3,1p4e24.16)') k,h0(k),ppoly0_e(k,1),u0(k),ppoly0_e(k,2)
        endif
      enddo
      write(0,'(a3,2a24)') 'k','u0','Polynomial coefficients'
      do k = 1, n0
        write(0,'(i3,1p6e24.16)') k,u0(k),ppoly0_coefs(k,:)
      enddo
      write(0,'(a3,3a24,a3,2a24)') 'k','Sub-cell h','Sub-cell u','Sub-cell hu','i0','xa','xb'
      xa = 0.
      dh0_eff = 0.
      do k = 1, n0+n1+1
        dh = h_sub(k)
        i0 = isub_src(k)
        dh0_eff = dh0_eff + dh ! Cumulative thickness within the source cell
        xb = dh0_eff / h0_eff(i0) ! This expression yields xa <= xb <= 1.0
        xb = min(1., xb) ! This is only needed when the total target column is wider than the source column
        write(0,'(i3,1p3e24.16,i3,1p2e24.16)') k,h_sub(k),u_sub(k),uh_sub(k),i0,xa,xb
        if (k<=n0+n1) then
          if (isub_src(k+1) /= i0) then
            dh0_eff = 0.; xa = 0.
          else
            xa = xb
          endif
        endif
      enddo
      call mom_error( fatal, 'MOM_remapping, remap_via_sub_cells: '//&
             'Remapping result is inconsistent!' )
    endif
    endif ! method<5
  endif ! debug_bounds
 
  ! Include the error remapping from source to sub-cells in the estimate of total remapping error
  uh_err = uh_err + u02_err
 
  if (present(ah_sub)) ah_sub(1:n0+n1+1) = h_sub(1:n0+n1+1)
  if (present(aisub_src)) aisub_src(1:n0+n1+1) = isub_src(1:n0+n1+1)
  if (present(aiss)) aiss(1:n0) = isrc_start(1:n0)
  if (present(aise)) aise(1:n0) = isrc_end(1:n0)
 

remapbydeltaz()

subroutine mom_remapping::remapbydeltaz ( integer, intent(in)  n0, real, dimension(:), intent(in)  h0, real, dimension(:), intent(in)  u0, real, dimension(:,:), intent(in)  ppoly0_E, real, dimension(:,:), intent(in)  ppoly0_coefs, integer, intent(in)  n1, real, dimension(:), intent(in)  dx1, integer, intent(in)  method, real, dimension(:), intent(out)  u1, real, dimension(:), intent(out), optional  h1, real, intent(in), optional  h_neglect  )

private

Remaps column of values u0 on grid h0 to implied grid h1 where the interfaces of h1 differ from those of h0 by dx. The new grid is defined relative to the original grid by change dx1(:) = xNew(:) - xOld(:) and the remapping calculated so that hNew(k) qNew(k) = hOld(k) qOld(k) + F(k+1) - F(k) where F(k) = dx1(k) qAverage and where qAverage is the average qOld in the region zOld(k) to zNew(k).

Parameters

[in]	n0	Number of cells in source grid
[in]	h0	Source grid sizes (size n0)
[in]	u0	Source cell averages (size n0)
[in]	ppoly0_e	Edge value of polynomial
[in]	ppoly0_coefs	Coefficients of polynomial
[in]	n1	Number of cells in target grid
[in]	dx1	Target grid edge positions (size n1+1)
[in]	method	Remapping scheme to use
[out]	u1	Target cell averages (size n1)
[out]	h1	Target grid widths (size n1)
[in]	h_neglect	A negligibly small width for the purpose of cell reconstructions in the same units as h.

Definition at line 1146 of file MOM_remapping.F90.

  integer,              intent(in)  :: n0     !< Number of cells in source grid
  real, dimension(:),   intent(in)  :: h0     !< Source grid sizes (size n0)
  real, dimension(:),   intent(in)  :: u0     !< Source cell averages (size n0)
  real, dimension(:,:), intent(in)  :: ppoly0_E !< Edge value of polynomial
  real, dimension(:,:), intent(in)  :: ppoly0_coefs !< Coefficients of polynomial
  integer,              intent(in)  :: n1     !< Number of cells in target grid
  real, dimension(:),   intent(in)  :: dx1    !< Target grid edge positions (size n1+1)
  integer,              intent(in)  :: method !< Remapping scheme to use
  real, dimension(:),   intent(out) :: u1     !< Target cell averages (size n1)
  real, dimension(:), &
              optional, intent(out) :: h1     !< Target grid widths (size n1)
  real,       optional, intent(in)  :: h_neglect !< A negligibly small width for the
                                           !! purpose of cell reconstructions
                                           !! in the same units as h.
  ! Local variables
  integer :: iTarget
  real    :: xL, xR    ! coordinates of target cell edges
  real    :: xOld, hOld, uOld
  real    :: xNew, hNew, h_err
  real    :: uhNew, hFlux, uAve, fluxL, fluxR
  integer :: jStart ! Used by integrateReconOnInterval()
  real    :: xStart ! Used by integrateReconOnInterval()
 
  ! Loop on cells in target grid. For each cell, iTarget, the left flux is
  ! the right flux of the cell to the left, iTarget-1.
  ! The left flux is initialized by started at iTarget=0 to calculate the
  ! right flux which can take into account the target left boundary being
  ! in the interior of the source domain.
  fluxr = 0.
  h_err = 0. ! For measuring round-off error
  jstart = 1
  xstart = 0.
  do itarget = 0,n1
    fluxl = fluxr ! This does nothing for iTarget=0
 
    if (itarget == 0) then
      xold = 0.     ! Left boundary is at x=0
      hold = -1.e30 ! Should not be used for iTarget = 0
      uold = -1.e30 ! Should not be used for iTarget = 0
    elseif (itarget <= n0) then
      xold = xold + h0(itarget) ! Position of right edge of cell
      hold = h0(itarget)
      uold = u0(itarget)
      h_err = h_err + epsilon(hold) * max(hold, xold)
    else
      hold = 0.       ! as if for layers>n0, they were vanished
      uold = 1.e30    ! and the initial value should not matter
    endif
    xnew = xold + dx1(itarget+1)
    xl = min( xold, xnew )
    xr = max( xold, xnew )
 
    ! hFlux is the positive width of the remapped volume
    hflux = abs(dx1(itarget+1))
    call integraterecononinterval( n0, h0, u0, ppoly0_e, ppoly0_coefs, method, &
                                   xl, xr, hflux, uave, jstart, xstart )
    ! uAve is the average value of u, independent of sign of dx1
    fluxr = dx1(itarget+1)*uave ! Includes sign of dx1
 
    if (itarget>0) then
      hnew = hold + ( dx1(itarget+1) - dx1(itarget) )
      hnew = max( 0., hnew )
      uhnew = ( uold * hold ) + ( fluxr - fluxl )
      if (hnew>0.) then
        u1(itarget) = uhnew / hnew
      else
        u1(itarget) = uave
      endif
      if (present(h1)) h1(itarget) = hnew
    endif
 
  enddo ! end iTarget loop on target grid cells
 

remapbyprojection()

subroutine mom_remapping::remapbyprojection ( integer, intent(in)  n0, real, dimension(:), intent(in)  h0, real, dimension(:), intent(in)  u0, real, dimension(:,:), intent(in)  ppoly0_E, real, dimension(:,:), intent(in)  ppoly0_coefs, integer, intent(in)  n1, real, dimension(:), intent(in)  h1, integer, intent(in)  method, real, dimension(:), intent(out)  u1, real, intent(in), optional  h_neglect  )

private

Remaps column of values u0 on grid h0 to grid h1 by integrating over the projection of each h1 cell onto the h0 grid.

Parameters

[in]	n0	Number of cells in source grid
[in]	h0	Source grid widths (size n0)
[in]	u0	Source cell averages (size n0)
[in]	ppoly0_e	Edge value of polynomial
[in]	ppoly0_coefs	Coefficients of polynomial
[in]	n1	Number of cells in target grid
[in]	h1	Target grid widths (size n1)
[in]	method	Remapping scheme to use
[out]	u1	Target cell averages (size n1)
[in]	h_neglect	A negligibly small width for the purpose of cell reconstructions in the same units as h.

Definition at line 1098 of file MOM_remapping.F90.

  integer,       intent(in)    :: n0     !< Number of cells in source grid
  real,          intent(in)    :: h0(:)  !< Source grid widths (size n0)
  real,          intent(in)    :: u0(:)  !< Source cell averages (size n0)
  real,          intent(in)    :: ppoly0_E(:,:)     !< Edge value of polynomial
  real,          intent(in)    :: ppoly0_coefs(:,:) !< Coefficients of polynomial
  integer,       intent(in)    :: n1     !< Number of cells in target grid
  real,          intent(in)    :: h1(:)  !< Target grid widths (size n1)
  integer,       intent(in)    :: method !< Remapping scheme to use
  real,          intent(out)   :: u1(:)  !< Target cell averages (size n1)
  real, optional, intent(in)   :: h_neglect !< A negligibly small width for the
                                           !! purpose of cell reconstructions
                                           !! in the same units as h.
  ! Local variables
  integer       :: iTarget
  real          :: xL, xR       ! coordinates of target cell edges
  integer       :: jStart ! Used by integrateReconOnInterval()
  real          :: xStart ! Used by integrateReconOnInterval()
 
  ! Loop on cells in target grid (grid1). For each target cell, we need to find
  ! in which source cells the target cell edges lie. The associated indexes are
  ! noted j0 and j1.
  xr = 0. ! Left boundary is at x=0
  jstart = 1
  xstart = 0.
  do itarget = 1,n1
    ! Determine the coordinates of the target cell edges
    xl = xr
    xr = xl + h1(itarget)
 
    call integraterecononinterval( n0, h0, u0, ppoly0_e, ppoly0_coefs, method, &
                                   xl, xr, h1(itarget), u1(itarget), jstart, xstart, h_neglect )
 
  enddo ! end iTarget loop on target grid cells
 

remapping_core_h()

subroutine, public mom_remapping::remapping_core_h	(	type(remapping_cs), intent(in)	CS,
		integer, intent(in)	n0,
		real, dimension(n0), intent(in)	h0,
		real, dimension(n0), intent(in)	u0,
		integer, intent(in)	n1,
		real, dimension(n1), intent(in)	h1,
		real, dimension(n1), intent(out)	u1,
		real, intent(in), optional	h_neglect,
		real, intent(in), optional	h_neglect_edge
	)

Remaps column of values u0 on grid h0 to grid h1 assuming the top edge is aligned.

Parameters

[in]	cs	Remapping control structure
[in]	n0	Number of cells on source grid
[in]	h0	Cell widths on source grid
[in]	u0	Cell averages on source grid
[in]	n1	Number of cells on target grid
[in]	h1	Cell widths on target grid
[out]	u1	Cell averages on target grid
[in]	h_neglect	A negligibly small width for the purpose of cell reconstructions in the same units as h0.
[in]	h_neglect_edge	A negligibly small width for the purpose of edge value calculations in the same units as h0.

Definition at line 189 of file MOM_remapping.F90.

  type(remapping_CS),  intent(in)  :: CS !< Remapping control structure
  integer,             intent(in)  :: n0 !< Number of cells on source grid
  real, dimension(n0), intent(in)  :: h0 !< Cell widths on source grid
  real, dimension(n0), intent(in)  :: u0 !< Cell averages on source grid
  integer,             intent(in)  :: n1 !< Number of cells on target grid
  real, dimension(n1), intent(in)  :: h1 !< Cell widths on target grid
  real, dimension(n1), intent(out) :: u1 !< Cell averages on target grid
  real, optional,      intent(in)  :: h_neglect !< A negligibly small width for the
                                         !! purpose of cell reconstructions
                                         !! in the same units as h0.
  real, optional,      intent(in)  :: h_neglect_edge !< A negligibly small width
                                         !! for the purpose of edge value
                                         !! calculations in the same units as h0.
  ! Local variables
  integer :: iMethod
  real, dimension(n0,2)           :: ppoly_r_E            !Edge value of polynomial
  real, dimension(n0,2)           :: ppoly_r_S            !Edge slope of polynomial
  real, dimension(n0,CS%degree+1) :: ppoly_r_coefs !Coefficients of polynomial
  integer :: k
  real :: eps, h0tot, h0err, h1tot, h1err, u0tot, u0err, u0min, u0max, u1tot, u1err, u1min, u1max, uh_err
  real :: hNeglect, hNeglect_edge
 
  hneglect = 1.0e-30 ; if (present(h_neglect)) hneglect = h_neglect
  hneglect_edge = 1.0e-10 ; if (present(h_neglect_edge)) hneglect_edge = h_neglect_edge
 
  call build_reconstructions_1d( cs, n0, h0, u0, ppoly_r_coefs, ppoly_r_e, ppoly_r_s, imethod, &
                                 hneglect, hneglect_edge )
 
  if (cs%check_reconstruction) call check_reconstructions_1d(n0, h0, u0, cs%degree, &
                                   cs%boundary_extrapolation, ppoly_r_coefs, ppoly_r_e, ppoly_r_s)
 
 
  call remap_via_sub_cells( n0, h0, u0, ppoly_r_e, ppoly_r_coefs, n1, h1, imethod, &
                            cs%force_bounds_in_subcell, u1, uh_err )
 
  if (cs%check_remapping) then
    ! Check errors and bounds
    call measure_input_bounds( n0, h0, u0, ppoly_r_e, h0tot, h0err, u0tot, u0err, u0min, u0max )
    call measure_output_bounds( n1, h1, u1, h1tot, h1err, u1tot, u1err, u1min, u1max )
    if (imethod<5) then ! We except PQM until we've debugged it
    if ( (abs(u1tot-u0tot)>(u0err+u1err)+uh_err .and. abs(h1tot-h0tot)<h0err+h1err) &
        .or. (u1min<u0min .or. u1max>u0max) ) then
      write(0,*) 'iMethod = ',imethod
      write(0,*) 'H: h0tot=',h0tot,'h1tot=',h1tot,'dh=',h1tot-h0tot,'h0err=',h0err,'h1err=',h1err
      if (abs(h1tot-h0tot)>h0err+h1err) &
        write(0,*) 'H non-conservation difference=',h1tot-h0tot,'allowed err=',h0err+h1err,' <-----!'
      write(0,*) 'UH: u0tot=',u0tot,'u1tot=',u1tot,'duh=',u1tot-u0tot,'u0err=',u0err,'u1err=',u1err,'uh_err=',uh_err
      if (abs(u1tot-u0tot)>(u0err+u1err)+uh_err) &
        write(0,*) 'U non-conservation difference=',u1tot-u0tot,'allowed err=',u0err+u1err+uh_err,' <-----!'
      write(0,*) 'U: u0min=',u0min,'u1min=',u1min
      if (u1min<u0min) write(0,*) 'U minimum overshoot=',u1min-u0min,' <-----!'
      write(0,*) 'U: u0max=',u0max,'u1max=',u1max
      if (u1max>u0max) write(0,*) 'U maximum overshoot=',u1max-u0max,' <-----!'
      write(0,'(a3,6a24)') 'k','h0','left edge','u0','right edge','h1','u1'
      do k = 1, max(n0,n1)
        if (k<=min(n0,n1)) then
          write(0,'(i3,1p6e24.16)') k,h0(k),ppoly_r_e(k,1),u0(k),ppoly_r_e(k,2),h1(k),u1(k)
        elseif (k>n0) then
          write(0,'(i3,96x,1p2e24.16)') k,h1(k),u1(k)
        else
          write(0,'(i3,1p4e24.16)') k,h0(k),ppoly_r_e(k,1),u0(k),ppoly_r_e(k,2)
        endif
      enddo
      write(0,'(a3,2a24)') 'k','u0','Polynomial coefficients'
      do k = 1, n0
        write(0,'(i3,1p6e24.16)') k,u0(k),ppoly_r_coefs(k,:)
      enddo
      call mom_error( fatal, 'MOM_remapping, remapping_core_h: '//&
             'Remapping result is inconsistent!' )
    endif
    endif ! method<5
  endif
 

remapping_core_w()

subroutine, public mom_remapping::remapping_core_w	(	type(remapping_cs), intent(in)	CS,
		integer, intent(in)	n0,
		real, dimension(n0), intent(in)	h0,
		real, dimension(n0), intent(in)	u0,
		integer, intent(in)	n1,
		real, dimension(n1+1), intent(in)	dx,
		real, dimension(n1), intent(out)	u1,
		real, intent(in), optional	h_neglect,
		real, intent(in), optional	h_neglect_edge
	)

Remaps column of values u0 on grid h0 to implied grid h1 where the interfaces of h1 differ from those of h0 by dx.

Parameters

[in]	cs	Remapping control structure
[in]	n0	Number of cells on source grid
[in]	h0	Cell widths on source grid
[in]	u0	Cell averages on source grid
[in]	n1	Number of cells on target grid
[in]	dx	Cell widths on target grid
[out]	u1	Cell averages on target grid
[in]	h_neglect	A negligibly small width for the purpose of cell reconstructions in the same units as h0.
[in]	h_neglect_edge	A negligibly small width for the purpose of edge value calculations in the same units as h0.

Definition at line 267 of file MOM_remapping.F90.

  type(remapping_CS),  intent(in)  :: CS !< Remapping control structure
  integer,             intent(in)  :: n0 !< Number of cells on source grid
  real, dimension(n0), intent(in)  :: h0 !< Cell widths on source grid
  real, dimension(n0), intent(in)  :: u0 !< Cell averages on source grid
  integer,             intent(in)  :: n1 !< Number of cells on target grid
  real, dimension(n1+1), intent(in) :: dx !< Cell widths on target grid
  real, dimension(n1), intent(out) :: u1 !< Cell averages on target grid
  real, optional,      intent(in)  :: h_neglect !< A negligibly small width for the
                                         !! purpose of cell reconstructions
                                         !! in the same units as h0.
  real, optional,      intent(in)  :: h_neglect_edge !< A negligibly small width
                                         !! for the purpose of edge value
                                         !! calculations in the same units as h0.
  ! Local variables
  integer :: iMethod
  real, dimension(n0,2)           :: ppoly_r_E            !Edge value of polynomial
  real, dimension(n0,2)           :: ppoly_r_S            !Edge slope of polynomial
  real, dimension(n0,CS%degree+1) :: ppoly_r_coefs !Coefficients of polynomial
  integer :: k
  real :: eps, h0tot, h0err, h1tot, h1err
  real :: u0tot, u0err, u0min, u0max, u1tot, u1err, u1min, u1max, uh_err
  real, dimension(n1) :: h1 !< Cell widths on target grid
  real :: hNeglect, hNeglect_edge
 
  hneglect = 1.0e-30 ; if (present(h_neglect)) hneglect = h_neglect
  hneglect_edge = 1.0e-10 ; if (present(h_neglect_edge)) hneglect_edge = h_neglect_edge
 
  call build_reconstructions_1d( cs, n0, h0, u0, ppoly_r_coefs, ppoly_r_e, ppoly_r_s, imethod,&
                                  hneglect, hneglect_edge )
 
  if (cs%check_reconstruction) call check_reconstructions_1d(n0, h0, u0, cs%degree, &
                                   cs%boundary_extrapolation, ppoly_r_coefs, ppoly_r_e, ppoly_r_s)
 
  ! This is a temporary step prior to switching to remapping_core_h()
  do k = 1, n1
    if (k<=n0) then
      h1(k) = max( 0., h0(k) + ( dx(k+1) - dx(k) ) )
    else
      h1(k) = max( 0., dx(k+1) - dx(k) )
    endif
  enddo
  call remap_via_sub_cells( n0, h0, u0, ppoly_r_e, ppoly_r_coefs, n1, h1, imethod, &
                            cs%force_bounds_in_subcell,u1, uh_err )
! call remapByDeltaZ( n0, h0, u0, ppoly_r_E, ppoly_r_coefs, n1, dx, iMethod, u1, hNeglect )
! call remapByProjection( n0, h0, u0, CS%ppoly_r, n1, h1, iMethod, u1, hNeglect )
 
  if (cs%check_remapping) then
    ! Check errors and bounds
    call measure_input_bounds( n0, h0, u0, ppoly_r_e, h0tot, h0err, u0tot, u0err, u0min, u0max )
    call measure_output_bounds( n1, h1, u1, h1tot, h1err, u1tot, u1err, u1min, u1max )
    if (imethod<5) then ! We except PQM until we've debugged it
    if ( (abs(u1tot-u0tot)>(u0err+u1err)+uh_err .and. abs(h1tot-h0tot)<h0err+h1err) &
        .or. (u1min<u0min .or. u1max>u0max) ) then
      write(0,*) 'iMethod = ',imethod
      write(0,*) 'H: h0tot=',h0tot,'h1tot=',h1tot,'dh=',h1tot-h0tot,'h0err=',h0err,'h1err=',h1err
      if (abs(h1tot-h0tot)>h0err+h1err) &
        write(0,*) 'H non-conservation difference=',h1tot-h0tot,'allowed err=',h0err+h1err,' <-----!'
      write(0,*) 'UH: u0tot=',u0tot,'u1tot=',u1tot,'duh=',u1tot-u0tot,'u0err=',u0err,'u1err=',u1err,'uh_err=',uh_err
      if (abs(u1tot-u0tot)>(u0err+u1err)+uh_err) &
        write(0,*) 'U non-conservation difference=',u1tot-u0tot,'allowed err=',u0err+u1err+uh_err,' <-----!'
      write(0,*) 'U: u0min=',u0min,'u1min=',u1min
      if (u1min<u0min) write(0,*) 'U minimum overshoot=',u1min-u0min,' <-----!'
      write(0,*) 'U: u0max=',u0max,'u1max=',u1max
      if (u1max>u0max) write(0,*) 'U maximum overshoot=',u1max-u0max,' <-----!'
      write(0,'(a3,6a24)') 'k','h0','left edge','u0','right edge','h1','u1'
      do k = 1, max(n0,n1)
        if (k<=min(n0,n1)) then
          write(0,'(i3,1p6e24.16)') k,h0(k),ppoly_r_e(k,1),u0(k),ppoly_r_e(k,2),h1(k),u1(k)
        elseif (k>n0) then
          write(0,'(i3,96x,1p2e24.16)') k,h1(k),u1(k)
        else
          write(0,'(i3,1p4e24.16)') k,h0(k),ppoly_r_e(k,1),u0(k),ppoly_r_e(k,2)
        endif
      enddo
      write(0,'(a3,2a24)') 'k','u0','Polynomial coefficients'
      do k = 1, n0
        write(0,'(i3,1p6e24.16)') k,u0(k),ppoly_r_coefs(k,:)
      enddo
      call mom_error( fatal, 'MOM_remapping, remapping_core_w: '//&
             'Remapping result is inconsistent!' )
    endif
    endif ! method<5
  endif
 

remapping_set_param()

subroutine, public mom_remapping::remapping_set_param	(	type(remapping_cs), intent(inout)	CS,
		character(len=*), intent(in), optional	remapping_scheme,
		logical, intent(in), optional	boundary_extrapolation,
		logical, intent(in), optional	check_reconstruction,
		logical, intent(in), optional	check_remapping,
		logical, intent(in), optional	force_bounds_in_subcell,
		logical, intent(in), optional	answers_2018
	)

Set parameters within remapping object.

Parameters

[in,out]	cs	Remapping control structure
[in]	remapping_scheme	Remapping scheme to use
[in]	boundary_extrapolation	Indicate to extrapolate in boundary cells
[in]	check_reconstruction	Indicate to check reconstructions
[in]	check_remapping	Indicate to check results of remapping
[in]	force_bounds_in_subcell	Force subcells values to be bounded
[in]	answers_2018	If true use older, less acccurate expressions.

Definition at line 90 of file MOM_remapping.F90.

  type(remapping_CS),         intent(inout) :: CS !< Remapping control structure
  character(len=*), optional, intent(in)    :: remapping_scheme !< Remapping scheme to use
  logical, optional,          intent(in)    :: boundary_extrapolation !< Indicate to extrapolate in boundary cells
  logical, optional,          intent(in)    :: check_reconstruction !< Indicate to check reconstructions
  logical, optional,          intent(in)    :: check_remapping !< Indicate to check results of remapping
  logical, optional,          intent(in)    :: force_bounds_in_subcell !< Force subcells values to be bounded
  logical, optional,          intent(in)    :: answers_2018 !< If true use older, less acccurate expressions.
 
  if (present(remapping_scheme)) then
    call setreconstructiontype( remapping_scheme, cs )
  endif
  if (present(boundary_extrapolation)) then
    cs%boundary_extrapolation = boundary_extrapolation
  endif
  if (present(check_reconstruction)) then
    cs%check_reconstruction = check_reconstruction
  endif
  if (present(check_remapping)) then
    cs%check_remapping = check_remapping
  endif
  if (present(force_bounds_in_subcell)) then
    cs%force_bounds_in_subcell = force_bounds_in_subcell
  endif
  if (present(answers_2018)) then
    cs%answers_2018 = answers_2018
  endif

remapping_unit_tests()

logical function, public mom_remapping::remapping_unit_tests

(

logical, intent(in)

verbose

)

Runs unit tests on remapping functions. Should only be called from a single/root thread Returns True if a test fails, otherwise False.

Parameters

[in]

verbose

If true, write results to stdout

Definition at line 1616 of file MOM_remapping.F90.

  logical, intent(in) :: verbose !< If true, write results to stdout
  ! Local variables
  integer, parameter :: n0 = 4, n1 = 3, n2 = 6
  real :: h0(n0), x0(n0+1), u0(n0)
  real :: h1(n1), x1(n1+1), u1(n1), hn1(n1), dx1(n1+1)
  real :: h2(n2), x2(n2+1), u2(n2), hn2(n2), dx2(n2+1)
  data u0 /9., 3., -3., -9./   ! Linear profile, 4 at surface to -4 at bottom
  data h0 /4*0.75/ ! 4 uniform layers with total depth of 3
  data h1 /3*1./   ! 3 uniform layers with total depth of 3
  data h2 /6*0.5/  ! 6 uniform layers with total depth of 3
  type(remapping_CS) :: CS !< Remapping control structure
  real, allocatable, dimension(:,:) :: ppoly0_E, ppoly0_S, ppoly0_coefs
  logical :: answers_2018 !  If true use older, less acccurate expressions.
  integer :: i
  real :: err, h_neglect, h_neglect_edge
  logical :: thisTest, v
 
  v = verbose
  answers_2018 = .false. ! .true.
  h_neglect = hneglect_dflt
  h_neglect_edge = hneglect_dflt ; if (answers_2018) h_neglect_edge = 1.0e-10
 
  write(*,*) '==== MOM_remapping: remapping_unit_tests ================='
  remapping_unit_tests = .false. ! Normally return false
 
  thistest = .false.
  call buildgridfromh(n0, h0, x0)
  do i=1,n0+1
    err=x0(i)-0.75*real(i-1)
    if (abs(err)>real(i-1)*epsilon(err)) thistest = .true.
  enddo
  if (thistest) write(*,*) 'remapping_unit_tests: Failed buildGridFromH() 1'
  remapping_unit_tests = remapping_unit_tests .or. thistest
  call buildgridfromh(n1, h1, x1)
  do i=1,n1+1
    err=x1(i)-real(i-1)
    if (abs(err)>real(i-1)*epsilon(err)) thistest = .true.
  enddo
  if (thistest) write(*,*) 'remapping_unit_tests: Failed buildGridFromH() 2'
  remapping_unit_tests = remapping_unit_tests .or. thistest
 
  thistest = .false.
  call initialize_remapping(cs, 'PPM_H4', answers_2018=answers_2018)
  if (verbose) write(*,*) 'h0 (test data)'
  if (verbose) call dumpgrid(n0,h0,x0,u0)
 
  call dzfromh1h2( n0, h0, n1, h1, dx1 )
  call remapping_core_w( cs, n0, h0, u0, n1, dx1, u1, h_neglect, h_neglect_edge)
  do i=1,n1
    err=u1(i)-8.*(0.5*real(1+n1)-real(i))
    if (abs(err)>real(n1-1)*epsilon(err)) thistest = .true.
  enddo
  if (verbose) write(*,*) 'h1 (by projection)'
  if (verbose) call dumpgrid(n1,h1,x1,u1)
  if (thistest) write(*,*) 'remapping_unit_tests: Failed remapping_core_w()'
  remapping_unit_tests = remapping_unit_tests .or. thistest
 
  thistest = .false.
  allocate(ppoly0_e(n0,2))
  allocate(ppoly0_s(n0,2))
  allocate(ppoly0_coefs(n0,cs%degree+1))
 
  ppoly0_e(:,:) = 0.0
  ppoly0_s(:,:) = 0.0
  ppoly0_coefs(:,:) = 0.0
 
  call edge_values_explicit_h4( n0, h0, u0, ppoly0_e, h_neglect=1e-10, answers_2018=answers_2018 )
  call ppm_reconstruction( n0, h0, u0, ppoly0_e, ppoly0_coefs, h_neglect, answers_2018=answers_2018 )
  call ppm_boundary_extrapolation( n0, h0, u0, ppoly0_e, ppoly0_coefs, h_neglect )
  u1(:) = 0.
  call remapbyprojection( n0, h0, u0, ppoly0_e, ppoly0_coefs, &
                          n1, h1, integration_ppm, u1, h_neglect )
  do i=1,n1
    err=u1(i)-8.*(0.5*real(1+n1)-real(i))
    if (abs(err)>2.*epsilon(err)) thistest = .true.
  enddo
  if (thistest) write(*,*) 'remapping_unit_tests: Failed remapByProjection()'
  remapping_unit_tests = remapping_unit_tests .or. thistest
 
  thistest = .false.
  u1(:) = 0.
  call remapbydeltaz( n0, h0, u0, ppoly0_e, ppoly0_coefs, &
                      n1, x1-x0(1:n1+1), &
                      integration_ppm, u1, hn1, h_neglect )
  if (verbose) write(*,*) 'h1 (by delta)'
  if (verbose) call dumpgrid(n1,h1,x1,u1)
  hn1=hn1-h1
  do i=1,n1
    err=u1(i)-8.*(0.5*real(1+n1)-real(i))
    if (abs(err)>2.*epsilon(err)) thistest = .true.
  enddo
  if (thistest) write(*,*) 'remapping_unit_tests: Failed remapByDeltaZ() 1'
  remapping_unit_tests = remapping_unit_tests .or. thistest
 
  thistest = .false.
  call buildgridfromh(n2, h2, x2)
  dx2(1:n0+1) = x2(1:n0+1) - x0
  dx2(n0+2:n2+1) = x2(n0+2:n2+1) - x0(n0+1)
  call remapbydeltaz( n0, h0, u0, ppoly0_e, ppoly0_coefs, &
                      n2, dx2, &
                      integration_ppm, u2, hn2, h_neglect )
  if (verbose) write(*,*) 'h2'
  if (verbose) call dumpgrid(n2,h2,x2,u2)
  if (verbose) write(*,*) 'hn2'
  if (verbose) call dumpgrid(n2,hn2,x2,u2)
 
  do i=1,n2
    err=u2(i)-8./2.*(0.5*real(1+n2)-real(i))
    if (abs(err)>2.*epsilon(err)) thistest = .true.
  enddo
  if (thistest) write(*,*) 'remapping_unit_tests: Failed remapByDeltaZ() 2'
  remapping_unit_tests = remapping_unit_tests .or. thistest
 
  if (verbose) write(*,*) 'Via sub-cells'
  thistest = .false.
  call remap_via_sub_cells( n0, h0, u0, ppoly0_e, ppoly0_coefs, &
                            n2, h2, integration_ppm, .false., u2, err )
  if (verbose) call dumpgrid(n2,h2,x2,u2)
 
  do i=1,n2
    err=u2(i)-8./2.*(0.5*real(1+n2)-real(i))
    if (abs(err)>2.*epsilon(err)) thistest = .true.
  enddo
  if (thistest) write(*,*) 'remapping_unit_tests: Failed remap_via_sub_cells() 2'
  remapping_unit_tests = remapping_unit_tests .or. thistest
 
  call remap_via_sub_cells( n0, h0, u0, ppoly0_e, ppoly0_coefs, &
                            6, (/.125,.125,.125,.125,.125,.125/), integration_ppm, .false., u2, err )
  if (verbose) call dumpgrid(6,h2,x2,u2)
 
  call remap_via_sub_cells( n0, h0, u0, ppoly0_e, ppoly0_coefs, &
                            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(*,*) 'Pass'
 
  write(*,*) '===== 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))
 
  call pcm_reconstruction(3, (/1.,2.,4./), ppoly0_e(1:3,:), &
                          ppoly0_coefs(1:3,:) )
  remapping_unit_tests = remapping_unit_tests .or. &
    test_answer(v, 3, ppoly0_e(:,1), (/1.,2.,4./), 'PCM: left edges')
  remapping_unit_tests = remapping_unit_tests .or. &
    test_answer(v, 3, ppoly0_e(:,2), (/1.,2.,4./), 'PCM: right edges')
  remapping_unit_tests = remapping_unit_tests .or. &
    test_answer(v, 3, ppoly0_coefs(:,1), (/1.,2.,4./), 'PCM: P0')
 
  call plm_reconstruction(3, (/1.,1.,1./), (/1.,3.,5./), ppoly0_e(1:3,:), &
                          ppoly0_coefs(1:3,:), h_neglect )
  remapping_unit_tests = remapping_unit_tests .or. &
    test_answer(v, 3, ppoly0_e(:,1), (/1.,2.,5./), 'Unlim PLM: left edges')
  remapping_unit_tests = remapping_unit_tests .or. &
    test_answer(v, 3, ppoly0_e(:,2), (/1.,4.,5./), 'Unlim PLM: right edges')
  remapping_unit_tests = remapping_unit_tests .or. &
    test_answer(v, 3, ppoly0_coefs(:,1), (/1.,2.,5./), 'Unlim PLM: P0')
  remapping_unit_tests = remapping_unit_tests .or. &
    test_answer(v, 3, ppoly0_coefs(:,2), (/0.,2.,0./), 'Unlim PLM: P1')
 
  call plm_reconstruction(3, (/1.,1.,1./), (/1.,2.,7./), ppoly0_e(1:3,:), &
                          ppoly0_coefs(1:3,:), h_neglect )
  remapping_unit_tests = remapping_unit_tests .or. &
    test_answer(v, 3, ppoly0_e(:,1), (/1.,1.,7./), 'Left lim PLM: left edges')
  remapping_unit_tests = remapping_unit_tests .or. &
    test_answer(v, 3, ppoly0_e(:,2), (/1.,3.,7./), 'Left lim PLM: right edges')
  remapping_unit_tests = remapping_unit_tests .or. &
    test_answer(v, 3, ppoly0_coefs(:,1), (/1.,1.,7./), 'Left lim PLM: P0')
  remapping_unit_tests = remapping_unit_tests .or. &
    test_answer(v, 3, ppoly0_coefs(:,2), (/0.,2.,0./), 'Left lim PLM: P1')
 
  call plm_reconstruction(3, (/1.,1.,1./), (/1.,6.,7./), ppoly0_e(1:3,:), &
                          ppoly0_coefs(1:3,:), h_neglect )
  remapping_unit_tests = remapping_unit_tests .or. &
    test_answer(v, 3, ppoly0_e(:,1), (/1.,5.,7./), 'Right lim PLM: left edges')
  remapping_unit_tests = remapping_unit_tests .or. &
    test_answer(v, 3, ppoly0_e(:,2), (/1.,7.,7./), 'Right lim PLM: right edges')
  remapping_unit_tests = remapping_unit_tests .or. &
    test_answer(v, 3, ppoly0_coefs(:,1), (/1.,5.,7./), 'Right lim PLM: P0')
  remapping_unit_tests = remapping_unit_tests .or. &
    test_answer(v, 3, ppoly0_coefs(:,2), (/0.,2.,0./), 'Right lim PLM: P1')
 
  call plm_reconstruction(3, (/1.,2.,3./), (/1.,4.,9./), ppoly0_e(1:3,:), &
                          ppoly0_coefs(1:3,:), h_neglect )
  remapping_unit_tests = remapping_unit_tests .or. &
    test_answer(v, 3, ppoly0_e(:,1), (/1.,2.,9./), 'Non-uniform line PLM: left edges')
  remapping_unit_tests = remapping_unit_tests .or. &
    test_answer(v, 3, ppoly0_e(:,2), (/1.,6.,9./), 'Non-uniform line PLM: right edges')
  remapping_unit_tests = remapping_unit_tests .or. &
    test_answer(v, 3, ppoly0_coefs(:,1), (/1.,2.,9./), 'Non-uniform line PLM: P0')
  remapping_unit_tests = remapping_unit_tests .or. &
    test_answer(v, 3, ppoly0_coefs(:,2), (/0.,4.,0./), 'Non-uniform line PLM: P1')
 
  call edge_values_explicit_h4( 5, (/1.,1.,1.,1.,1./), (/1.,3.,5.,7.,9./), ppoly0_e, &
                                h_neglect=1e-10, answers_2018=answers_2018 )
  ! The next two tests currently fail due to roundoff, but pass when given a reasonable tolerance.
  thistest = test_answer(v, 5, ppoly0_e(:,1), (/0.,2.,4.,6.,8./), 'Line H4: left edges', tol=8.0e-15)
  remapping_unit_tests = remapping_unit_tests .or. thistest
  thistest = test_answer(v, 5, ppoly0_e(:,2), (/2.,4.,6.,8.,10./), 'Line H4: right edges', tol=1.0e-14)
  remapping_unit_tests = remapping_unit_tests .or. thistest
  ppoly0_e(:,1) = (/0.,2.,4.,6.,8./)
  ppoly0_e(:,2) = (/2.,4.,6.,8.,10./)
  call ppm_reconstruction(5, (/1.,1.,1.,1.,1./), (/1.,3.,5.,7.,9./), ppoly0_e(1:5,:), &
                              ppoly0_coefs(1:5,:), h_neglect, answers_2018=answers_2018 )
  remapping_unit_tests = remapping_unit_tests .or. &
    test_answer(v, 5, ppoly0_coefs(:,1), (/1.,2.,4.,6.,9./), 'Line PPM: P0')
  remapping_unit_tests = remapping_unit_tests .or. &
    test_answer(v, 5, ppoly0_coefs(:,2), (/0.,2.,2.,2.,0./), 'Line PPM: P1')
  remapping_unit_tests = remapping_unit_tests .or. &
    test_answer(v, 5, ppoly0_coefs(:,3), (/0.,0.,0.,0.,0./), 'Line PPM: P2')
 
  call edge_values_explicit_h4( 5, (/1.,1.,1.,1.,1./), (/1.,1.,7.,19.,37./), ppoly0_e, &
                                h_neglect=1e-10, answers_2018=answers_2018 )
  ! The next two tests are now passing when answers_2018 = .false., but otherwise only work to roundoff.
  thistest = test_answer(v, 5, ppoly0_e(:,1), (/3.,0.,3.,12.,27./), 'Parabola H4: left edges', tol=2.7e-14)
  remapping_unit_tests = remapping_unit_tests .or. thistest
  thistest = test_answer(v, 5, ppoly0_e(:,2), (/0.,3.,12.,27.,48./), 'Parabola H4: right edges', tol=4.8e-14)
  remapping_unit_tests = remapping_unit_tests .or. thistest
  ppoly0_e(:,1) = (/0.,0.,3.,12.,27./)
  ppoly0_e(:,2) = (/0.,3.,12.,27.,48./)
  call ppm_reconstruction(5, (/1.,1.,1.,1.,1./), (/0.,1.,7.,19.,37./), ppoly0_e(1:5,:), &
                          ppoly0_coefs(1:5,:), h_neglect, answers_2018=answers_2018 )
  remapping_unit_tests = remapping_unit_tests .or. &
    test_answer(v, 5, ppoly0_e(:,1), (/0.,0.,3.,12.,37./), 'Parabola PPM: left edges')
  remapping_unit_tests = remapping_unit_tests .or. &
    test_answer(v, 5, ppoly0_e(:,2), (/0.,3.,12.,27.,37./), 'Parabola PPM: right edges')
  remapping_unit_tests = remapping_unit_tests .or. &
    test_answer(v, 5, ppoly0_coefs(:,1), (/0.,0.,3.,12.,37./), 'Parabola PPM: P0')
  remapping_unit_tests = remapping_unit_tests .or. &
    test_answer(v, 5, ppoly0_coefs(:,2), (/0.,0.,6.,12.,0./), 'Parabola PPM: P1')
  remapping_unit_tests = remapping_unit_tests .or. &
    test_answer(v, 5, ppoly0_coefs(:,3), (/0.,3.,3.,3.,0./), 'Parabola PPM: P2')
 
  ppoly0_e(:,1) = (/0.,0.,6.,10.,15./)
  ppoly0_e(:,2) = (/0.,6.,12.,17.,15./)
  call ppm_reconstruction(5, (/1.,1.,1.,1.,1./), (/0.,5.,7.,16.,15./), ppoly0_e(1:5,:), &
                          ppoly0_coefs(1:5,:), h_neglect, answers_2018=answers_2018 )
  remapping_unit_tests = remapping_unit_tests .or. &
    test_answer(v, 5, ppoly0_e(:,1), (/0.,3.,6.,16.,15./), 'Limits PPM: left edges')
  remapping_unit_tests = remapping_unit_tests .or. &
    test_answer(v, 5, ppoly0_e(:,2), (/0.,6.,9.,16.,15./), 'Limits PPM: right edges')
  remapping_unit_tests = remapping_unit_tests .or. &
    test_answer(v, 5, ppoly0_coefs(:,1), (/0.,3.,6.,16.,15./), 'Limits PPM: P0')
  remapping_unit_tests = remapping_unit_tests .or. &
    test_answer(v, 5, ppoly0_coefs(:,2), (/0.,6.,0.,0.,0./), 'Limits PPM: P1')
  remapping_unit_tests = remapping_unit_tests .or. &
    test_answer(v, 5, ppoly0_coefs(:,3), (/0.,-3.,3.,0.,0./), 'Limits PPM: P2')
 
  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. &
    test_answer(v, 4, ppoly0_e(1:4,1), (/5.,5.,3.,1./), 'PPM: left edges h=0110')
  remapping_unit_tests = remapping_unit_tests .or. &
    test_answer(v, 4, ppoly0_e(1:4,2), (/5.,3.,1.,1./), 'PPM: right edges h=0110')
  call remap_via_sub_cells( 4, (/0.,1.,1.,0./), (/5.,4.,2.,1./), ppoly0_e(1:4,:), &
                            ppoly0_coefs(1:4,:), &
                            2, (/1.,1./), integration_plm, .false., u2, err )
  remapping_unit_tests = remapping_unit_tests .or. &
    test_answer(v, 2, u2, (/4.,2./), 'PLM: remapped  h=0110->h=11')
 
  deallocate(ppoly0_e, ppoly0_s, ppoly0_coefs)
 
  if (.not. remapping_unit_tests) write(*,*) 'Pass'
 

setreconstructiontype()

subroutine mom_remapping::setreconstructiontype ( character(len=*), intent(in)  string, type(remapping_cs), intent(inout)  CS  )

private

Changes the method of reconstruction Use this routine to parse a string parameter specifying the reconstruction and re-allocates work arrays appropriately. It is called from initialize_remapping but can be called from an external module too.

Parameters

[in]	string	String to parse for method
[in,out]	cs	Remapping control structure

Definition at line 1571 of file MOM_remapping.F90.

  character(len=*),   intent(in)    :: string !< String to parse for method
  type(remapping_CS), intent(inout) :: CS !< Remapping control structure
  ! Local variables
  integer :: degree
  degree = -99
  select case ( uppercase(trim(string)) )
    case ("PCM")
      cs%remapping_scheme = remapping_pcm
      degree = 0
    case ("PLM")
      cs%remapping_scheme = remapping_plm
      degree = 1
    case ("PPM_H4")
      cs%remapping_scheme = remapping_ppm_h4
      degree = 2
    case ("PPM_IH4")
      cs%remapping_scheme = remapping_ppm_ih4
      degree = 2
    case ("PQM_IH4IH3")
      cs%remapping_scheme = remapping_pqm_ih4ih3
      degree = 4
    case ("PQM_IH6IH5")
      cs%remapping_scheme = remapping_pqm_ih6ih5
      degree = 4
    case default
      call mom_error(fatal, "setReconstructionType: "//&
       "Unrecognized choice for REMAPPING_SCHEME ("//trim(string)//").")
  end select
 
  cs%degree = degree
 

test_answer()

logical function mom_remapping::test_answer ( logical, intent(in)  verbose, integer, intent(in)  n, real, dimension(n), intent(in)  u, real, dimension(n), intent(in)  u_true, character(len=*), intent(in)  label, real, intent(in), optional  tol  )

private

Returns true if any cell of u and u_true are not identical. Returns false otherwise.

Parameters

[in]	verbose	If true, write results to stdout
[in]	n	Number of cells in u
[in]	u	Values to test
[in]	u_true	Values to test against (correct answer)
[in]	label	Message
[in]	tol	The tolerance for differences between u and u_true

Definition at line 1887 of file MOM_remapping.F90.

  logical,            intent(in) :: verbose !< If true, write results to stdout
  integer,            intent(in) :: n      !< Number of cells in u
  real, dimension(n), intent(in) :: u      !< Values to test
  real, dimension(n), intent(in) :: u_true !< Values to test against (correct answer)
  character(len=*),   intent(in) :: label  !< Message
  real,     optional, intent(in) :: tol    !< The tolerance for differences between u and u_true
  ! Local variables
  real :: tolerance ! The tolerance for differences between u and u_true
  integer :: k
 
  tolerance = 0.0 ; if (present(tol)) tolerance = tol
  test_answer = .false.
  do k = 1, n
    if (abs(u(k) - u_true(k)) > tolerance) test_answer = .true.
  enddo
  if (test_answer .or. verbose) then
    write(stdout,'(a4,2a24,x,a)') 'k','Calculated value','Correct value',label
    do k = 1, n
      if (abs(u(k) - u_true(k)) > tolerance) then
        write(stdout,'(i4,1p2e24.16,a,1pe24.16,a)') k,u(k),u_true(k),' err=',u(k)-u_true(k),' < wrong'
        write(stderr,'(i4,1p2e24.16,a,1pe24.16,a)') k,u(k),u_true(k),' err=',u(k)-u_true(k),' < wrong'
      else
        write(stdout,'(i4,1p2e24.16)') k,u(k),u_true(k)
      endif
    enddo
  endif