module coordinate

  use tide_param, only: pi
  integer, save :: ny
  integer, save :: maxhgh
  integer, save :: nhgh = 10
  integer, save :: nxx=100
  real(8), save :: xx_min=0.0d0, xx_max=-1d0! xx_max = -log(p_top/p_bottom)
  real(8), allocatable, save :: xlat  (:)   ! latitudes
  real(8), allocatable, save :: theta (:)   ! latitudes
  real(8), allocatable, save :: sinth (:)   ! sin(theta)
  real(8), allocatable, save :: costh (:)   ! cos(theta)
  real(8), allocatable, save :: cotth (:)   ! cos(theta)
  real(8), allocatable, save :: wtheta(:)   ! gaussian weight in radian
  real(8), allocatable, save :: wmu   (:)   ! delta mu = delta ( costh )
  real(8), allocatable, save :: hough (:,:) ! latitudinal function shape, even mode
  real(8), allocatable, save :: eqdep (:)   ! equivalent depth [m]
  real(8), allocatable, save :: xx    (:)
  real(8), allocatable, save :: zz    (:)

  interface d_dtheta
    module procedure d_dtheta_real8, d_dtheta_cmplx8
  end interface
  interface d_dz
    module procedure d_dz_real8, d_dz_cmplx8
  end interface

contains

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

  subroutine coordinate__set_lat
    implicit none
    integer :: mtin, irc
    integer :: i
    character(len=128) :: file_hough='../../Hough/hough_s2_even_modes_tl95.txt'

    open(newunit=mtin,file=trim(file_hough),status='old',iostat=irc)
    if( irc/=0 )then
      write(*,*) 'Failed to open hough_s2.txt'
      stop 99
    endif

    read(mtin,*) maxhgh
    read(mtin,*) ny

    nhgh = min(maxhgh,nhgh)

    allocate( xlat(ny), wmu(ny), wtheta(ny) )
    allocate( hough(ny,maxhgh), eqdep(maxhgh) )

    read(mtin,*) xlat
    read(mtin,*) wmu

    do i=1, maxhgh
      read(mtin,*) eqdep(i)
      read(mtin,*) hough(:,i)
    enddo

    allocate( theta(ny), sinth(ny), costh(ny), cotth(ny) )
    do i=1, ny
      theta(i) = (90.0d0-xlat(i))*pi/180d0
      sinth(i) = sin(theta(i))
      costh(i) = cos(theta(i))
      cotth(i) = costh(i)/sinth(i)
      wtheta(i)= wmu(i)/sinth(i)
    enddo

    close(mtin)
  end subroutine coordinate__set_lat

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

  subroutine coordinate__set_xx( ps, pt )
    implicit none

    real(8), intent(in) :: ps, pt
    integer :: i

    ! determiine zz later
    allocate( xx(nxx), zz(nxx) )

    xx_max = -log(pt/ps)
    write(*,*) 'coordinate Xmax is : ', xx_max

    do i=1, nxx
      xx(i) = (xx_max - xx_min) * (i-1)/(nxx-1)
    enddo

  end subroutine coordinate__set_xx

  !--------------------------------------------------------------------
  function d_dtheta_real8( vvin ) result( vvout )
    implicit none
    real(8), intent(in) :: vvin(:,:)
    real(8)             :: vvout(size(vvin,1),size(vvin,2))
    integer :: iy, iz, nz
    real(8) :: dff, dtheta, dthetap, dthetam, dffp, dffm

    if( ny/=size(vvin,1) )then
      write(*,*) 'Y size error'
      stop 10
    endif

    nz = size(vvin,2)

    ! iy = 1
    dthetam = theta(2) - theta(1)
    dthetap = theta(3) - theta(1)
    do iz=1, nz
      dffm = vvin(2,iz) - vvin(1,iz)
      dffp = vvin(3,iz) - vvin(1,iz)
      vvout(1,iz) = (dthetap*dffm/dthetam - dthetam*dffp/dthetap ) &
                  &/(dthetap-dthetam)
    enddo

    ! iy = ny
    dthetam = theta(ny-1) - theta(ny)
    dthetap = theta(ny-2) - theta(ny)
    do iz=1, nz
      dffm = vvin(ny-1,iz) - vvin(ny,iz)
      dffp = vvin(ny-2,iz) - vvin(ny,iz)
      vvout(ny,iz) = (dthetap*dffm/dthetam - dthetam*dffp/dthetap ) &
                   &/(dthetap-dthetam)
    enddo

    do iz=1, nz
      do iy=2,ny-1

        dthetam = theta(iy-1) - theta(iy)
        dthetap = theta(iy+1) - theta(iy)
        dffm    = vvin(iy-1,iz) - vvin(iy,iz)
        dffp    = vvin(iy+1,iz) - vvin(iy,iz)

        vvout(iy,iz) = (dthetap*dffm/dthetam - dthetam*dffp/dthetap ) &
                     &/(dthetap-dthetam)
      enddo
    enddo

  end function d_dtheta_real8
  !--------------------------------------------------------------------
  function d_dtheta_cmplx8( vvin ) result( vvout )
    implicit none
    complex(8), intent(in) :: vvin(:,:)
    complex(8)             :: vvout(size(vvin,1),size(vvin,2))

    vvout = cmplx( d_dtheta(real(vvin)), d_dtheta(imag(vvin)) )

  end function d_dtheta_cmplx8
  !--------------------------------------------------------------------
  function d_dz_real8( vvin ) result( vvout )
    implicit none
    real(8), intent(in) :: vvin(:,:)
    real(8)             :: vvout(size(vvin,1),size(vvin,2))
    integer :: iy, iz, nz
    real(8) :: dff, dzz, dzzp, dzzm, dffp, dffm

    if( ny/=size(vvin,1) )then
      write(*,*) 'Y size error'
      stop 10
    endif

    nz = size(vvin,2)

    ! iz = 1
