module interpolate

  use vtype_module

  implicit none

  private

  public :: InterpolateHorVar


  interface InterpolateHorVar
    module procedure             &
      InterpolateHorVarDouble3D, &
      InterpolateHorVarDouble2D, &
      InterpolateHorVarReal3D, &
      InterpolateHorVarReal2D, &
      InterpolateHorVarInt3D, &
      InterpolateHorVarInt2D
  end interface InterpolateHorVar


contains

  !----------------------------------------------------------------------------

  real(dp) function InterpolateLinear2DDouble(   &
    & a_x, a_y, aa_z, x, y  &
    & ) result( Res )

    real(dp)    , intent(in ) :: a_x(2)
    real(dp)    , intent(in ) :: a_y(2)
    real(dp)    , intent(in ) :: aa_z(2,2)
    real(dp)    , intent(in ) :: x
    real(dp)    , intent(in ) :: y


    !
    ! Local variables
    !
    real(dp) :: a_z(2)
    integer(i4b) :: l


    do l = 1, 2
      a_z(l) = ( aa_z(2,l) - aa_z(1,l) ) &
        &    / ( a_x(2)    - a_x(1)    ) &
        &    * ( x         - a_x(1)    ) &
        &    + aa_z(1,l)
    end do
    Res = ( a_z(2) - a_z(1) ) / ( a_y(2) - a_y(1) ) * ( y - a_y(1) ) + a_z(1)


  end function InterpolateLinear2DDouble

  !----------------------------------------------------------------------------

  subroutine InterpolateHorVarDouble3D( &
    & imaxOrg, jmaxOrg, kmaxOrg,     &
    & x_LonOrg, y_LatOrg, & !z_SigOrg,  &
    & xyz_VarOrg,                    &
    & imax,    jmax,    kmax,        &
    & x_Lon,    y_Lat,    & !z_Sig,     &
    & xyz_Var                        &
    & )

    integer(i4b), intent(in ) :: imaxOrg
    integer(i4b), intent(in ) :: jmaxOrg
    integer(i4b), intent(in ) :: kmaxOrg
    real(dp)    , intent(in ) :: x_LonOrg(0:imaxOrg-1)
    real(dp)    , intent(in ) :: y_LatOrg(1:jmaxOrg  )
!!$    real(dp)    , intent(in ) :: z_SigOrg(kmaxOrg)
    real(dp)    , intent(in ) :: xyz_VarOrg(0:imaxOrg-1, 1:jmaxOrg, 1:kmaxOrg)
    integer(i4b), intent(in ) :: imax
    integer(i4b), intent(in ) :: jmax
    integer(i4b), intent(in ) :: kmax
    real(dp)    , intent(in ) :: x_Lon(0:imax-1)
    real(dp)    , intent(in ) :: y_Lat(1:jmax  )
!!$    real(dp)    , intent(in ) :: z_Sig(kmax)
    real(dp)    , intent(out) :: xyz_Var(0:imax-1, 1:jmax, 1:kmax)


    !
    ! Local variables
    !
    integer :: i
    integer :: j
    integer :: k
    integer :: ii
    integer :: jj
    integer :: kk

    real(dp) :: x_ExtLonOrg(0:imaxOrg  )
    real(dp) :: y_ExtLatOrg(0:jmaxOrg+1)
    real(dp) :: xyz_ExtVarOrg(0:imaxOrg, 0:jmaxOrg+1, 1:kmaxOrg)

    real(dp) :: Var1
    real(dp) :: Var2


    x_ExtLonOrg(0:imaxOrg-1) = x_LonOrg(0:imaxOrg-1)
    x_ExtLonOrg(imaxOrg    ) = 360.0
    y_ExtLatOrg(1:jmaxOrg  ) = y_LatOrg(1:jmaxOrg  )
    y_ExtLatOrg(0          ) = -90.0
    y_ExtLatOrg(jmaxOrg+1  ) =  90.0

    xyz_ExtVarOrg (0:imaxOrg-1, 1:jmaxOrg, 1:kmaxOrg) = &
      & xyz_VarOrg(0:imaxOrg-1, 1:jmaxOrg, 1:kmaxOrg)
    !
