        parameter (nl=40,nt=100)
        real ss(nl+1),ff(nl+1),qq(nt+1),rr(nt+1)

        real hdpth(3)
        data hdpth / 12., 25., 50./
        open (11,file='kelvin',form='unformatted',
     *    access='direct',recl=(nl+1)*4,status='unknown')

        open (12,file='mrg',form='unformatted',
     *    access='direct',recl=(nt+1)*4,status='unknown')

        open (13,file='er',form='unformatted',
     *    access='direct',recl=(nl+1)*4,status='unknown')

        open (14,file='ig1',form='unformatted',
     *    access='direct',recl=(nl+1)*4,status='unknown')

        open (15,file='ig2',form='unformatted',
     *    access='direct',recl=(nl+1)*4,status='unknown')

        do i=1,nl+1
         ss(i)=float(i-1-nl/2)
        enddo 

        do irec=1,3
          call disperse_Kelvin (ss,nl+1,hdpth(irec),ff)
          write(11,rec=irec) ff
        enddo

        do j=1,nt+1
         qq(j)=(j-1)*0.005
        enddo 

c------ n = 0

        do irec=1,3
          call disperse_MRG(qq,nt+1,hdpth(irec),rr)
          write(12,rec=irec) rr
        enddo

         n = 1

         iwave =0

         do irec=1,3
           call disperse_ER_IG (ss,nl+1,hdpth(irec),ff,n,iwave)
           do i=nl/2+2,nl+1
              ff(i)=-999.
           enddo
           write(13,rec=irec) ff
         enddo

         iwave =1

         do irec=1,3
           call disperse_ER_IG (ss,nl+1,hdpth(irec),ff,n,iwave)
           write(14,rec=irec) ff
         enddo

         n = 2

         do irec=1,3
           call disperse_ER_IG (ss,nl+1,hdpth(irec),ff,n,iwave)
           write(15,rec=irec) ff
         enddo


        stop
        end
c-------------------------------------------------------------------------
        subroutine disperse_Kelvin (wn0,ii,H,f0)
        parameter (nx=40)

c--     input
        real H        ! equivalent depth
        real wn0(nx+1)  ! wavenumber (dimensionless)

c--     output
        real f0 (nx+1)  ! dispersion frequency [s-1]

c--     constants
        real g        ! gravity
        data g /9.8/

        real Re       ! Earth radius
        data Re /6.37e+6/

c--     local variables
        real k,omega,c
        real undef
        data undef /-999./
        real pi
        data pi /3.14159/

        c = (g*H)**0.5

        f0 = undef

        do i=ii/2+1,ii
          k=wn0(i)/Re
          omega=k*c
          f0(i)=omega/(2.0*pi)*3600.*24.
        enddo

        return
        end
c-------------------------------------------------------------------------
        subroutine disperse0 (H,n,wn,f,y,jj)
        parameter (nx=100)
c--     input
        real H        ! equivalent depth
        real wn       
        integer n 
        real f(nx)         

c--     output
        real y(nx)

c--     local variables
        real g,Re,pi,secday,beta
        real c,k,omega

        g = 9.8
        Re = 6.37e+6
        pi = 3.14159
        secday=86400.
        beta = 2*7.292E-5/Re 

        c = (g*H)**0.5

        do j=1,jj
          k=wn/Re
          omega=2.0*pi*f(j)/secday
          y(j)=(omega/c)**2.0-k**2.0-beta*k/omega-(2.0*n+1.0)*beta/c
        enddo

        return
        end
c-------------------------------------------------------------------------
        subroutine disperse_ER_IG (wn0,ii,H,f0,n,iwave)

c------- iwave = 0 ER
c------- iwave = 1 IG

        parameter (nx=40,nt=100)

c--     input
        real H        ! equivalent depth
        real wn0(nx+1)  ! wavenumber (dimensionless)
        integer n, iwave

c--     output
        real f0 (nx+1)  ! dispersion frequency [s-1]

c--     local variables
        real f00 (nt), y00 (nt)
        real k,omega,c
        integer kk(1)

        pi = 3.14159
        Re = 6.37e+6
        g = 9.8
        secday=86400.
        beta = 2*7.292E-5/Re 

        c = (g*H)**0.5

        do 10 j=1,nt
          f00(j)=float(j)*0.005
10      continue

        do 1000 i=1,ii

          wn=wn0(i)

          if (iwave .eq. 0 .and. wn .ge. 0 ) goto 90

          call disperse0 (H,n,wn,f00,y00,nt)

          kk = minloc (y00,dim=1)
          fa0 = f00 ( kk(1) )
          kk = maxloc (y00,dim=1)
          fb0 = f00 ( kk(1) )

          if (iwave .eq. 0 ) then
               fa = 0.0000001
               fb = 0.2
          endif
          if (iwave .eq. 1 ) then
               fa = 0.2
               fb = 10.0
          endif

          call disperse0 (H,n,wn,fa,ya,1)
          call disperse0 (H,n,wn,fb,yb,1)

80        continue
       
          f=(fa+fb)*0.5
          call disperse0 (H,n,wn,f,y,1)
          
          if (y*ya .gt. 0.0 ) then
              fa=f
              ya=y
          else
              fb=f
              yb=y
          endif
          df=abs(fa-fb)
          df_max=0.001

          if (df .gt. df_max) goto 80
          
          f0(i)=f

1000    continue

90      continue

        return
        end
c-------------------------------------------------------------------------
        subroutine disperse_MRG (f0,jj,H,wn0)
        parameter (ny=100)

c--     input
        real H        ! equivalent depth
        real f0 (ny+1)  ! dispersion frequency [s-1]

c--     output
        real wn0(ny+1)  ! wavenumber (dimensionless)

c--     local variables
        real temp (ny+1)  
        real wn1  (ny+1)  
        real wn2  (ny+1) 
        integer kk(1)

        real k,omega,c
        integer n

        pi = 3.14159
        Re = 6.37e+6
        g = 9.8
        secday=86400.
        beta = 2.0*7.292E-5/Re 

        c = (g*H)**0.5

        n = 0

        do j=1,jj
          omega = 2.0*pi*f0(j)/secday
          con1 = beta/omega

          wn1(j)=Re*(-con1-(con1**2.0-4*((2*n+1)*beta/c
     *                      -(omega/c)**2.0))**0.5)/2.0
          wn2(j)=Re*(-con1+(con1**2.0-4*((2*n+1)*beta/c
     *                      -(omega/c)**2.0))**0.5)/2.0

          temp(j)=abs(wn1(j)-wn2(j))
        enddo
        
        rmin = minval(temp,dim=1)
        kk = minloc(temp,dim=1)
        i=kk(1)

        do j=1,i
           wn0(j)=wn1(j)
        enddo
        do j=i+1,jj
           wn0(j)=wn2(j)
        enddo
        wn0(1)=-999.
        return
        end
