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)
	Set parameters within remapping object. More...
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)
	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)
	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	build_reconstructions_1d (n0, h0, u0, remapping_scheme, deg, boundary_extrapolation, ppoly_r_coefficients, ppoly_r_E, ppoly_r_S, iMethod)
	Creates polynomial reconstructions of u0 on the source grid h0. More...
subroutine	check_reconstructions_1d (n0, h0, u0, deg, boundary_extrapolation, ppoly_r_coefficients, 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_coefficients, 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. More...
real function	average_value_ppoly (n0, u0, ppoly0_E, ppoly0_coefficients, 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, ppoly_E, 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_coefficients, n1, h1, method, u1)
	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_coefficients, n1, dx1, method, u1, h1)
	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_coefficients, method, xL, xR, hC, uAve, jStart, xStart)
	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)
	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)
	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. More...
integer, parameter	remapping_plm = 1
	O(h^2) remapping scheme. More...
integer, parameter	remapping_ppm_h4 = 2
	O(h^3) remapping scheme. More...
integer, parameter	remapping_ppm_ih4 = 3
	O(h^3) remapping scheme. More...
integer, parameter	remapping_pqm_ih4ih3 = 4
	O(h^4) remapping scheme. More...
integer, parameter	remapping_pqm_ih6ih5 = 5
	O(h^5) remapping scheme. More...
integer, parameter	integration_pcm = 0
	Piecewise Constant Method. More...
integer, parameter	integration_plm = 1
	Piecewise Linear Method. More...
integer, parameter	integration_ppm = 3
	Piecewise Parabolic Method. More...
integer, parameter	integration_pqm = 5
	Piecewise Quartic Method. More...
character(len=40)	mdl = "MOM_remapping"
	This module's name. More...
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. More...
character(len=3), public	remappingdefaultscheme = "PLM"
	Default remapping method. More...
real, parameter	h_neglect = 1.E-30
	A dimensional (H units) number that can be added to thicknesses in a denominator without changing the numerical result, except where a division by zero would otherwise occur. More...
logical, parameter	old_algorithm = .false.
	Use the old "broken" algorithm. This is a temporary measure to assist debugging until we delete the old algorithm. More...

Function/Subroutine Documentation

average_value_ppoly()

real function 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_coefficients, integer, intent(in)  method, integer, intent(in)  i0, real, intent(in)  xa, real, intent(in)  xb  )

private

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_coefficients	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 860 of file MOM_remapping.F90.

References integration_pcm, integration_plm, integration_ppm, integration_pqm, and mom_error_handler::mom_error().

