cubic_polynomials.f90

Go to the documentation of this file.

subroutine cubic_polynomials(i, j, r0, s0, r, s, h, hr, hs, &
 hrs, hrr, hss)
  
!------------------------------------------------------------------
! subroutine to calculate the value of the cubic polynomials and
! the derivatives of the corner marked by r0,s0 at the position r,s
!------------------------------------------------------------------

  use itm_types

  implicit none
  
  integer(itm_i4), intent(in) :: i, j
  real(r8), intent(in) :: r0, s0, r, s
  
  real(r8), intent(out) :: h, hr, hs, hrs, hrr, hss
  
  real(r8) :: hi, hri, hrri
  real(r8) :: hj, hsj, hssj
  
  if (i == 0) then
    hi = -(r + r0)**2 * (r * r0 - 2._r8) / 4._r8
    hri = -(r + r0) * (r * r0 - 2._r8) / 2._r8 - r0 * (r + r0)**2 &
     / 4._r8
    hrri = -1.5_r8 * r * r0
  else
    hi = r0 * (r + r0)**2 * (r * r0 - 1._r8) / 4._r8
    hri = r0 * (r + r0) * (r * r0 - 1._r8) / 2._r8 + (r + r0)**2 &
     / 4._r8
    hrri = 1.5_r8 * r + 0.5_r8 * r0
  end if
  if (j == 0) then
    hj = -(s + s0)**2 * (s * s0 - 2._r8) / 4._r8
    hsj = -(s + s0) * (s * s0 - 2._r8) / 2._r8 - s0 * (s + s0)**2 &
     / 4._r8
    hssj = -1.5_r8 * s * s0
  else
    hj = s0 * (s + s0)**2 * (s * s0 - 1._r8) / 4._r8
    hsj = s0 * (s + s0) * (s * s0 - 1._r8) / 2._r8 + (s + s0)**2 &
     / 4._r8
    hssj = 1.5_r8 * s + 0.5_r8 * s0
  end if
  
  h = hi * hj
  hr = hri * hj
  hs = hi * hsj
  hrs = hri * hsj
  hrr = hrri * hj
  hss = hi * hssj
  
  return
end subroutine cubic_polynomials