      subroutine DISPER(bigKh,kh)

c**********************************************************************
c**** DESCRIPTION: Subroutine which calculates the root of         ****
c****              x = y tanh(y) using Newman's results from       ****
c****              Applied Ocean Research 1990, Vol 12, No 1       ****
c****              pp14-18.                                        ****
c**** USAGE:       On input:                                       ****
c****              bigKh:         Value of $Kh$                    ****
c****               (double)                                       ****
c****              On output:                                      ****
c****              kh:            Solution of dispersion relation  ****
c****               (double)                                       ****
c**** DEPENDENCIES:                                                ****
c**** AUTHORS:     R Porter                                        ****
c**** CHANGE_LOG:  Written 13/08/99                                ****
c**********************************************************************
c                                                     column 72 here -->
      implicit none

      double precision pi
      double precision bigKh,kh
      double precision x,y
      double precision fn,fny,fnyy

      common /const/   pi

      x  = bigKh

      if (x.lt.2.0d0) then
        y = dsqrt(x)*(0.9994d0+0.1701d0*x+0.0305d0*x**2)
      else if (x.gt.10.0d0) then
        y = x
        goto 10
      else
        y = x+2.0d0*x*dexp((-2.0d0)*x)-6.0d0*x**2*dexp((-4.0d0)*x)
      endif

c**********************************************************************
c**** 2 iterations of the higher order Newton required for         ****
c**** double precision                                             ****
c**********************************************************************

      fn   = 0.5d0*dlog((y+x)/(y-x))-y
      fny  = (-x)/(y*y-x*x)-1.0d0
      fnyy = 2.0d0*x*y/(y*y-x*x)**2
      y    = y-(fn/fny)*(1.0d0+0.5d0*(fn*fnyy/fny**2))

      fn   = 0.5d0*dlog((y+x)/(y-x))-y
      fny  = (-x)/(y*y-x*x)-1.0d0
      fnyy = 2.0d0*x*y/(y*y-x*x)**2
      y    = y-(fn/fny)*(1.0d0+0.5d0*(fn*fnyy/fny**2))

 10   kh = y

      return

      end