Referenced by remap_via_sub_cells().

  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_coefficients(:,:) !< 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_coefficients(i0,1)                       &
          + ppoly0_coefficients(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_coefficients(i0,1)                                                         &
          + ( ppoly0_coefficients(i0,2) * 0.5 * ( xapxb )                                       &
          + ( ppoly0_coefficients(i0,3) * r_3 * ( xa2pxb2 + xa*xb )                             &
          + ( ppoly0_coefficients(i0,4) * 0.25* ( xa2pxb2 * xapxb )                             &
          +   ppoly0_coefficients(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_coefficients(i0,1)
      case ( integration_plm )
       !u_ave =        ppoly0_coefficients(i0,1)   &
       !      + xa *   ppoly0_coefficients(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_coefficients(i0,1)   &
       !      + xa * ( ppoly0_coefficients(i0,2)   &
       !      + xa *   ppoly0_coefficients(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_coefficients(i0,1)   &
              + xa * ( ppoly0_coefficients(i0,2)   &
              + xa * ( ppoly0_coefficients(i0,3)   &
              + xa * ( ppoly0_coefficients(i0,4)   &
              + xa *   ppoly0_coefficients(i0,5) ) ) )
      case default
        call mom_error( fatal,'The selected integration method is invalid' )
    end select
  endif
  average_value_ppoly = u_ave

Here is the call graph for this function:

Here is the caller graph for this function:

build_reconstructions_1d()

subroutine mom_remapping::build_reconstructions_1d ( integer, intent(in)  n0, real, dimension(n0), intent(in)  h0, real, dimension(n0), intent(in)  u0, integer, intent(in)  remapping_scheme, integer, intent(in)  deg, logical, intent(in)  boundary_extrapolation, real, dimension(n0,deg+1), intent(out)  ppoly_r_coefficients, real, dimension(n0,2), intent(out)  ppoly_r_E, real, dimension(n0,2), intent(out)  ppoly_r_S, integer, intent(out)  iMethod  )

private

Creates polynomial reconstructions of u0 on the source grid h0.

Parameters

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

Definition at line 306 of file MOM_remapping.F90.

References regrid_edge_slopes::edge_slopes_implicit_h3(), regrid_edge_slopes::edge_slopes_implicit_h5(), regrid_edge_values::edge_values_implicit_h4(), regrid_edge_values::edge_values_implicit_h6(), integration_pcm, integration_plm, integration_ppm, integration_pqm, mom_error_handler::mom_error(), pcm_functions::pcm_reconstruction(), plm_functions::plm_boundary_extrapolation(), plm_functions::plm_reconstruction(), ppm_functions::ppm_boundary_extrapolation(), ppm_functions::ppm_reconstruction(), pqm_functions::pqm_boundary_extrapolation_v1(), pqm_functions::pqm_reconstruction(), remapping_pcm, remapping_plm, remapping_ppm_h4, remapping_ppm_ih4, remapping_pqm_ih4ih3, and remapping_pqm_ih6ih5.

Referenced by remapping_core_h(), and remapping_core_w().

  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)  :: remapping_scheme !< Remapping scheme
  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_coefficients !< 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
  ! Local variables
  integer :: local_remapping_scheme

  ! Reset polynomial
  ppoly_r_e(:,:) = 0.0
  ppoly_r_s(:,:) = 0.0
  ppoly_r_coefficients(:,:) = 0.0
  imethod = -999

  local_remapping_scheme = 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_coefficients)
      imethod = integration_pcm
    case ( remapping_plm )
      call plm_reconstruction( n0, h0, u0, ppoly_r_e, ppoly_r_coefficients )
      if ( boundary_extrapolation) then
        call plm_boundary_extrapolation( n0, h0, u0, ppoly_r_e, ppoly_r_coefficients)
      end if
      imethod = integration_plm
    case ( remapping_ppm_h4 )
      call edge_values_explicit_h4( n0, h0, u0, ppoly_r_e )
      call ppm_reconstruction( n0, h0, u0, ppoly_r_e, ppoly_r_coefficients )
      if ( boundary_extrapolation) then
        call ppm_boundary_extrapolation( n0, h0, u0, ppoly_r_e, ppoly_r_coefficients )
      end if
      imethod = integration_ppm
    case ( remapping_ppm_ih4 )
      call edge_values_implicit_h4( n0, h0, u0, ppoly_r_e )
      call ppm_reconstruction( n0, h0, u0, ppoly_r_e, ppoly_r_coefficients )
      if ( boundary_extrapolation) then
        call ppm_boundary_extrapolation( n0, h0, u0, ppoly_r_e, ppoly_r_coefficients )
      end if
      imethod = integration_ppm
    case ( remapping_pqm_ih4ih3 )
      call edge_values_implicit_h4( n0, h0, u0, ppoly_r_e )
      call edge_slopes_implicit_h3( n0, h0, u0, ppoly_r_s )
      call pqm_reconstruction( n0, h0, u0, ppoly_r_e, ppoly_r_s, ppoly_r_coefficients )
      if ( boundary_extrapolation) then
        call pqm_boundary_extrapolation_v1( n0, h0, u0, ppoly_r_e, ppoly_r_s, ppoly_r_coefficients )
      end if
      imethod = integration_pqm
    case ( remapping_pqm_ih6ih5 )
      call edge_values_implicit_h6( n0, h0, u0, ppoly_r_e )
      call edge_slopes_implicit_h5( n0, h0, u0, ppoly_r_s )
      call pqm_reconstruction( n0, h0, u0, ppoly_r_e, ppoly_r_s, ppoly_r_coefficients )
      if ( boundary_extrapolation) then
        call pqm_boundary_extrapolation_v1( n0, h0, u0, ppoly_r_e, ppoly_r_s, ppoly_r_coefficients )
      end if
      imethod = integration_pqm
    case default
      call mom_error( fatal, 'MOM_remapping, build_reconstructions_1d: '//&
           'The selected remapping method is invalid' )
  end select

Here is the call graph for this function:

Here is the caller graph for this function:

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 114 of file MOM_remapping.F90.

Referenced by remapping_unit_tests().

  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)
  end do

Here is the caller graph for this function:

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_coefficients, 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_coefficients	Coefficients of polynomial
[out]	ppoly_r_e	Edge value of polynomial
[out]	ppoly_r_s	Edge slope of polynomial

Definition at line 383 of file MOM_remapping.F90.

References mom_error_handler::mom_error().

Referenced by remapping_core_h(), and remapping_core_w().

  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_coefficients !< 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_coefficients(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_coefficients(n,:)
      enddo
      call mom_error(fatal, 'MOM_remapping, check_reconstructions_1d: '// &
                   'Edge values or polynomial coefficients were inconsistent!')
    endif
  enddo

Here is the call graph for this function:

Here is the caller graph for this function:

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 1771 of file MOM_remapping.F90.

Referenced by remapping_unit_tests().

  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(*,'("i=",20i10)') (i,i=1,n+1)
  write(*,'("x=",20es10.2)') (x(i),i=1,n+1)
  write(*,'("i=",5x,20i10)') (i,i=1,n)
  write(*,'("h=",5x,20es10.2)') (h(i),i=1,n)
  write(*,'("u=",5x,20es10.2)') (u(i),i=1,n)

Here is the caller graph for this function:

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 1440 of file MOM_remapping.F90.

Referenced by mom_ale::ale_remap_scalar(), and remapping_unit_tests().

  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

Here is the caller graph for this function:

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 1520 of file MOM_remapping.F90.

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

  cs%degree = 0

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
	)

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

Definition at line 1465 of file MOM_remapping.F90.

References remapping_set_param().

Referenced by mom_ale_sponge::initialize_ale_sponge(), and remapping_unit_tests().

  ! 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

  ! Note that remapping_scheme is mandatory fir 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)

Here is the call graph for this function:

Here is the caller graph for this function:

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_coefficients, 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  )

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_coefficients	Coefficients of polynomial
[in]	method	Remapping scheme to use
[in]	xr	Left/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.

Definition at line 1148 of file MOM_remapping.F90.

References h_neglect, integration_pcm, integration_plm, integration_ppm, integration_pqm, and mom_error_handler::mom_error().

Referenced by remapbydeltaz(), and remapbyprojection().

  integer, intent(in)  :: n0       !< Number of cells in source grid
  real,    intent(in)  :: h0(:)    !< Source grid sizes (size n0)
  real,    intent(in)  :: u0(:)    !< Source cell averages
  real,    intent(in)  :: ppoly0_e(:,:)            !< Edge value of polynomial
  real,    intent(in)  :: ppoly0_coefficients(:,:) !< Coefficients of polynomial
  integer, intent(in)  :: method   !< Remapping scheme to use
  real,    intent(in)  :: xl, xr   !< Left/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.
  ! 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, parameter :: r_3 = 1.0/3.0 ! Used in evaluation of integrated polynomials

  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) + h_neglect ) - x0jll / ( h0(jl) + h_neglect )

      select case ( method )
        case ( integration_pcm )
          uave =         ppoly0_coefficients(jl,1)
        case ( integration_plm )
          uave =         ppoly0_coefficients(jl,1)   &
               + xi0 *   ppoly0_coefficients(jl,2)
        case ( integration_ppm )
          uave =         ppoly0_coefficients(jl,1)   &
               + xi0 * ( ppoly0_coefficients(jl,2)   &
               + xi0 *   ppoly0_coefficients(jl,3) )
        case ( integration_pqm )
          uave =         ppoly0_coefficients(jl,1)   &
               + xi0 * ( ppoly0_coefficients(jl,2)   &
               + xi0 * ( ppoly0_coefficients(jl,3)   &
               + xi0 * ( ppoly0_coefficients(jl,4)   &
               + xi0 *   ppoly0_coefficients(jl,5) ) ) )
        case default
          call mom_error( fatal,'The selected integration method is invalid' )
      end select

    end if ! 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) + h_neglect ) ) )
      xi1 = max( 0., min( 1., ( xr - x0jll ) / ( h0(jl) + h_neglect ) ) )
