      subroutine EIGVAL(bigKh,krh,rTrunc,maxR)

c**********************************************************************
c**** DESCRIPTION: Subroutine which calculates the first rTrunc    ****
c****              roots of -x = y tan(y) using Newman's results   ****
c****              from Applied Ocean Research 1990, Vol 12, No 1  ****
c****              pp14-18.                                        ****
c**** USAGE:       On input:                                       ****
c****              bigKh:         Value of $Kh$                    ****
c****               (double)                                       ****
c****              rTrunc:        Number of roots to be calculated ****
c****               (integer)                                      ****
c****              maxR:          Physical size of array krh       ****
c****               (integer)                                      ****
c****              On output:                                      ****
c****              krh:           Array of physical size maxR, the ****
c****               (integer)     first rTrunc entries contain the ****
c****                             roots of the dispersion relation ****
c**** DEPENDENCIES:                                                ****
c**** AUTHORS:     R Porter                                        ****
c**** CHANGE_LOG:  Written 13/08/99                                ****
c**********************************************************************
c                                                     column 72 here -->
      implicit none

      integer          r,rTrunc,maxR

      double precision pi,bigKh
      double precision krh(maxR)
      double precision x,y,u,delta0,delta1,deltan
      double precision fn,fny,fnyy

      pi = 3.141592653589793238462643d0
      x  = bigKh

c**********************************************************************
c**** now find $k_1h$, need.....                                   ****
c**********************************************************************

      u      = (3.0d0*x/(7.0d0+3.0d0*x))
      delta1 = 0.0159d0+0.1032d0*u+4.3152d0*u*u-2.8768d0*u*u*u
      y      = pi-delta1

c**********************************************************************
c**** calculate functions for the 1st iteration                    ****
c**********************************************************************

      fn     = datan(x/y)-delta1
      fny    = -1.0d0+x/(x*x+y*y)
      fnyy   = 2.0d0*x*y/((x*x+y*y)**2)
      deltan = delta1-(fn/fny)*(1.0d0+0.5d0*(fn*fnyy/fny**2))
      y      = pi-deltan

c**********************************************************************
c**** 2nd iteration                                                ****
c**********************************************************************

      fn     = datan(x/y)-deltan
      fny    = -1.0d0+x/(x*x+y*y)
      fnyy   = 2.0d0*x*y/((x*x+y*y)**2)
      deltan = deltan-(fn/fny)*(1.0d0+0.5d0*(fn*fnyy/fny**2))
      y      = pi-deltan
      krh(1) = y

c**********************************************************************
c**** now find the k_rh, r=2,3,...                                 ****
c**********************************************************************

      do 10 r=2,rTrunc
        delta0 = pi*dble(r-1)-krh(r-1)
        delta1 = delta0-pi*x/(x*x+pi*dble(r)*krh(r-1))
        y      = pi*dble(r)-delta1

c**********************************************************************
c**** 1st iteration                                                ****
c**********************************************************************

        fn     = datan(x/y)-delta1
        fny    = -1.0d0+x/(x*x+y*y)
        fnyy   = 2.0d0*x*y/((x*x+y*y)**2)
        deltan = delta1-(fn/fny)*(1.0d0+0.5d0*(fn*fnyy/fny**2))
        y      = pi*dble(r)-deltan

c**********************************************************************
c**** 2nd iteration                                                ****
c**********************************************************************

        fn     = datan(x/y)-deltan
        fny    = -1.0d0+x/(x*x+y*y)
        fnyy   = 2.0d0*x*y/((x*x+y*y)**2)
        deltan = deltan-(fn/fny)*(1.0d0+0.5d0*(fn*fnyy/fny**2))
        y      = pi*dble(r)-deltan
        krh(r) = y

 10   continue

      return

      end
