mod_interpolation.f90

Go to the documentation of this file.

module mod_interpolation

  use itm_types

  implicit none

  contains

!-----------------------------------------------------------------------------
  subroutine interpolation(type_interp, xn1, xn2, xn3, xn4, r, s, x, xr, xs, &
   xrs, xrr, xss, yn1, yn2, yn3, yn4, pn1, pn2, pn3, pn4, yr, ys, ps)
!-----------------------------------------------------------------------------
! subroutine calculates the interpolated value of the function x given
! by xi(1..4) at the four nodes using bi-cubic hermite elements
!
! input:     
!   type_interp
!     0 : function value
!     1 : funct. value + 1st derivatives
!     2 : funct. value + 1st and 2nd derivatives
!     3 : funct. value + 1st derivatives + yr, ys, ps
!   xn1, xn2, xn3, xn4
!   r, s
!   optional: yn1, yn2, yn3, yn4, pn1, pn2, pn3, pn4 (type_interp = 3)
!
! output:
!   x
!   optional: xr, xs, xs, xrs, xrr, xss, yr, ys, ps
!-----------------------------------------------------------------------------

    real(r8), dimension(4), intent(in) :: xn1, xn2, xn3, xn4
    real(r8), dimension(4), intent(in), optional :: yn1, yn2, yn3, yn4
    real(r8), dimension(4), intent(in), optional :: pn1, pn2, pn3, pn4
    real(r8), intent(in) :: r, s
    integer(itm_i4), intent(in) :: type_interp

    real(r8), intent(out) :: x
    real(r8), intent(out), optional :: xr, xs
    real(r8), intent(out), optional :: xrs, xrr, xss
    real(r8), intent(out), optional :: yr, ys, ps

    real(r8) :: hi0m, hri0m, hrri0m, hi1m, hri1m, hrri1m
    real(r8) :: hj0m, hsj0m, hssj0m, hj1m, hsj1m, hssj1m
    real(r8) :: hi0p, hri0p, hrri0p, hi1p, hri1p, hrri1p
    real(r8) :: hj0p, hsj0p, hssj0p, hj1p, hsj1p, hssj1p

    hi0m = -(r - 1._r8)**2 * (-r - 2._r8) * 0.25_r8
    hi1m = -(r - 1._r8)**2 * (-r - 1._r8) * 0.25_r8

    hj0m = -(s - 1._r8)**2 * (-s - 2._r8) * 0.25_r8
    hj1m = -(s - 1._r8)**2 * (-s - 1._r8) * 0.25_r8

    hi0p = -(r + 1._r8)**2 * (r - 2._r8) * 0.25_r8
    hi1p = (r + 1._r8)**2 * (r - 1._r8) * 0.25_r8

    hj0p = -(s + 1._r8)**2 * (s - 2._r8) * 0.25_r8
    hj1p = (s + 1._r8)**2 * (s - 1._r8) * 0.25_r8

    x = hi0m * hj0m * xn1(1) + hi1m * hj0m * xn1(2) &
     + hi0m * hj1m * xn1(3) + hi1m * hj1m * xn1(4) & 
     + hi0m * hj0p * xn2(1) + hi1m * hj0p * xn2(2) &
     + hi0m * hj1p * xn2(3) + hi1m * hj1p * xn2(4) & 
     + hi0p * hj0m * xn4(1) + hi1p * hj0m * xn4(2) &
     + hi0p * hj1m * xn4(3) + hi1p * hj1m * xn4(4) & 
     + hi0p * hj0p * xn3(1) + hi1p * hj0p * xn3(2) &
     + hi0p * hj1p * xn3(3) + hi1p * hj1p * xn3(4)

    if (type_interp == 0) return

    hri0m = -(r - 1._r8) * (-r - 2._r8) * 0.5_r8 &
     + (r - 1._r8)**2 * 0.25_r8
    hri1m = -(r - 1._r8) * (-r - 1._r8) * 0.5_r8 &
     + (r - 1._r8)**2 * 0.25_r8

    hsj0m = -(s - 1._r8) * (-s - 2._r8) * 0.5_r8 &
     + (s - 1._r8)**2 * 0.25_r8
    hsj1m = -(s - 1._r8) * (-s - 1._r8) * 0.5_r8 &
     + (s - 1._r8)**2 * 0.25_r8

    hri0p = -(r + 1._r8) * (r - 2._r8) * 0.5_r8 - (r + 1._r8)**2 &
     * 0.25_r8
    hri1p = (r + 1._r8) * (r - 1._r8) * 0.5_r8 + (r + 1._r8)**2 &
     * 0.25_r8

    hsj0p = -(s + 1._r8) * (s - 2._r8) * 0.5_r8 - (s + 1._r8)**2 &
     * 0.25_r8
    hsj1p = (s + 1._r8) * (s - 1._r8) * 0.5_r8 + (s + 1._r8)**2 &
     * 0.25_r8

    xr = hri0m * hj0m * xn1(1) + hri1m * hj0m * xn1(2) &
     + hri0m * hj1m * xn1(3) + hri1m * hj1m * xn1(4) &    
     + hri0m * hj0p * xn2(1) + hri1m * hj0p * xn2(2) &
     + hri0m * hj1p * xn2(3) + hri1m * hj1p * xn2(4) & 
     + hri0p * hj0m * xn4(1) + hri1p * hj0m * xn4(2) &
     + hri0p * hj1m * xn4(3) + hri1p * hj1m * xn4(4) & 
     + hri0p * hj0p * xn3(1) + hri1p * hj0p * xn3(2) &
     + hri0p * hj1p * xn3(3) + hri1p * hj1p * xn3(4)

    xs = hi0m * hsj0m * xn1(1) + hi1m * hsj0m * xn1(2) &
     + hi0m * hsj1m * xn1(3) + hi1m * hsj1m * xn1(4) & 
     + hi0m * hsj0p * xn2(1) + hi1m * hsj0p * xn2(2) &
     + hi0m * hsj1p * xn2(3) + hi1m * hsj1p * xn2(4) & 
     + hi0p * hsj0m * xn4(1) + hi1p * hsj0m * xn4(2) &
     + hi0p * hsj1m * xn4(3) + hi1p * hsj1m * xn4(4) & 
     + hi0p * hsj0p * xn3(1) + hi1p * hsj0p * xn3(2) &
     + hi0p * hsj1p * xn3(3) + hi1p * hsj1p * xn3(4)

    if (type_interp == 1) return

    if (type_interp == 2) then

      hrri0m = 1.5_r8 * r
      hrri1m = 1.5_r8 * r - 0.5_r8

      hssj0m = 1.5_r8 * s
      hssj1m = 1.5_r8 * s - 0.5_r8

      hrri0p = -1.5_r8 * r
      hrri1p = 1.5_r8 * r + 0.5_r8
  
      hssj0p = -1.5_r8 * s
      hssj1p = 1.5_r8 * s + 0.5_r8

      xrr = hrri0m * hj0m * xn1(1) + hrri1m * hj0m * xn1(2) &
       + hrri0m * hj1m * xn1(3) + hrri1m * hj1m * xn1(4) &
       + hrri0m * hj0p * xn2(1) + hrri1m * hj0p * xn2(2) &
       + hrri0m * hj1p * xn2(3) + hrri1m * hj1p * xn2(4) &
       + hrri0p * hj0m * xn4(1) + hrri1p * hj0m * xn4(2) &
       + hrri0p * hj1m * xn4(3) + hrri1p * hj1m * xn4(4) &
       + hrri0p * hj0p * xn3(1) + hrri1p * hj0p * xn3(2) &
       + hrri0p * hj1p * xn3(3) + hrri1p * hj1p * xn3(4)

      xss = hi0m * hssj0m * xn1(1) + hi1m * hssj0m * xn1(2) &
       + hi0m * hssj1m * xn1(3) + hi1m * hssj1m * xn1(4) &
       + hi0m * hssj0p * xn2(1) + hi1m * hssj0p * xn2(2) &
       + hi0m * hssj1p * xn2(3) + hi1m * hssj1p * xn2(4) &
       + hi0p * hssj0m * xn4(1) + hi1p * hssj0m * xn4(2) &
       + hi0p * hssj1m * xn4(3) + hi1p * hssj1m * xn4(4) &
       + hi0p * hssj0p * xn3(1) + hi1p * hssj0p * xn3(2) &
       + hi0p * hssj1p * xn3(3) + hi1p * hssj1p * xn3(4)

      xrs = hri0m * hsj0m * xn1(1) + hri1m * hsj0m * xn1(2) &
       + hri0m * hsj1m * xn1(3) + hri1m * hsj1m * xn1(4) &
       + hri0m * hsj0p * xn2(1) + hri1m * hsj0p * xn2(2) &
       + hri0m * hsj1p * xn2(3) + hri1m * hsj1p * xn2(4) &
       + hri0p * hsj0m * xn4(1) + hri1p * hsj0m * xn4(2) &
       + hri0p * hsj1m * xn4(3) + hri1p * hsj1m * xn4(4) &
       + hri0p * hsj0p * xn3(1) + hri1p * hsj0p * xn3(2) &
       + hri0p * hsj1p * xn3(3) + hri1p * hsj1p * xn3(4)

    else

      ps = hi0m * hj0m * pn1(1) + hi1m * hj0m * pn1(2) &
       + hi0m * hj1m * pn1(3) + hi1m * hj1m * pn1(4) &
       + hi0m * hj0p * pn2(1) + hi1m * hj0p * pn2(2) &
       + hi0m * hj1p * pn2(3) + hi1m * hj1p * pn2(4) &
       + hi0p * hj0m * pn4(1) + hi1p * hj0m * pn4(2) &
       + hi0p * hj1m * pn4(3) + hi1p * hj1m * pn4(4) &
       + hi0p * hj0p * pn3(1) + hi1p * hj0p * pn3(2) &
       + hi0p * hj1p * pn3(3) + hi1p * hj1p * pn3(4)

      yr = hri0m * hj0m * yn1(1) + hri1m * hj0m * yn1(2) &
       + hri0m * hj1m * yn1(3) + hri1m * hj1m * yn1(4) &
       + hri0m * hj0p * yn2(1) + hri1m * hj0p * yn2(2) &
       + hri0m * hj1p * yn2(3) + hri1m * hj1p * yn2(4) &
       + hri0p * hj0m * yn4(1) + hri1p * hj0m * yn4(2) &
       + hri0p * hj1m * yn4(3) + hri1p * hj1m * yn4(4) &
       + hri0p * hj0p * yn3(1) + hri1p * hj0p * yn3(2) &
       + hri0p * hj1p * yn3(3) + hri1p * hj1p * yn3(4)

      ys = hi0m * hsj0m * yn1(1) + hi1m * hsj0m * yn1(2) &
       + hi0m * hsj1m * yn1(3) + hi1m * hsj1m * yn1(4) &
       + hi0m * hsj0p * yn2(1) + hi1m * hsj0p * yn2(2) &
       + hi0m * hsj1p * yn2(3) + hi1m * hsj1p * yn2(4) &
       + hi0p * hsj0m * yn4(1) + hi1p * hsj0m * yn4(2) &
       + hi0p * hsj1m * yn4(3) + hi1p * hsj1m * yn4(4) &
       + hi0p * hsj0p * yn3(1) + hi1p * hsj0p * yn3(2) &
       + hi0p * hsj1p * yn3(3) + hi1p * hsj1p * yn3(4)
    end if

    return

  end subroutine interpolation