!!$    xyz_ExtVarOrg (0:imaxOrg  , 0        , 1:kmaxOrg) = &
!!$      & xyz_VarOrg(0:imaxOrg  , 1        , 1:kmaxOrg)
!!$    xyz_ExtVarOrg (0:imaxOrg  , jmaxOrg+1, 1:kmaxOrg) = &
!!$      & xyz_VarOrg(0:imaxOrg  , jmaxOrg  , 1:kmaxOrg)
    !
    xyz_ExtVarOrg(0:imaxOrg-1, 0        , 1:kmaxOrg) =        &
      &   (   xyz_VarOrg(0:imaxOrg-1, 2        , 1:kmaxOrg)   &
      &     - xyz_VarOrg(0:imaxOrg-1, 1        , 1:kmaxOrg) ) &
      & / ( y_LatOrg(2)    - y_LatOrg(1) )                    &
      & * ( y_ExtLatOrg(0) - y_LatOrg(1) )                    &
      & + xyz_VarOrg(0:imaxOrg-1, 1        , 1:kmaxOrg)
    do ii = 1, imaxOrg-1
      xyz_ExtVarOrg(0 , 0        , 1:kmaxOrg) =     &
        &   xyz_ExtVarOrg(0 , 0        , 1:kmaxOrg) &
        & + xyz_ExtVarOrg(ii, 0        , 1:kmaxOrg)
    end do
    xyz_ExtVarOrg(0, 0        , 1:kmaxOrg) =      &
      & xyz_ExtVarOrg(0, 0        , 1:kmaxOrg) / dble( imaxOrg )
    do ii = 0, imaxOrg
      xyz_ExtVarOrg (ii, 0        , 1:kmaxOrg) = &
        xyz_ExtVarOrg (0, 0        , 1:kmaxOrg)
    end do
    !
    xyz_ExtVarOrg (0:imaxOrg-1, jmaxOrg+1, 1:kmaxOrg) =       &
      &   (   xyz_VarOrg(0:imaxOrg-1, jmaxOrg  , 1:kmaxOrg)   &
      &     - xyz_VarOrg(0:imaxOrg-1, jmaxOrg-1, 1:kmaxOrg) ) &
      & / ( y_LatOrg   (jmaxOrg  ) - y_LatOrg(jmaxOrg-1) )  &
      & * ( y_ExtLatOrg(jmaxOrg+1) - y_LatOrg(jmaxOrg-1) )  &
      & + xyz_VarOrg(0:imaxOrg-1, jmaxOrg-1, 1:kmaxOrg)
    do ii = 1, imaxOrg-1
      xyz_ExtVarOrg(0 , jmaxOrg+1, 1:kmaxOrg) =       &
        &   xyz_ExtVarOrg(0 , jmaxOrg+1, 1:kmaxOrg)   &
        & + xyz_ExtVarOrg(ii, jmaxOrg+1, 1:kmaxOrg)
    end do
    xyz_ExtVarOrg(0 , jmaxOrg+1, 1:kmaxOrg) =       &
      & xyz_ExtVarOrg(0 , jmaxOrg+1, 1:kmaxOrg) / dble( imaxOrg )
    do ii = 0, imaxOrg
      xyz_ExtVarOrg(ii, jmaxOrg+1, 1:kmaxOrg) =       &
        & xyz_ExtVarOrg(0 , jmaxOrg+1, 1:kmaxOrg)
    end do
    !
    xyz_ExtVarOrg(imaxOrg    , 1:jmaxOrg, 1:kmaxOrg) = &
      & xyz_ExtVarOrg(0          , 1:jmaxOrg, 1:kmaxOrg)

    do k = 1, kmax
      do j = 1, jmax
        do i = 0, imax-1

          find_i : do ii = 1, imaxOrg-1
            if ( x_ExtLonOrg(ii) >= x_Lon(i) ) exit find_i
          end do find_i
          find_j : do jj = 1, jmaxOrg
            if ( y_ExtLatOrg(jj) >= y_Lat(j) ) exit find_j
          end do find_j
          kk = k


          Var1 = &
            &   ( xyz_ExtVarOrg(ii  ,jj-1,kk) - xyz_ExtVarOrg(ii-1,jj-1,kk) ) &
            & / ( x_ExtLonOrg  (ii  )         - x_ExtLonOrg  (ii-1)         ) &
            & * ( x_Lon        (i   )         - x_ExtLonOrg  (ii-1)         ) &
            & + xyz_ExtVarOrg(ii-1,jj-1,kk)
          Var2 = &
            &   ( xyz_ExtVarOrg(ii  ,jj  ,kk) - xyz_ExtVarOrg(ii-1,jj  ,kk) ) &
            & / ( x_ExtLonOrg  (ii  )         - x_ExtLonOrg  (ii-1)         ) &
            & * ( x_Lon        (i   )         - x_ExtLonOrg  (ii-1)         ) &
            & + xyz_ExtVarOrg(ii-1,jj  ,kk)

          xyz_Var(i,j,k) = &
            &   ( Var2 - Var1 )                           &
            & / ( y_ExtLatOrg(jj  ) - y_ExtLatOrg(jj-1) ) &
            & * ( y_Lat      (j   ) - y_ExtLatOrg(jj-1) ) &
            & + Var1
        end do
      end do
    end do


  end subroutine InterpolateHorVarDouble3D

  !----------------------------------------------------------------------------

  subroutine InterpolateHorVarReal3D( &
    & imaxOrg, jmaxOrg, kmaxOrg,     &
    & x_LonOrg, y_LatOrg, & !z_SigOrg,  &
    & xyz_VarOrg,                    &
    & imax,    jmax,    kmax,        &
    & x_Lon,    y_Lat,    & !z_Sig,     &
    & xyz_Var                        &
    & )

    integer(i4b), intent(in ) :: imaxOrg
    integer(i4b), intent(in ) :: jmaxOrg
    integer(i4b), intent(in ) :: kmaxOrg
    real(dp)    , intent(in ) :: x_LonOrg(imaxOrg)
    real(dp)    , intent(in ) :: y_LatOrg(jmaxOrg)
