interpolations.f90

Go to the documentation of this file.

MODULE interpolations

  USE l3interps
  USE ids_schemas, ONLY: ids_real

  IMPLICIT NONE

CONTAINS

  SUBROUTINE lininterp( YIN , XIN , NIN , YOUT , XOUT , NOUT )
    REAL(8), INTENT(IN)  :: yin(nin)
    REAL(8), INTENT(IN)  :: xin(nin)
    REAL(8), INTENT(OUT) :: yout(nout)
    REAL(8), INTENT(IN)  :: xout(nout)
    INTEGER, INTENT(IN)  :: nin
    INTEGER, INTENT(IN)  :: nout

    ! Internal
    REAL(8), PARAMETER :: very_small = sqrt(tiny(1.0d0))
    REAL(8) :: frac
    REAL(8) :: direction
    INTEGER :: jin, jout

    ! Check for identical grids:
    IF (nin == nout) THEN
       IF ( maxval(abs(xin-xout)) < very_small) THEN
          yout = yin
          RETURN
       END IF
    END IF

    yout = -9.0d40

    direction = sign(1d0, xin(nin) - xin(1))

    DO jout = 1,nout
       ! Increasing XIN (DIRECTION=+1), then XOUT-XIN(1) < 0 mean that XOUT < XIN(1)
       ! Decreasing XIN (DIRECTION=-1), then -XOUT+XIN(1) < 0 mean that XOUT > XIN(1)
       IF (direction * (xout(jout) - xin(1)) <= 0.0d0) THEN
          yout(jout) = yin(1)
       ELSEIF ( direction * (xout(jout) - xin(nin)) >= 0.0d0 ) THEN
          yout(jout) = yin(nin)
       ELSE
          DO jin = 1,nin
             IF (direction * (xout(jout) - xin(jin+1)) <= 0.0d0) THEN
                frac = (xout(jout) - xin(jin)) / (xin(jin+1) - xin(jin))
                yout(jout) = (1.0d0 - frac) * yin(jin) + frac * yin(jin+1) 
                EXIT
             END IF
          END DO
       END IF
    END DO
  END SUBROUTINE lininterp


  SUBROUTINE linderiv( YIN , XIN , NIN , YOUT , XOUT , NOUT )
    REAL(8), INTENT(IN)  :: yin(nin)
    REAL(8), INTENT(IN)  :: xin(nin)
    REAL(8), INTENT(OUT) :: yout(nout)
    REAL(8), INTENT(IN)  :: xout(nout)
    INTEGER, INTENT(IN)  :: nin
    INTEGER, INTENT(IN)  :: nout

    ! Internal
    REAL(8) :: dx
    REAL(8) :: direction
    INTEGER :: jin, jout

    yout = -9.0d40

    direction = sign(1d0, xin(nin) - xin(1))

    DO jout = 1,nout
       ! Increasing XIN (DIRECTION=+1), then XOUT-XIN(1) < 0 mean that XOUT < XIN(1)
       ! Decreasing XIN (DIRECTION=-1), then -XOUT+XIN(1) < 0 mean that XOUT > XIN(1)
       IF (direction * (xout(jout) - xin(1)) <= 0.0d0) THEN
          yout(jout) = (yin(2) - yin(1)) / (xin(2) - xin(1))
       ELSEIF ( direction * (xout(jout) - xin(nin)) >= 0.0d0 ) THEN
          yout(jout) = (yin(nin) - yin(nin-1)) / (xin(nin) - xin(nin-1))
       ELSE
          DO jin = 1,nin
             IF (direction * (xout(jout) - xin(jin+1)) <= 0.0d0) THEN
                yout(jout) = (yin(jin+1) - yin(jin)) / (xin(jin+1) - xin(jin))
                EXIT
             END IF
          END DO
       END IF
    END DO
  END SUBROUTINE linderiv


  SUBROUTINE linderiv_pointer( YIN , XIN , YOUT , XOUT )
    REAL(8), POINTER, INTENT(IN)  :: yin(:)
    REAL(8), POINTER, INTENT(IN)  :: xin(:)
    REAL(8),          INTENT(OUT) :: yout(:)
    REAL(8),          INTENT(IN)  :: xout(:)

    IF ( (.NOT. ASSOCIATED(xin)) .OR. (.NOT. ASSOCIATED(yin)) ) THEN
       yout(:) = -9d40
       RETURN
    END IF
    IF ( SIZE(xin) /= SIZE(yin) ) THEN
       yout(:) = -9d40
       RETURN
    END IF

    CALL linderiv( yin , xin , SIZE(xin) , yout , xout , SIZE(xout) )

  END SUBROUTINE linderiv_pointer


  SUBROUTINE l3interp_pointer( YIN , XIN , YOUT , XOUT )
    REAL(8), POINTER, INTENT(IN)  ::  yin(:)
    REAL(8), POINTER, INTENT(IN)  ::  xin(:)
    REAL(8),          INTENT(OUT) ::  yout(:)
    REAL(8),          INTENT(IN)  ::  xout(:)

    IF ( (.NOT. ASSOCIATED(xin)) .OR. (.NOT. ASSOCIATED(yin)) ) THEN
       yout(:) = -9d40
       RETURN
    END IF
    IF ( SIZE(xin) /= SIZE(yin) ) THEN
       yout(:) = -9d40
       RETURN
    END IF

    CALL l3interp( yin , xin , SIZE(xin) , yout , xout , SIZE(xout) )

  END SUBROUTINE l3interp_pointer


  SUBROUTINE l3interp_alloc( YIN , XIN , YOUT , XOUT )
    REAL(8), ALLOCATABLE, INTENT(IN)  ::  yin(:)
    REAL(8), ALLOCATABLE, INTENT(IN)  ::  xin(:)
    REAL(8),              INTENT(OUT) ::  yout(:)
    REAL(8),              INTENT(IN)  ::  xout(:)

    IF ( (.NOT. ALLOCATED(xin)) .OR. (.NOT. ALLOCATED(yin)) ) THEN
       yout(:) = -9d40
       RETURN
    END IF
    IF ( SIZE(xin) /= SIZE(yin) ) THEN
       yout(:) = -9d40
       RETURN
    END IF

    CALL l3interp( yin , xin , SIZE(xin) , yout , xout , SIZE(xout) )

  END SUBROUTINE l3interp_alloc


  SUBROUTINE l3deriv_pointer( YIN , XIN , YOUT , XOUT )
    REAL(8), POINTER, INTENT(IN)  ::  yin(:)
    REAL(8), POINTER, INTENT(IN)  ::  xin(:)
    REAL(8),          INTENT(OUT) ::  yout(:)
    REAL(8),          INTENT(IN)  ::  xout(:)

    IF ( (.NOT. ASSOCIATED(xin)) .OR. (.NOT. ASSOCIATED(yin)) ) THEN
       yout(:) = -9d40
       RETURN
    END IF
    IF ( SIZE(xin) /= SIZE(yin) ) THEN
       yout(:) = -9d40
       RETURN
    END IF

    CALL l3deriv( yin , xin , SIZE(xin) , yout , xout , SIZE(xout) )

  END SUBROUTINE l3deriv_pointer


  SUBROUTINE l3deriv_alloc( YIN , XIN , YOUT , XOUT )
    REAL(8), ALLOCATABLE, INTENT(IN)  ::  yin(:)
    REAL(8), ALLOCATABLE, INTENT(IN)  ::  xin(:)
    REAL(8),              INTENT(OUT) ::  yout(:)
    REAL(8),              INTENT(IN)  ::  xout(:)

    IF ( (.NOT. ALLOCATED(xin)) .OR. (.NOT. ALLOCATED(yin)) ) THEN
       yout(:) = -9d40
       RETURN
    END IF
    IF ( SIZE(xin) /= SIZE(yin) ) THEN
       yout(:) = -9d40
       RETURN
    END IF

    CALL l3deriv( yin , xin , SIZE(xin) , yout , xout , SIZE(xout) )

  END SUBROUTINE l3deriv_alloc


  REAL(IDS_REAL) FUNCTION integral(XIN, YIN)
    USE prec_rkind, ONLY: rkind
    USE interpos_module, ONLY: interpos
    REAL(IDS_REAL), INTENT(IN) :: xin(:)
    REAL(IDS_REAL), INTENT(IN) :: yin(:)

    ! Internal
    REAL(RKIND) :: cum_integral(size(xin))
    REAL(RKIND) :: xin_kind(size(xin))
    REAL(RKIND) :: yin_kind(size(xin))
    REAL(RKIND) :: ybc(2) = (/ 0.0_rkind , 0.0_rkind /) ! Both boundary condition are y''=0.
    INTEGER     :: nbc(2) = (/ 0 , 0 /)                 ! Specify y'' at both boundaries using YBC.
    INTEGER     :: nrho

    nrho=SIZE(xin)
    xin_kind = REAL(xin, rkind)
    yin_kind = REAL(yin, rkind)
    CALL interpos(xin_kind, yin_kind, nrho, &
         nrho, xout=xin_kind, youtint=cum_integral, nbc=nbc, ybc=ybc)
    integral = REAL(CUM_INTEGRAL(NRHO), ids_real)

  END FUNCTION integral

END MODULE interpolations