!----------------------------------------------------------------------------
  subroutine transfer_1d_profile(x_in, f_in, x_out, f_out, df_out, &
   spline_type, monotonic)
!----------------------------------------------------------------------------
! This subroutine transfers a 1D profile f_in given on x_in to the profile
! f_out given on x_out by means of spline interpolation.
! df_out holds the derivatives df/dx on output if specified.
! for spline_type == 1 the first derivatives are extrapolated onto the
! boundaries
!----------------------------------------------------------------------------

    use mod_profiles
    use helena_spline

    real(r8), dimension(:) :: x_in, f_in
    real(r8), dimension(:) :: x_out, f_out
    real(r8), dimension(:), intent(inout), optional :: df_out
    integer(itm_i4), optional :: spline_type
    logical, intent(in), optional :: monotonic

    type(spline_coefficients) :: f_spline
    real(r8), dimension(3) :: abltg
    real(r8) :: df1 = 0._r8, dfn = 0._r8
    real(r8) :: df2, df3
    real(r8) :: dfn1, dfn2
    integer(itm_i4) :: s_type
    integer(itm_i4) :: n
    integer(itm_i4) :: i
    logical :: inverted

!-- natural spline as default
    if (.not. present(spline_type)) then
      s_type = 2
    else
      s_type = spline_type
    end if

    if (size(x_in) /= size(f_in) .or. size(x_out) /= size(f_out)) &
     stop 'Fatal error in transfer_1d_profile: size mismatch!'

    if (present(df_out)) then
      if (size(x_out) /= size(df_out)) &
       stop 'Fatal error in transfer_1d_profile: size mismatch!'
    end if

    inverted = .false.
    if (x_in(size(x_in)) < x_in(1)) then
      inverted = .true.    ! inverted abscissa
      x_in = -x_in
      x_out = -x_out
    end if