#else
      xi0 = xl / h0(jl) - x0jll / ( h0(jl) + h_neglect )
      xi1 = xr / h0(jl) - x0jll / ( h0(jl) + h_neglect )
#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_coefficients(jl,1)
        case ( integration_plm )
          q = ( xr - xl ) * (                                 &
              ppoly0_coefficients(jl,1)                       &
            + ppoly0_coefficients(jl,2) * 0.5 * ( xi1 + xi0 ) )
        case ( integration_ppm )
          q = ( xr - xl ) * (                                   &
                ppoly0_coefficients(jl,1)                       &
            + ( ppoly0_coefficients(jl,2) * 0.5 * ( xi1 + xi0 ) &
            +   ppoly0_coefficients(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_coefficients(jl,1)                                                         &
            + ( ppoly0_coefficients(jl,2) * 0.5 * ( xi1 + xi0 )                                   &
            + ( ppoly0_coefficients(jl,3) * r_3 * ( x02px12 + xi0*xi1 )                           &
            +   ppoly0_coefficients(jl,4) * 0.25* ( x02px12 * x0px1 )                             &
            +   ppoly0_coefficients(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) + h_neglect ) ) )
#else
      xi0 = (xl - x0jll) / ( h0(jl) + h_neglect )
#endif
      xi1 = 1.0

      hact = h0(jl) * ( xi1 - xi0 )

      select case ( method )
        case ( integration_pcm )
          q = q + ( x0jlr - xl ) * ppoly0_coefficients(jl,1)
        case ( integration_plm )
          q = q + ( x0jlr - xl ) * (                          &
              ppoly0_coefficients(jl,1)                       &
            + ppoly0_coefficients(jl,2) * 0.5 * ( xi1 + xi0 ) )
        case ( integration_ppm )
          q = q + ( x0jlr - xl ) * (                            &
                ppoly0_coefficients(jl,1)                       &
            + ( ppoly0_coefficients(jl,2) * 0.5 * ( xi1 + xi0 ) &
            +   ppoly0_coefficients(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_coefficients(jl,1)                                                         &
            + ( ppoly0_coefficients(jl,2) * 0.5 * ( xi1 + xi0 )                                   &
            + ( ppoly0_coefficients(jl,3) * r_3 * ( x02px12 + xi0*xi1 )                           &
            +   ppoly0_coefficients(jl,4) * 0.25* ( x02px12 * x0px1 )                             &
            +   ppoly0_coefficients(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)
        end do
      end if

      ! 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) + h_neglect ) ) )
#else
      xi1 = (xr - x0jrl) / ( h0(jr) + h_neglect )
#endif

      hact = hact + h0(jr) * ( xi1 - xi0 )

      select case ( method )
        case ( integration_pcm )
          q = q + ( xr - x0jrl ) * ppoly0_coefficients(jr,1)
        case ( integration_plm )
          q = q + ( xr - x0jrl ) * (                          &
              ppoly0_coefficients(jr,1)                       &
            + ppoly0_coefficients(jr,2) * 0.5 * ( xi1 + xi0 ) )
        case ( integration_ppm )
          q = q + ( xr - x0jrl ) * (                            &
                ppoly0_coefficients(jr,1)                       &
            + ( ppoly0_coefficients(jr,2) * 0.5 * ( xi1 + xi0 ) &
            +   ppoly0_coefficients(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_coefficients(jr,1)                                                         &
            + ( ppoly0_coefficients(jr,2) * 0.5 * ( xi1 + xi0 )                                   &
            + ( ppoly0_coefficients(jr,3) * r_3 * ( x02px12 + xi0*xi1 )                           &
            +   ppoly0_coefficients(jr,4) * 0.25* ( x02px12 * x0px1 )                             &
            +   ppoly0_coefficients(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

    end if ! 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_coefficients(jl,1)
else
    uave = q / hact
endif
#else
    uave = q / hc
#endif

  end if ! end if clause to check if cell is vanished

Here is the call graph for this function:

Here is the caller graph for this function:

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 136 of file MOM_remapping.F90.

References mom_error_handler::mom_error().

  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

Here is the call graph for this function:

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)  ppoly_E, 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]	ppoly_e	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 961 of file MOM_remapping.F90.

Referenced by remap_via_sub_cells(), remapping_core_h(), and remapping_core_w().

  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)  :: ppoly_e !< 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( ppoly_e(1,1), ppoly_e(1,2) )
  u0max = max( ppoly_e(1,1), ppoly_e(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, ppoly_e(k,1), ppoly_e(k,2) )
    u0max = max( u0max, ppoly_e(k,1), ppoly_e(k,2) )
  enddo

Here is the caller graph for this function:

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 995 of file MOM_remapping.F90.

Referenced by remap_via_sub_cells(), remapping_core_h(), and remapping_core_w().

  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

Here is the caller graph for this function:

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_coefficients, 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.

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_coefficients	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 457 of file MOM_remapping.F90.

References average_value_ppoly(), measure_input_bounds(), measure_output_bounds(), mom_error_handler::mom_error(), and old_algorithm.

Referenced by remapping_core_h(), remapping_core_w(), and remapping_unit_tests().

  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_coefficients(:,:) !< 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_coefficients, 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_coefficients(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_coefficients(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)

Here is the call graph for this function:

Here is the caller graph for this function:

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_coefficients, integer, intent(in)  n1, real, dimension(:), intent(in)  dx1, integer  method, real, dimension(:), intent(out)  u1, real, dimension(:), intent(out), optional  h1  )

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 widths (size n0)
[in]	u0	Source cell averages (size n0)
[in]	ppoly0_e	Edge value of polynomial
[in]	ppoly0_coefficients	Coefficients of polynomial
[in]	n1	Number of cells in target grid
[in]	dx1	Target grid edge positions (size n1+1)
	method	Remapping scheme to use
[out]	u1	Target cell averages (size n1)
[out]	h1	Target grid widths (size n1)

Definition at line 1073 of file MOM_remapping.F90.

References integraterecononinterval().

Referenced by remapping_unit_tests().

  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_coefficients(:,:) !< Coefficients of polynomial
  integer,        intent(in)  :: n1     !< Number of cells in target grid
  real,           intent(in)  :: dx1(:) !< Target grid edge positions (size n1+1)
  integer                     :: method !< Remapping scheme to use
  real,           intent(out) :: u1(:)  !< Target cell averages (size n1)
  real, optional, intent(out) :: h1(:)  !< Target grid widths (size n1)
  ! 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_coefficients, 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

  end do ! end iTarget loop on target grid cells

Here is the call graph for this function:

Here is the caller graph for this function:

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_coefficients, integer, intent(in)  n1, real, dimension(:), intent(in)  h1, integer, intent(in)  method, real, dimension(:), intent(out)  u1  )

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_coefficients	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)

Definition at line 1029 of file MOM_remapping.F90.

References integraterecononinterval().

Referenced by remapping_unit_tests().

  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_coefficients(:,:) !< 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)
  ! 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_coefficients, method, &
                                   xl, xr, h1(itarget), u1(itarget), jstart, xstart )

  end do ! end iTarget loop on target grid cells

Here is the call graph for this function:

Here is the caller graph for this function:

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
	)

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

Definition at line 163 of file MOM_remapping.F90.

References build_reconstructions_1d(), check_reconstructions_1d(), measure_input_bounds(), measure_output_bounds(), mom_error_handler::mom_error(), and remap_via_sub_cells().

Referenced by coord_rho::build_rho_column(), mom_state_initialization::cut_off_column_top(), mom_diag_remap::diag_remap_do_remap(), and mom_wave_speed::wave_speed().

  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
  ! 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_coefficients !Coefficients of polynomial
  integer :: k
  real :: eps, h0tot, h0err, h1tot, h1err, u0tot, u0err, u0min, u0max, u1tot, u1err, u1min, u1max, uh_err

  call build_reconstructions_1d( n0, h0, u0, cs%remapping_scheme, cs%degree, cs%boundary_extrapolation, &
                                   ppoly_r_coefficients, ppoly_r_e, ppoly_r_s, imethod )

  if (cs%check_reconstruction) call check_reconstructions_1d(n0, h0, u0, cs%degree, &
                                   cs%boundary_extrapolation, ppoly_r_coefficients, ppoly_r_e, ppoly_r_s)


  call remap_via_sub_cells( n0, h0, u0, ppoly_r_e, ppoly_r_coefficients, 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_coefficients(k,:)
      enddo
      call mom_error( fatal, 'MOM_remapping, remapping_core_h: '//&
             'Remapping result is inconsistent!' )
    endif
    endif ! method<5
  endif

Here is the call graph for this function:

Here is the caller graph for this function:

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
	)

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

Definition at line 229 of file MOM_remapping.F90.

References build_reconstructions_1d(), check_reconstructions_1d(), measure_input_bounds(), measure_output_bounds(), mom_error_handler::mom_error(), and remap_via_sub_cells().

Referenced by remapping_unit_tests().

  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
  ! 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_coefficients !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

  call build_reconstructions_1d( n0, h0, u0, cs%remapping_scheme, cs%degree, cs%boundary_extrapolation, &
                                   ppoly_r_coefficients, ppoly_r_e, ppoly_r_s, imethod )

  if (cs%check_reconstruction) call check_reconstructions_1d(n0, h0, u0, cs%degree, &
                                   cs%boundary_extrapolation, ppoly_r_coefficients, 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_coefficients, n1, h1, imethod, &
                            cs%force_bounds_in_subcell,u1, uh_err )
! call remapByDeltaZ( n0, h0, u0, ppoly_r_E, ppoly_r_coefficients, n1, dx, iMethod, u1 )
! call remapByProjection( n0, h0, u0, CS%ppoly_r, n1, h1, iMethod, u1 )

  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_coefficients(k,:)
      enddo
      call mom_error( fatal, 'MOM_remapping, remapping_core_w: '//&
             'Remapping result is inconsistent!' )
    endif
    endif ! method<5
  endif

Here is the call graph for this function:

Here is the caller graph for this function:

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
	)

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

Definition at line 88 of file MOM_remapping.F90.

References setreconstructiontype().

Referenced by initialize_remapping().

  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

  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

Here is the call graph for this function:

Here is the caller graph for this function:

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 1530 of file MOM_remapping.F90.

References buildgridfromh(), dumpgrid(), dzfromh1h2(), initialize_remapping(), integration_plm, integration_ppm, pcm_functions::pcm_reconstruction(), plm_functions::plm_reconstruction(), ppm_functions::ppm_boundary_extrapolation(), ppm_functions::ppm_reconstruction(), remap_via_sub_cells(), remapbydeltaz(), remapbyprojection(), remapping_core_w(), and test_answer().

Referenced by mom_unit_tests::unit_tests().

  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_coefficients
  integer :: i
  real :: err
  logical :: thistest, v

  v = verbose

  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')
  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 )
  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_coefficients(n0,cs%degree+1))

  ppoly0_e(:,:) = 0.0
  ppoly0_s(:,:) = 0.0
  ppoly0_coefficients(:,:) = 0.0

  call edge_values_explicit_h4( n0, h0, u0, ppoly0_e )
  call ppm_reconstruction( n0, h0, u0, ppoly0_e, ppoly0_coefficients )
  call ppm_boundary_extrapolation( n0, h0, u0, ppoly0_e, ppoly0_coefficients )
  u1(:) = 0.
  call remapbyprojection( n0, h0, u0, ppoly0_e, ppoly0_coefficients, &
                          n1, h1, integration_ppm, u1 )
  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_coefficients, &
                      n1, x1-x0(1:n1+1), &
                      integration_ppm, u1, hn1 )
  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_coefficients, &
                      n2, dx2, &
                      integration_ppm, u2, hn2 )
  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_coefficients, &
                            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_coefficients, &
                            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_coefficients, &
                            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_coefficients)
  allocate(ppoly0_coefficients(5,6))
  allocate(ppoly0_e(5,2))
  allocate(ppoly0_s(5,2))

  call pcm_reconstruction(3, (/1.,2.,4./), ppoly0_e(1:3,:), ppoly0_coefficients(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_coefficients(:,1), (/1.,2.,4./), 'PCM: P0')

  call plm_reconstruction(3, (/1.,1.,1./), (/1.,3.,5./), ppoly0_e(1:3,:), ppoly0_coefficients(1:3,:) )
  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_coefficients(:,1), (/1.,2.,5./), 'Unlim PLM: P0')
  remapping_unit_tests = remapping_unit_tests .or. test_answer(v, 3, ppoly0_coefficients(:,2), (/0.,2.,0./), 'Unlim PLM: P1')

  call plm_reconstruction(3, (/1.,1.,1./), (/1.,2.,7./), ppoly0_e(1:3,:), ppoly0_coefficients(1:3,:) )
  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_coefficients(:,1), (/1.,1.,7./), 'Left lim PLM: P0')
  remapping_unit_tests = remapping_unit_tests .or. test_answer(v, 3, ppoly0_coefficients(:,2), (/0.,2.,0./), 'Left lim PLM: P1')

  call plm_reconstruction(3, (/1.,1.,1./), (/1.,6.,7./), ppoly0_e(1:3,:), ppoly0_coefficients(1:3,:) )
  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_coefficients(:,1), (/1.,5.,7./), 'Right lim PLM: P0')
  remapping_unit_tests = remapping_unit_tests .or. test_answer(v, 3, ppoly0_coefficients(:,2), (/0.,2.,0./), 'Right lim PLM: P1')

  call plm_reconstruction(3, (/1.,2.,3./), (/1.,4.,9./), ppoly0_e(1:3,:), ppoly0_coefficients(1:3,:) )
  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_coefficients(:,1), (/1.,2.,9./), 'Non-uniform line PLM: P0')
  remapping_unit_tests = remapping_unit_tests .or. test_answer(v, 3, ppoly0_coefficients(:,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 )
  thistest = test_answer(v, 5, ppoly0_e(:,1), (/0.,2.,4.,6.,8./), 'Line H4: left edges') ! Currently fails due to roundoff
  thistest = test_answer(v, 5, ppoly0_e(:,2), (/2.,4.,6.,8.,10./), 'Line H4: right edges') ! Currently fails due to roundoff
  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_coefficients(1:5,:) )
  remapping_unit_tests = remapping_unit_tests .or. test_answer(v, 5, ppoly0_coefficients(:,1), (/1.,2.,4.,6.,9./), 'Line PPM: P0')
  remapping_unit_tests = remapping_unit_tests .or. test_answer(v, 5, ppoly0_coefficients(:,2), (/0.,2.,2.,2.,0./), 'Line PPM: P1')
  remapping_unit_tests = remapping_unit_tests .or. test_answer(v, 5, ppoly0_coefficients(:,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 )
  thistest = test_answer(v, 5, ppoly0_e(:,1), (/3.,0.,3.,12.,27./), 'Parabola H4: left edges') ! Currently fails due to roundoff
  thistest = test_answer(v, 5, ppoly0_e(:,2), (/0.,3.,12.,27.,48./), 'Parabola H4: right edges') ! Currently fails due to roundoff
  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_coefficients(1:5,:) )
  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_coefficients(:,1), (/0.,0.,3.,12.,37./), 'Parabola PPM: P0')
  remapping_unit_tests = remapping_unit_tests .or. test_answer(v, 5, ppoly0_coefficients(:,2), (/0.,0.,6.,12.,0./), 'Parabola PPM: P1')
  remapping_unit_tests = remapping_unit_tests .or. test_answer(v, 5, ppoly0_coefficients(:,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_coefficients(1:5,:) )
  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_coefficients(:,1), (/0.,3.,6.,16.,15./), 'Limits PPM: P0')
  remapping_unit_tests = remapping_unit_tests .or. test_answer(v, 5, ppoly0_coefficients(:,2), (/0.,6.,0.,0.,0./), 'Limits PPM: P1')
  remapping_unit_tests = remapping_unit_tests .or. test_answer(v, 5, ppoly0_coefficients(:,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_coefficients(1:4,:) )
  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_coefficients(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_coefficients)

  if (.not. remapping_unit_tests) write(*,*) 'Pass'

Here is the call graph for this function:

Here is the caller graph for this function:

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 1485 of file MOM_remapping.F90.

References mom_error_handler::mom_error(), remapping_pcm, remapping_plm, remapping_ppm_h4, remapping_ppm_ih4, remapping_pqm_ih4ih3, remapping_pqm_ih6ih5, and mom_string_functions::uppercase().

Referenced by remapping_set_param().

  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

Here is the call graph for this function:

Here is the caller graph for this function:

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  )

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

