pqm_end_value_h4 Subroutine

private pure subroutine pqm_end_value_h4(dz, u, csys)

One-sided 4th-order polynomial fit of the cell averages u to the four boundary layers dz (thicknesses, must be positive), returning the four coefficients csys of the fit (White & Adcroft 2008, appendix; roundoff-safe closed form). csys(1) is the edge VALUE at the boundary interface and csys(2) is the edge SLOPE there.

Arguments

Type IntentOptional Attributes Name
real(kind=wp), intent(in) :: dz(4)

Thicknesses of the 4 boundary layers, starting at the edge
real(kind=wp), intent(in) :: u(4)

Cell averages of the 4 boundary layers, starting at the edge
real(kind=wp), intent(out) :: csys(4)

Coefficients of the 4th-order fit polynomial in z

Called by

proc~~pqm_end_value_h4~~CalledByGraph proc~pqm_end_value_h4 pqm_end_value_h4 proc~remap_column_pqm remap_column_pqm proc~remap_column_pqm->proc~pqm_end_value_h4 proc~remap_column remap_column proc~remap_column->proc~remap_column_pqm proc~ocean_remap_tracer_column ocean_remap_tracer_column proc~ocean_remap_tracer_column->proc~remap_column proc~ocean_remap_tracer_field ocean_remap_tracer_field proc~ocean_remap_tracer_field->proc~remap_column proc~remap_layer_to_density_impl remap_layer_to_density_impl proc~remap_layer_to_density_impl->proc~remap_column proc~remap_layer_to_vcoord_impl remap_layer_to_vcoord_impl proc~remap_layer_to_vcoord_impl->proc~remap_column proc~remap_tracer_grounded remap_tracer_grounded proc~remap_tracer_grounded->proc~remap_column proc~remap_x_face_grounded remap_x_face_grounded proc~remap_x_face_grounded->proc~remap_column proc~remap_x_face_velocity remap_x_face_velocity proc~remap_x_face_velocity->proc~remap_column proc~remap_y_face_grounded remap_y_face_grounded proc~remap_y_face_grounded->proc~remap_column proc~remap_y_face_velocity remap_y_face_velocity proc~remap_y_face_velocity->proc~remap_column proc~ocean_apply_ale_remap_centres ocean_apply_ale_remap_centres proc~ocean_apply_ale_remap_centres->proc~ocean_remap_tracer_field proc~ocean_apply_ale_remap_faces ocean_apply_ale_remap_faces proc~ocean_apply_ale_remap_faces->proc~remap_x_face_velocity proc~ocean_apply_ale_remap_faces->proc~remap_y_face_velocity proc~ocean_apply_ale_remap_step ocean_apply_ale_remap_step proc~ocean_apply_ale_remap_step->proc~ocean_remap_tracer_field proc~ocean_apply_ale_remap_step->proc~remap_x_face_velocity proc~ocean_apply_ale_remap_step->proc~remap_y_face_velocity proc~ocean_apply_conservative_min_thickness ocean_apply_conservative_min_thickness proc~ocean_apply_conservative_min_thickness->proc~remap_tracer_grounded proc~ocean_apply_conservative_min_thickness->proc~remap_x_face_grounded proc~ocean_apply_conservative_min_thickness->proc~remap_y_face_grounded proc~remap_layer_to_density remap_layer_to_density proc~remap_layer_to_density->proc~remap_layer_to_density_impl proc~remap_layer_to_sigma remap_layer_to_sigma proc~remap_layer_to_sigma->proc~remap_layer_to_vcoord_impl proc~remap_layer_to_z remap_layer_to_z proc~remap_layer_to_z->proc~remap_layer_to_vcoord_impl proc~remap_layer_to_zstar remap_layer_to_zstar proc~remap_layer_to_zstar->proc~remap_layer_to_vcoord_impl proc~ocean_dyn_step_split ocean_dyn_step_split proc~ocean_dyn_step_split->proc~ocean_apply_ale_remap_step proc~run_continuity_chain run_continuity_chain proc~run_continuity_chain->proc~ocean_apply_conservative_min_thickness

Variables

Type Visibility Attributes Name Initial
real(kind=wp), private :: du1
real(kind=wp), private :: du2
real(kind=wp), private :: du3
real(kind=wp), private :: h1
real(kind=wp), private :: h12
real(kind=wp), private :: h123
real(kind=wp), private :: h1234
real(kind=wp), private :: h2
real(kind=wp), private :: h23
real(kind=wp), private :: h234
real(kind=wp), private :: h3
real(kind=wp), private :: h34
real(kind=wp), private :: h4
real(kind=wp), private :: i_denb3
real(kind=wp), private :: i_denom
real(kind=wp), private :: i_h12
real(kind=wp), private :: i_h123
real(kind=wp), private :: i_h1234
real(kind=wp), private :: i_h23
real(kind=wp), private :: i_h234
real(kind=wp), private :: i_h34
real(kind=wp), private :: wt(3,4)