!    dzz = zz(2)-zz(1)
!    do iy=1, ny
!      dff = vvin(iy,2)-vvin(iy,1)
!      vvout(iy,1) = dff/dzz
!    enddo
    dzzm = zz(2)-zz(1)
    dzzp = zz(3)-zz(1)
    do iy=1, ny
      dffm = vvin(iy,2) - vvin(iy,1)
      dffp = vvin(iy,3) - vvin(iy,1)
      vvout(iy,1) = (dzzp*dffm/dzzm - dzzm*dffp/dzzp ) &
                   &/(dzzp-dzzm)
    enddo

    ! iz = nz
!    dzz = zz(nz)-zz(nz-1)
!    do iy=1, ny
!      dff = vvin(iy,nz)-vvin(iy,nz-1)
!      vvout(iy,nz) = dff/dzz
!    enddo
    dzzm = zz(nz-1)-zz(nz)
    dzzp = zz(nz-2)-zz(nz)
    do iy=1, ny
      dffm = vvin(iy,nz-1) - vvin(iy,nz)
      dffp = vvin(iy,nz-2) - vvin(iy,nz)
      vvout(iy,nz) = (dzzp*dffm/dzzm - dzzm*dffp/dzzp ) &
                   &/(dzzp-dzzm)
    enddo

    do iz=2, nz-1
      dzzm = zz(iz-1) - zz(iz)
      dzzp = zz(iz+1) - zz(iz)
      do iy=1,ny

        dffm = vvin(iy,iz-1) - vvin(iy,iz)
        dffp = vvin(iy,iz+1) - vvin(iy,iz)

        vvout(iy,iz) = (dzzp*dffm/dzzm - dzzm*dffp/dzzp ) &
                     &/(dzzp-dzzm)
      enddo
    enddo

  end function d_dz_real8
  !--------------------------------------------------------------------
  function d_dz_cmplx8( vvin ) result( vvout )
    implicit none
    complex(8), intent(in) :: vvin(:,:)
    complex(8)             :: vvout(size(vvin,1),size(vvin,2))

    vvout = cmplx( d_dz(real(vvin)), d_dz(imag(vvin)) )

  end function d_dz_cmplx8
  !--------------------------------------------------------------------
  subroutine zintp_spline_real8( m, ni, xxi, yyi, n, xx, yy  )
    implicit none
    ! 'i' for interpolation
    integer, intent(in)  :: m, ni
    real(8), intent(in)  :: xxi(ni)
    real(8), intent(out) :: yyi(m,ni)
    integer, intent(in)  :: n
    real(8), intent(in)  :: xx(n)
    real(8), intent(in)  :: yy(m,n)

    real(8) :: ztmp(n)
    integer :: k, j

    do j=1, m
      ztmp = yy(j,:)
      do k=1, ni
        if( xxi(k)>xx(n) )then
          yyi(j,k) = ztmp(n)
        else
          call spline_b_val( n, xx, ztmp, xxi(k), yyi(j,k) )
        endif
      enddo
    enddo

  end subroutine zintp_spline_real8
  !--------------------------------------------------------------------
  subroutine zintp_spline_cmplx8( m, ni, xxi, yyi, n, xx, yy  )
    implicit none
    ! 'i' for interpolation
    integer   , intent(in)  :: m, ni
    real(8)   , intent(in)  :: xxi(ni)
    complex(8), intent(out) :: yyi(m,ni)
    integer   , intent(in)  :: n
    real(8)   , intent(in)  :: xx(n)
    complex(8), intent(in)  :: yy(m,n)

    integer :: k, j
    real(8) :: zzr, zzi, ztmpr(n), ztmpi(n)

    do j=1, m
      ztmpr = real(yy(j,:))
      ztmpi = imag(yy(j,:))

      do k=1, ni
        call spline_b_val( n, xx, ztmpr, xxi(k), zzr )
        call spline_b_val( n, xx, ztmpi, xxi(k), zzi )
        yyi(j,k) = cmplx(zzr,zzi)
      enddo
    enddo

  end subroutine zintp_spline_cmplx8

end module coordinate