!==================================
! Linear x1-interpolation of f(x1)

SUBROUTINE INTERP1D(nx1_in,x1_in,var_in,nx1_out,x1_out,var_out)

use swim_global_data_mod, only: rspec

implicit none

integer          :: nx1_in, nx1_out, i_x1, i_left, j_x1
real             :: x1_in(nx1_in), var_in(nx1_in)
REAL(KIND=rspec) :: x1_out(nx1_out), var_out(nx1_out)

j_x1=2
DO WHILE( (x1_in(j_x1) .LE. x1_out(1)) .and. (j_x1 .le. nx1_in) )
  j_x1=j_x1+1
ENDDO
i_left=max(1,j_x1-1)

DO j_x1=1,nx1_out
  IF(x1_out(j_x1) .LT. x1_in(1)) THEN
    var_out(j_x1)=real(var_in(1),kind=rspec)
  ELSE IF(x1_out(j_x1) .GT. x1_in(nx1_in)) THEN
    var_out(j_x1)=real(var_in(nx1_in),kind=rspec)
  ELSE
    var_out(j_x1) = real ( var_in(i_left) + &
        (var_in(i_left+1) - var_in(i_left)) * &
        (x1_out(j_x1) - x1_in(i_left)) / &
        (x1_in(i_left+1) - x1_in(i_left)) &
         ,kind=rspec)
    IF (j_x1 .NE. nx1_out) THEN
      i_x1=i_left
      DO WHILE ( (x1_in(i_x1) .LE. x1_out(j_x1+1)) .and. (i_x1 .le. nx1_in) )
        i_x1=i_x1+1
      ENDDO
      i_left=max(i_left,i_x1-1)
    ENDIF
  ENDIF
ENDDO

RETURN
END SUBROUTINE INTERP1D

!=====================================
! Linear x1-interpolation of f(x1,x2)

SUBROUTINE INTERP2D(nx2,nx1_in,x1_in,var_in,nx1_out,x1_out,var_out)

use swim_global_data_mod, only: rspec

implicit none

integer          :: nx2, nx1_in,nx1_out, j_x2
real             :: x1_in(nx1_in), var_in(nx1_in,nx2)
REAL(KIND=rspec) :: x1_out(nx1_out), var_out(nx1_out,nx2)

DO j_x2=1,nx2
  CALL INTERP1D(nx1_in,x1_in,var_in(1:nx1_in,j_x2),nx1_out,x1_out,var_out(1:nx1_out,j_x2))
ENDDO

RETURN
END SUBROUTINE INTERP2D

!========================================
! Linear x1-interpolation of f(x1,x2,x3)

SUBROUTINE INTERP3D(nx2,nx3,nx1_in,x1_in,var_in,nx1_out,x1_out,var_out)

use swim_global_data_mod, only: rspec

implicit none

integer          :: nx2, nx3, nx1_in, nx1_out, j_x2, j_x3
real             :: x1_in(nx1_in), var_in(nx1_in,nx2,nx3)
REAL(KIND=rspec) :: x1_out(nx1_out), var_out(nx1_out,nx2,nx3)

DO j_x2=1,nx2
  DO j_x3=1,nx3
    CALL INTERP1D(nx1_in,x1_in,var_in(1:nx1_in,j_x2,j_x3),nx1_out,x1_out,var_out(1:nx1_out,j_x2,j_x3))
  ENDDO
ENDDO

RETURN
END SUBROUTINE INTERP3D

!========================================
! Linear x1-interpolation of f(x1,x2,x3)

SUBROUTINE APPROX3D(nx2,nx3,nx1_in,x1_in,var_in,nx1_out,x1_out,var_out)

use swim_global_data_mod, only: rspec

implicit none

integer          :: nx2, nx3, nx1_in, nx1_out, j_x1, j1, jt(nx1_in)
real             :: x1_in(nx1_in), var_in(nx1_in,nx2,nx3),tim(nx1_in)
REAL(KIND=rspec) :: x1_out(nx1_out), var_out(nx1_out,nx2,nx3)

DO j_x1=1,nx1_out
  DO j1=1,nx1_in
    tim(j1)=ABS(x1_in(j1) - real(x1_out(j_x1)))
  ENDDO
  jt=MINLOC(tim)
  write(6,*) 'JT1'
  write(6,*) jt,tim(jt(1)),x1_in(jt(1)),x1_out(j_x1)
  var_out(j_x1,:,:)=var_in(jt(1),:,:)
ENDDO

RETURN
END SUBROUTINE APPROX3D