Definition at line 1744 of file MOM_remapping.F90.

Referenced by remapping_unit_tests().

  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
  ! Local variables
  integer :: k

  test_answer = .false.
  do k = 1, n
    if (u(k) /= u_true(k)) test_answer = .true.
  enddo
  if (test_answer .or. verbose) then
    write(*,'(a4,2a24,x,a)') 'k','Calculated value','Correct value',label
    do k = 1, n
      if (u(k) /= u_true(k)) then
        write(*,'(i4,1p2e24.16,a,1pe24.16,a)') k,u(k),u_true(k),' err=',u(k)-u_true(k),' < wrong'
      else
        write(*,'(i4,1p2e24.16)') k,u(k),u_true(k)
      endif
    enddo
  endif

Here is the caller graph for this function:

Variable Documentation

h_neglect

real, parameter mom_remapping::h_neglect = 1.E-30

private

A dimensional (H units) number that can be added to thicknesses in a denominator without changing the numerical result, except where a division by zero would otherwise occur.

Definition at line 74 of file MOM_remapping.F90.

Referenced by integraterecononinterval().

real, parameter :: h_neglect = 1.e-30 !< A dimensional (H units) number that can be

integration_pcm

integer, parameter mom_remapping::integration_pcm = 0

private

Piecewise Constant Method.

