module tide_profile

  use tide_param, only: pi, gasr, grav, cp
  ! migrating semidiurnal tide
  real(8) :: h1eq=0.035d0, h2eq=0.36d0 ! ideal forcing
  real(8) :: j2tau=0.21d0*1.241329067669135   ! model tendency
  ! intg( J0*max(cos(theta),0) dtheta ) = J0/pi, therefore:
  real(8) :: mean_to_snapshot_J=pi

  integer             ,save :: nz
  real(8), allocatable,save :: zs(:), btt(:), buu(:), heat(:)
  real(8), allocatable,save :: btti(:), buui(:), heati(:), kn2(:)

contains

  !--------------------------------------------------------------------
  subroutine profile__init
    use tide_param, only: sgm
    use coordinate, only: nzi
    implicit none

    character(len=128) :: fin='../../Vert/lindzen_like_tjprof_normal.txt'

    integer :: mtin
    integer :: iz, irc
    real(8) :: dummy, v1, v2
    real(8), allocatable :: plev(:)

    !---------------------------------------------------------
    ! HEAT needs to be converted to Tau as in Jannssen (1999)
    !---------------------------------------------------------

    open(newunit=mtin,file=trim(fin),status='old',iostat=irc)
    if( irc/=0 ) stop 99

    call get_num_rec( mtin, nz )
    allocate( btt(nz), zs(nz), plev(nz), heat(nz), buu(nz) )

    if( index(fin,'lindzen')>0 )then ! ideal heating with O3 and H20

      do iz=1, nz
        !                z,     pres        tmp, H1, H2
!        read(mtin,*) gph(iz), plev(iz), btt(iz), v1, v2
        read(mtin,*)   dummy, plev(iz), btt(iz), v1, v2
        ! heati(=tau) = H2O absorption + O3 absorption
        heat(iz) = h1eq*v1 + h2eq*v2
      enddo

    else ! ESM2

      do iz=1, nz

!        read(mtin,*) dummy, gph(iz), plev(iz), btt(iz), heat(iz)
        read(mtin,*) dummy,   dummy, plev(iz), btt(iz), heat(iz)

      enddo

      ! convert from MEAN J to snapshot Tau
      ! see EQ144 in Lindzen(1970)
      heat(:) = heat(:)*mean_to_snapshot_J*j2tau/sgm

    endif

    if( plev(1)<800d2 .or. plev(1)>1050d2 )then
      write(*,*) 'P(1) has to be in Pa.'
      stop 99
    endif

    zs (:) = -log( plev(:)/plev(1) )
    buu(:) = 0.0d0

    close(mtin)

    ! variable allocation
    allocate( btti(nzi), buui(nzi), heati(nzi) )
    allocate( kn2(nzi) )

  end subroutine profile__init
  !--------------------------------------------------------------------
  subroutine calc_z_from_zst( ny, nz, zs, btt, zz )
    implicit none
    integer, intent(in)  :: ny, nz
    real(8), intent(in)  :: zs(ny), btt(ny,nz)
    real(8), intent(out) :: zz(ny,nz)

    integer :: iy, iz
    real(8) :: bttint, dzz

    do iy=1, ny
      zz(iy,1) = 0d0
    enddo

    do iz=2, nz-1
      do iy=1, ny
        bttint = 0.5d0*(btt(iy,iz)+btt(iy,iz-1))
        dzz    = zs(iz)-zs(iz-1)
        zz(iy,iz) = zz(iy,iz-1) + gasr*bttint/grav*dzz
      enddo
    enddo

  end subroutine
  !--------------------------------------------------------------------
  subroutine get_num_rec( mtin, nrec )
    implicit none
    integer, intent(in)  :: mtin
    integer, intent(out) :: nrec

    character(len=128)   :: header
    integer              :: irc

    ! read header
    read(mtin,'(A)',iostat=irc) header
    if(irc/=0) stop 88

    ! count num. of record
    nrec = 0
    do
      read(mtin,*,iostat=irc)
      if(irc/=0) exit
      nrec = nrec + 1
    enddo

    ! skip header
    rewind(mtin)
    read(mtin,*)
  end subroutine get_num_rec
  !--------------------------------------------------------------------
  subroutine profile__vint
    use coordinate, only: nzi, zsi
    integer :: i

    do i=1, nzi
      call spline_b_val( nz, zs, btt, zsi(i), btti(i) )
      call spline_b_val( nz, zs, heat , zsi(i), heati (i) )
    enddo

  end subroutine profile__vint
  !--------------------------------------------------------------------
  subroutine make_kn2(imode)
    use coordinate, only: nzi, zsi, eqdep
    implicit none
    integer, intent(in) :: imode
    integer :: i, i0, i1
    real(8) :: dx, ddx, hn

    dx = zsi(2)-zsi(1)

    hn = eqdep(imode)

    do i = 1, nzi
      i1 = min(i+1, nzi)
      i0 = max(i-1,   1)

      ddx = (i1-i0)*dx

      ! EQ7
      KN2(i) = gasr/grav/hn*( (btti(i1)-btti(i0))/ddx + gasr*btti(i)/cp ) - 0.25d0

    enddo

  end subroutine make_kn2
  !--------------------------------------------------------------------
end module tide_profile