!!$    real(dp)    , intent(in ) :: z_SigOrg(kmaxOrg)
    real(sp)    , intent(in ) :: xyz_VarOrg(imaxOrg, jmaxOrg, kmaxOrg)
    integer(i4b), intent(in ) :: imax
    integer(i4b), intent(in ) :: jmax
    integer(i4b), intent(in ) :: kmax
    real(dp)    , intent(in ) :: x_Lon(imax)
    real(dp)    , intent(in ) :: y_Lat(jmax)
!!$    real(dp)    , intent(in ) :: z_Sig(kmax)
    real(sp)    , intent(out) :: xyz_Var(imax, jmax, kmax)


    !
    ! Local variables
    !
    real(dp) :: xyz_dVarOrg(imaxOrg, jmaxOrg, kmaxOrg)
    real(dp) :: xyz_dVar(imax, jmax, kmax)


    xyz_dVarOrg = xyz_VarOrg


    call InterpolateHorVar( &
      & imaxOrg, jmaxOrg, kmaxOrg,     &
      & x_LonOrg, y_LatOrg, & !z_SigOrg,  &
      & xyz_dVarOrg,                    &
      & imax,    jmax,    kmax,        &
      & x_Lon,    y_Lat,    & !z_Sig,     &
      & xyz_dVar                        &
      & )

    xyz_Var = xyz_dVar


  end subroutine InterpolateHorVarReal3D

  !----------------------------------------------------------------------------

  subroutine InterpolateHorVarInt3D( &
    & imaxOrg, jmaxOrg, kmaxOrg,     &
    & x_LonOrg, y_LatOrg, &!z_SigOrg,  &
    & xyz_VarOrg,                    &
    & imax,    jmax,    kmax,        &
    & x_Lon,    y_Lat,    &!z_Sig,     &
    & xyz_Var                        &
    & )

    integer(i4b), intent(in ) :: imaxOrg
    integer(i4b), intent(in ) :: jmaxOrg
    integer(i4b), intent(in ) :: kmaxOrg
    real(dp)    , intent(in ) :: x_LonOrg(imaxOrg)
    real(dp)    , intent(in ) :: y_LatOrg(jmaxOrg)