Definition at line 52 of file MOM_remapping.F90.

Referenced by average_value_ppoly(), build_reconstructions_1d(), and integraterecononinterval().

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

integration_plm

integer, parameter mom_remapping::integration_plm = 1

private

Piecewise Linear Method.

Definition at line 53 of file MOM_remapping.F90.

Referenced by average_value_ppoly(), build_reconstructions_1d(), integraterecononinterval(), and remapping_unit_tests().

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

integration_ppm

integer, parameter mom_remapping::integration_ppm = 3

private

Piecewise Parabolic Method.

Definition at line 54 of file MOM_remapping.F90.

Referenced by average_value_ppoly(), build_reconstructions_1d(), integraterecononinterval(), and remapping_unit_tests().

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

integration_pqm

integer, parameter mom_remapping::integration_pqm = 5

private

Piecewise Quartic Method.

Definition at line 55 of file MOM_remapping.F90.

Referenced by average_value_ppoly(), build_reconstructions_1d(), and integraterecononinterval().

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

mdl

character(len=40) mom_remapping::mdl = "MOM_remapping"

private

This module's name.

Definition at line 57 of file MOM_remapping.F90.

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

old_algorithm

logical, parameter mom_remapping::old_algorithm = .false.

private

Use the old "broken" algorithm. This is a temporary measure to assist debugging until we delete the old algorithm.