!-- extrapolate first derivatives for s_type == 1
    if (s_type == 1) then
      n = size(f_in)
      df2 = (f_in(3) - f_in(1)) / (x_in(3) - x_in(1))
      df3 = (f_in(4) - f_in(2)) / (x_in(4) - x_in(2))
      dfn1 = (f_in(n) - f_in(n - 2)) / (x_in(n) - x_in(n - 2))
      dfn2 = (f_in(n - 1) - f_in(n - 3)) / (x_in(n - 1) - x_in(n - 3))

      df1 = (x_in(3) - x_in(1)) / (x_in(3) - x_in(2)) * df2 &
       - (x_in(2) - x_in(1)) / (x_in(3) - x_in(2)) * df3
      dfn = (x_in(n) - x_in(n - 2)) / (x_in(n - 1) - x_in(n - 2)) * dfn1 &
       - (x_in(n) - x_in(n - 1)) / (x_in(n - 1) - x_in(n - 2)) * dfn2
      if(present(monotonic)) then
         write(*,*) 'DPC: dfN, monotonic = ', dfn, monotonic
         if(dfn .lt. 0.0_r8) then
            dfn = 0.0_r8
            write(*,*) 'DPC: dfN changed to ', dfn
         endif
      else
         write(*,*) 'DPC: dfN = ', dfn
      endif
    end if

    call allocate_spline_coefficients(f_spline, size(f_in))
    call spline(size(x_in), x_in, f_in, df1, dfn, s_type, &
     f_spline)

    do i = 1, size(x_out)
      f_out(i) = spwert(size(x_in), x_out(i), f_spline, x_in, abltg, 1)
      if (present(df_out)) df_out(i) = abltg(1)
    end do

!-- enforce monotonicity if monotonic == .true.
    if (present(monotonic)) then
      if (monotonic) then
        do i = 2, size(x_out) - 1
          if ((f_out(i) - f_out(i - 1)) * (f_out(i + 1) - f_out(i)) < 0._r8) &
           then
            write(*,*) 'ERR?: monotonicity fix applied to ', &
                  f_out(i - 1), f_out(i), f_out(i + 1)
            f_out(i) = f_out(i - 1) + (x_out(i) - x_out(i - 1)) &
             * (f_out(i + 1) - f_out(i - 1)) / (x_out(i + 1) - x_out(i - 1))
            write(*,*) 'ERR?: resulting in ', &
                 f_out(i - 1), f_out(i), f_out(i + 1)
            if (present(df_out)) &
             df_out(i) = (f_out(i + 1) - f_out(i - 1)) / (x_out(i + 1) &
             - x_out(i - 1))
          end if
        end do
      end if
    end if

    call deallocate_spline_coefficients(f_spline)

    if (inverted) then
      x_in = -x_in
      x_out = -x_out
      if (present(df_out)) df_out = -df_out
    end if

  end subroutine transfer_1d_profile

end module mod_interpolation