!!$    real(dp)    , intent(in ) :: z_SigOrg(kmaxOrg)
    integer(i4b), intent(in ) :: xyz_VarOrg(imaxOrg, jmaxOrg, kmaxOrg)
    integer(i4b), intent(in ) :: imax
    integer(i4b), intent(in ) :: jmax
    integer(i4b), intent(in ) :: kmax
    real(dp)    , intent(in ) :: x_Lon(imax)
    real(dp)    , intent(in ) :: y_Lat(jmax)
!!$    real(dp)    , intent(in ) :: z_Sig(kmax)
    integer(i4b), intent(out) :: xyz_Var(imax, jmax, kmax)


    !
    ! Local variables
    !
    integer :: i
    integer :: j
    integer :: k
    integer :: ii
    integer :: jj
    integer :: kk

    real(dp) :: xyz_dVarOrg(imaxOrg, jmaxOrg, kmaxOrg)
    real(dp) :: xyz_dVar(imax, jmax, kmax)


    xyz_dVarOrg = xyz_VarOrg


    call InterpolateHorVar( &
      & imaxOrg, jmaxOrg, kmaxOrg,     &
      & x_LonOrg, y_LatOrg, & !z_SigOrg,  &
      & xyz_dVarOrg,                    &
      & imax,    jmax,    kmax,        &
      & x_Lon,    y_Lat,    & !z_Sig,     &
      & xyz_dVar                        &
      & )

    xyz_Var = xyz_dVar


  end subroutine InterpolateHorVarInt3D

  !----------------------------------------------------------------------------

  subroutine InterpolateHorVarDouble2D( &
    & imaxOrg, jmaxOrg,              &
    & x_LonOrg, y_LatOrg,            &
    & xy_VarOrg,                     &
    & imax,    jmax,                 &
    & x_Lon,   y_Lat,                &
    & xy_Var                         &
    & )

    integer(i4b), intent(in ) :: imaxOrg
    integer(i4b), intent(in ) :: jmaxOrg
    real(dp)    , intent(in ) :: x_LonOrg(imaxOrg)
    real(dp)    , intent(in ) :: y_LatOrg(jmaxOrg)
    real(dp)    , intent(in ) :: xy_VarOrg(imaxOrg, jmaxOrg)
    integer(i4b), intent(in ) :: imax
    integer(i4b), intent(in ) :: jmax
    real(dp)    , intent(in ) :: x_Lon(imax)
    real(dp)    , intent(in ) :: y_Lat(jmax)
    real(dp)    , intent(out) :: xy_Var(imax, jmax)


    !
    ! Local variables
    !
    integer(i4b), parameter :: kmaxOrg = 1
!!$    real(dp) :: z_SigOrg(kmaxOrg)
    real(dp) :: xyz_VarOrg(imaxOrg, jmaxOrg, kmaxOrg)
    integer(i4b), parameter :: kmax    = 1
!!$    real(dp) :: z_Sig(kmax)
    real(dp) :: xyz_Var   (imax   , jmax   , kmax   )


!!$    z_SigOrg(1) = 0.0d0
!!$    z_Sig   (1) = 0.0d0

    xyz_VarOrg(:,:,1) = xy_VarOrg

    call InterpolateHorVar( &
      & imaxOrg, jmaxOrg, kmaxOrg,     &
      & x_LonOrg, y_LatOrg, & !z_SigOrg,  &
      & xyz_VarOrg,                    &
      & imax,    jmax,    kmax,        &
      & x_Lon,    y_Lat,    & !z_Sig,     &
      & xyz_Var                        &
      & )

    xy_Var = xyz_Var(:,:,1)


  end subroutine InterpolateHorVarDouble2D

  !----------------------------------------------------------------------------

  subroutine InterpolateHorVarReal2D( &
    & imaxOrg, jmaxOrg,              &
    & x_LonOrg, y_LatOrg,            &
    & xy_VarOrg,                     &
    & imax,    jmax,                 &
    & x_Lon,   y_Lat,                &
    & xy_Var                         &
    & )

    integer(i4b), intent(in ) :: imaxOrg
    integer(i4b), intent(in ) :: jmaxOrg
    real(dp)    , intent(in ) :: x_LonOrg(imaxOrg)
    real(dp)    , intent(in ) :: y_LatOrg(jmaxOrg)
    real(sp)    , intent(in ) :: xy_VarOrg(imaxOrg, jmaxOrg)
    integer(i4b), intent(in ) :: imax
    integer(i4b), intent(in ) :: jmax
    real(dp)    , intent(in ) :: x_Lon(imax)
    real(dp)    , intent(in ) :: y_Lat(jmax)
    real(sp)    , intent(out) :: xy_Var(imax, jmax)


    !
    ! Local variables
    !
    integer(i4b), parameter :: kmaxOrg = 1