Definition at line 79 of file MOM_remapping.F90.

Referenced by remap_via_sub_cells().

logical, parameter :: old_algorithm = .false. !< Use the old "broken" algorithm.

remapping_pcm

integer, parameter mom_remapping::remapping_pcm = 0

O(h^1) remapping scheme.

Definition at line 45 of file MOM_remapping.F90.

Referenced by build_reconstructions_1d(), and setreconstructiontype().

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

remapping_plm

integer, parameter mom_remapping::remapping_plm = 1

private

O(h^2) remapping scheme.

Definition at line 46 of file MOM_remapping.F90.

Referenced by build_reconstructions_1d(), and setreconstructiontype().

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

remapping_ppm_h4

integer, parameter mom_remapping::remapping_ppm_h4 = 2

private

O(h^3) remapping scheme.

Definition at line 47 of file MOM_remapping.F90.

Referenced by build_reconstructions_1d(), and setreconstructiontype().

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

remapping_ppm_ih4

integer, parameter mom_remapping::remapping_ppm_ih4 = 3

private

O(h^3) remapping scheme.

Definition at line 48 of file MOM_remapping.F90.

Referenced by build_reconstructions_1d(), and setreconstructiontype().

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

remapping_pqm_ih4ih3

integer, parameter mom_remapping::remapping_pqm_ih4ih3 = 4

private

O(h^4) remapping scheme.

Definition at line 49 of file MOM_remapping.F90.

Referenced by build_reconstructions_1d(), and setreconstructiontype().

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

remapping_pqm_ih6ih5

integer, parameter mom_remapping::remapping_pqm_ih6ih5 = 5

private

O(h^5) remapping scheme.

Definition at line 50 of file MOM_remapping.F90.

Referenced by build_reconstructions_1d(), and setreconstructiontype().

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

remappingdefaultscheme

character(len=3), public mom_remapping::remappingdefaultscheme = "PLM"

Default remapping method.

Definition at line 67 of file MOM_remapping.F90.

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

remappingschemesdoc

character(len=256), public mom_remapping::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.

Definition at line 60 of file MOM_remapping.F90.

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"