Source Code

   pure subroutine pqm_end_value_h4(dz, u, csys)
      !$acc routine seq
      !! One-sided 4th-order polynomial fit of the cell averages `u` to the
      !! four boundary layers `dz` (thicknesses, must be positive), returning
      !! the four coefficients `csys` of the fit (White & Adcroft 2008,
      !! appendix; roundoff-safe closed form).  `csys(1)` is the edge VALUE at
      !! the boundary interface and `csys(2)` is the edge SLOPE there.
      real(wp), intent(in) :: dz(4)
         !! Thicknesses of the 4 boundary layers, starting at the edge
      real(wp), intent(in) :: u(4)
         !! Cell averages of the 4 boundary layers, starting at the edge
      real(wp), intent(out) :: csys(4)
         !! Coefficients of the 4th-order fit polynomial in z

      real(wp) :: wt(3, 4)
      real(wp) :: h1, h2, h3, h4
      real(wp) :: h12, h23, h34, h123, h234, h1234
      real(wp) :: i_h12, i_h23, i_h34, i_h123, i_h234, i_h1234
      real(wp) :: i_denom, i_denb3
      real(wp) :: du1, du2, du3

      h1 = dz(1)
      h2 = dz(2)
      h3 = dz(3)
      h4 = dz(4)
      ! Bound the thickness ratios so property differences at the level of
      ! roundoff are not amplified to order one.
      if ((h2 + h3) < PQM_MIN_FRAC*h1) h3 = PQM_MIN_FRAC*h1 - h2
      if ((h3 + h4) < PQM_MIN_FRAC*h1) h4 = PQM_MIN_FRAC*h1 - h3

      h12 = h1 + h2
      h23 = h2 + h3
      h34 = h3 + h4
      h123 = h12 + h3
      h234 = h2 + h34
      h1234 = h12 + h34
      ! Three reciprocals from a single division each, for efficiency.
      i_denb3 = 1.0_wp/(h123*h12*h23)
      i_h12 = (h123*h23)*i_denb3
      i_h23 = (h12*h123)*i_denb3
      i_h123 = (h12*h23)*i_denb3
      i_denom = 1.0_wp/(h1234*(h234*h34))
      i_h34 = (h1234*h234)*i_denom
      i_h234 = (h1234*h34)*i_denom
      i_h1234 = (h234*h34)*i_denom

      wt(1, 1) = -h1*(i_h1234 + i_h123 + i_h12)
      wt(2, 1) = h1*h12*(i_h234*i_h1234 + i_h23*(i_h234 + i_h123))
      wt(3, 1) = -h1*h12*h123*i_denom

      wt(1, 2) = 2.0_wp*(i_h12*(1.0_wp + (h1 + h12)*(i_h1234 + i_h123)) + h1*i_h1234*i_h123)
      wt(2, 2) = -2.0_wp*((h1*h12*i_h1234)*(i_h23*(i_h234 + i_h123)) + &
                          (h1 + h12)*(i_h1234*i_h234 + i_h23*(i_h234 + i_h123)))
      wt(3, 2) = 2.0_wp*((h1 + h12)*h123 + h1*h12)*i_denom

      wt(1, 3) = -3.0_wp*i_h12*i_h123*(1.0_wp + i_h1234*((h1 + h12) + h123))
      wt(2, 3) = 3.0_wp*i_h23*(i_h123 + i_h1234*((h1 + h12) + h123)*(i_h123 + i_h234))
      wt(3, 3) = -3.0_wp*((h1 + h12) + h123)*i_denom

      wt(1, 4) = 4.0_wp*i_h1234*i_h123*i_h12
      wt(2, 4) = -4.0_wp*i_h1234*(i_h23*(i_h123 + i_h234))
      wt(3, 4) = 4.0_wp*i_denom

      du1 = u(2) - u(1)
      du2 = u(3) - u(2)
      du3 = u(4) - u(3)
      csys(1) = ((u(1) + (wt(1, 1)*du1)) + (wt(2, 1)*du2)) + (wt(3, 1)*du3)
      csys(2) = ((wt(1, 2)*du1) + (wt(2, 2)*du2)) + (wt(3, 2)*du3)
      csys(3) = ((wt(1, 3)*du1) + (wt(2, 3)*du2)) + (wt(3, 3)*du3)
      csys(4) = ((wt(1, 4)*du1) + (wt(2, 4)*du2)) + (wt(3, 4)*du3)
   end subroutine pqm_end_value_h4