!!$    real(dp) :: z_SigOrg(kmaxOrg)
    real(sp) :: xyz_VarOrg(imaxOrg, jmaxOrg, kmaxOrg)
    integer(i4b), parameter :: kmax    = 1
!!$    real(dp) :: z_Sig(kmax)
    real(sp) :: xyz_Var   (imax   , jmax   , kmax   )


!!$    z_SigOrg(1) = 0.0d0
!!$    z_Sig   (1) = 0.0d0

    xyz_VarOrg(:,:,1) = xy_VarOrg

    call InterpolateHorVar( &
      & imaxOrg, jmaxOrg, kmaxOrg,     &
      & x_LonOrg, y_LatOrg, & !z_SigOrg,  &
      & xyz_VarOrg,                    &
      & imax,    jmax,    kmax,        &
      & x_Lon,    y_Lat,    & !z_Sig,     &
      & xyz_Var                        &
      & )

    xy_Var = xyz_Var(:,:,1)


  end subroutine InterpolateHorVarReal2D

  !----------------------------------------------------------------------------

  subroutine InterpolateHorVarInt2D( &
    & imaxOrg, jmaxOrg,              &
    & x_LonOrg, y_LatOrg,            &
    & xy_VarOrg,                     &
    & imax,    jmax,                 &
    & x_Lon,   y_Lat,                &
    & xy_Var                         &
    & )

    integer(i4b), intent(in ) :: imaxOrg
    integer(i4b), intent(in ) :: jmaxOrg
    real(dp)    , intent(in ) :: x_LonOrg(imaxOrg)
    real(dp)    , intent(in ) :: y_LatOrg(jmaxOrg)
    integer(i4b), intent(in ) :: xy_VarOrg(imaxOrg, jmaxOrg)
    integer(i4b), intent(in ) :: imax
    integer(i4b), intent(in ) :: jmax
    real(dp)    , intent(in ) :: x_Lon(imax)
    real(dp)    , intent(in ) :: y_Lat(jmax)
    integer(i4b), intent(out) :: xy_Var(imax, jmax)


    !
    ! Local variables
    !
    integer(i4b), parameter :: kmaxOrg = 1
!!$    real(dp) :: z_SigOrg(kmaxOrg)
    integer(i4b) :: xyz_VarOrg(imaxOrg, jmaxOrg, kmaxOrg)
    integer(i4b), parameter :: kmax    = 1
!!$    real(dp) :: z_Sig(kmax)
    integer(i4b) :: xyz_Var   (imax   , jmax   , kmax   )


!!$    z_SigOrg(1) = 0.0d0
!!$    z_Sig   (1) = 0.0d0

    xyz_VarOrg(:,:,1) = xy_VarOrg

    call InterpolateHorVar( &
      & imaxOrg, jmaxOrg, kmaxOrg,     &
      & x_LonOrg, y_LatOrg, & !z_SigOrg,  &
      & xyz_VarOrg,                    &
      & imax,    jmax,    kmax,        &
      & x_Lon,    y_Lat,    & !z_Sig,     &
      & xyz_Var                        &
      & )

    xy_Var = xyz_Var(:,:,1)


  end subroutine InterpolateHorVarInt2D

  !----------------------------------------------------------------------------

end module interpolate
