c234567
      program testana

C2  CALL      CALL TANAKA(fmax,n,m,alpha,Hd,X,Z,PHI,NF,ISAVE,IO6)       RAD00130
C3  CALL ARG. fmax                = Max. x/d to be considered (e.g. 40) RAD00140
C3            n                   = Number of nodes of FS (e.g. 80)     RAD00150
C3            alpha,m             = Stretching coef. near crest (5,0.01)RAD00160
C3            Hd                  = Wave amplitude H/d                  RAD00160
C3            TOLH                = Tolerance in the crest acc.(.00001)
C3            dXW                 = Step on the FS of eq. sp. BEM nodes
C3            NF                  = Actual number of BEM FS nodes
C3            N4=NF+4             = Working number of FS nodes
C3            ISAVE               = (0/1) existing stored grid (no/yes)
C3            ISF                 = File number for storing grid analys.
C3            IO6                 = File for results printing
C3            WK1,WK2,WK3         = Working variables
      PARAMETER (NOM=400,NN=80)
      IMPLICIT REAL*8(A-H,O-Z)                                          

      DIMENSION XG(NOM),ZG(NOM),PHI(NOM),PHIN(NOM)      

      COMMON /MAILLE/ NOMM,IDV1(3),NBPS(4),IBCOND(4),IBLIM(8),IDU2(23)
      COMMON /BOUNDC/ DUM3(8),GE,DE,DUM4(2),HSOL,EPSZ,IDW2(2) 
      COMMON /TANAK/ fmax,alpha,TOLH,Hd,dXW,h,m,ISAVE,ISF,IO6,NW


      DE = .310d0
      HSOL = .424*DE
c      HSOL = .82388*DE
      nbps(1)=80
      iblim(1)=1
      iblim(2)=nbps(1)
      xg(iblim(1))=0
      xg(iblim(2))=39
      ge=9.81
      epsz=.01
      call solwa(xg,zg,phi,phin,iodat)

      end

C----------------------------------------------------------------------
      	SUBROUTINE SOLWA(XG,ZG,PHI,PHIN,IODAT)
C----------------------------------------------------------------------
C
      PARAMETER (NOM=400,NN=80)
C
      IMPLICIT REAL*8(A-H,O-Z)                                          
c     IMPLICIT REAL*8(a-h,o-z)                                          
C
      CHARACTER *8 TEXTE
C
      COMMON /MAILLE/ NOMM,IDV1(3),NBPS(4),IBCOND(4),IBLIM(8),IDU2(23)
      COMMON /BOUNDC/ DUM3(8),GE,DE,DUM4(2),HSOL,EPSZ,IDW2(2) 
      COMMON /TANAK/ fmax,alpha,TOLH,Hd,dXW,h,m,ISAVE,ISF,IO6,NW
     .            
C
      DIMENSION XG(NOMM),ZG(NOMM),PHI(NOMM),PHIN(NOMM)
      DIMENSION XW(NOM),ZW(NOM),PHIW(NOM),PHINW(NOM),UW(NOM),WW(NOM), 
     .          WK3(NOM)
      DIMENSION x(NOM),y(NOM),dy(NOM),fai(NOM),q(NOM),th(NOM)
      DIMENSION ao(NN,NN),be(NN,NN),WK1(NN,NN),IWK(NN)
C
      DATA TEXTE/'SOLWAV01'/,ZERO/0.D0/
C
C.....Open data file for matrix saving. Determines if file already exists
C
      ISAVE = 0
      ISF   = 98
      OPEN(UNIT=ISF,STATUS='UNKNOWN',FORM='UNFORMATTED',FILE='mat.dat',
     .     ERR=5)
      GOTO 6
    5    ISAVE = 1
         OPEN(UNIT=ISF,STATUS='OLD',FORM='UNFORMATTED',FILE='mat.dat')
C
C.....Initialize dimensions and input solitary wave data
C
    6 NF     = NOM - 10
      N4     = NF  + 4
      NX     = 0
      Hd     = HSOL/DE
      d      = DE
      
      MX     = NBPS(1) - 1
      ZL0    = XG(IBLIM(2)) - XG(IBLIM(1))
      dXW    = ZL0/MX
      Xstart = XG(IBLIM(1))
      NFS    = MX + 1
      IO6    = IODAT
C
C     Fixed TANAKA's data
C
      n      = NN
      fmax   = 40.D0
      alpha  = 0.01D0
      m      = 5
      TOLH   = 1.D-05
      iflr   = 1
C
C.....Computation of s.w. profile and kinematics by Tanaka's method,
C     Interpolation of the results in a fixed x-grid of step dXW.
C     The final wave propagates right-left (x<0) and is max at x=0.
C
C      write(*,*) HSOL,EPSZ,GE,NOMM,NX
      CALL TANAKA(XW,ZW,PHIW,PHINW,ao,be,WK1,WW,WK3,IWK,   
C     -----------
     .            x,y,UW,dy,fai,q,th,NN,NF,N4,NX,n)
C
C.....Dimensionalize the results
C
      vchar = DSQRT(GE*d)
C
      DO i=1,NW
         XW(i)    = XW(i)*d
         ZW(i)    = ZW(i)*d
         UW(i)    = UW(i)*vchar
         WW(i)    = WW(i)*vchar
         PHINW(i) = PHINW(i)*vchar
         PHIW(i)  = PHIW(i)*d*vchar
      END DO
C
C.....Shift and truncate results where z< EPSZ*H
C
      Xmax  = Xstart + (NFS-1)*dXW
      WLmax = Xmax - Xstart
      Zmin  = EPSZ * Hd * d
C
C     Locate i of Zmin
C
      i=0
   10 i=i+1
         IF((i.LE.NW).AND.(ZW(i).LT.Zmin)) GOTO 10
      il=i-1
      ir=NW-il+1
C
C     Check horizontal extension
C
      WLENG = XW(ir) - XW(il)
c      IF((WLENG.GT.WLmax).OR.((ir-il+1).GT.NFS)) CALL ERRORS(TEXTE)
      IF((WLENG.GT.WLmax).OR.((ir-il+1).GT.NFS)) then
         write(*,*) TEXTE
         stop
      endif
C                                                -----------
C.....Calculate shift left in both cases (l-r, r-l)
C
      IF(Iflr.EQ.1) THEN
         jm = -1
	 ig = il
	 id = NW 
      ELSE
         jm = 1
	 ig = 1
	 id = ir
      END IF
C
      j=1
      Xshift = XW(ig) - Xstart
      DO i=ig,id
         XW(j)    = XW(i) - Xshift
         ZW(j)    = ZW(i)
         PHIW(j)  = PHIW(i) *jm
         PHINW(j) = PHINW(i)*jm
         UW(j)    = UW(i)*jm
         WW(j)    = WW(i)*jm
C
         j = j + 1
      END DO
C
      NW = id - ig + 1
C
C     Shift right if necessary (r-l)
C
      IF(Iflr.EQ.0) THEN
         j = NFS
         Xshift = XW(NW) - Xmax
         DO i=NW,1,-1
            XW(j)    = XW(i) - Xshift
            ZW(j)    = ZW(i)
            PHIW(j)  = PHIW(i)
            PHINW(j) = PHINW(i)
            UW(j)    = UW(i)
            WW(j)    = WW(i)
C
            j = j - 1
         END DO
      END IF
C
C.....Clear unused part of the wave in both cases
C
      IF(Iflr.EQ.1) THEN
         i1 = NW + 1
         i2 = NFS
         ic = NW
      ELSE
         i1 = 1
         i2 = NFS - NW
         ic = NFS - NW + 1
      END IF
C
      DO i=i1,i2
         XW(i)    = Xstart + (i-1)*dXW
         ZW(i)    = ZERO
         PHIW(i)  = PHIW(ic)
         PHINW(i) = ZERO
         UW(i)    = ZERO
         WW(i)    = ZERO
      END DO
C
C.....Write results
C
C
c      WRITE(IO6,2005)
C
      DO i=1,NFS
c         WRITE(IO6,2000) i,XW(i),ZW(i),PHIW(i),PHINW(i)
      END DO
C
c      WRITE(IO6,2010)
C
      DO i=1,NFS
c         WRITE(IO6,2000) i,XW(i),ZW(i),UW(i),WW(i)
      END DO
C
C.....Specify free surface data for BEM computations
C
      DO I=1,NFS
         XG(I)   = XW(I)
         ZG(I)   = ZW(I)
	 PHI(I)  = PHIW(I)
	 PHIN(I) = PHINW(I)
      END DO
C
      CLOSE(ISF)
C      
      RETURN
C
 2000 FORMAT(I5,4(F16.11,1X))
 2005 FORMAT(//'Print i,x,z,Phi,Phin'/)
 2010 FORMAT(//'Print i,x,z,u,v'/)
 2015 FORMAT(4F20.13)
C
      END
C---------------------------------------------------------------------- 
      SUBROUTINE TANAKA(XW,ZW,PHIW,PHINW,ao,be,WK1,WK2,WK3,IWK,     
     .                  x,y,dx,dy,fai,q,th,NN,NF,N4,NX,n)
C---------------------------------------------------------------------- 
C0  TANAKAF   TANAKA                                                    
C1  PURPOSE   Calculates steady solitary waves shape and potential
C1            by an iterative method based on Cauchy's integral theorem.
C1            (see Author Mitsuhiro Tanaka, Phys. Fluids 29(3), 1986).
C1               The program has been modified by S. Grilli for its
C1            introducing in the BEM nonlinear software. Full restru-
C1            ring and vector optimization of the program has been done.
CS  SYNOPSIS  :
CS  This program uses the above data to generate a solitary wave prof.
CS  with a given amplitude Hd. From that, an initial guess for qc, the
CS  crest velocity is found in table ATAB. qc lies between 0 and 1 when
CS  Hd varies between 0.833 and 0 (which gives a flat surface).
CS  Prior to calculating the waves, the "grid" for the points must be
CS  made. The point separation is governed by sep=a*x**m where a=0.01
CS  and m=5 in general, in the given above data.
CS  The grid data is saved in 'file solmat' for possible later reuse.
C1                                                                      RAD00120
C2  CALL      CALL TANAKA(fmax,n,m,alpha,Hd,X,Z,PHI,NF,ISAVE,IO6)       RAD00130
C3  CALL ARG. fmax                = Max. x/d to be considered (e.g. 40) RAD00140
C3            n                   = Number of nodes of FS (e.g. 80)     RAD00150
C3            alpha,m             = Stretching coef. near crest (5,0.01)RAD00160
C3            Hd                  = Wave amplitude H/d                  RAD00160
C3            TOLH                = Tolerance in the crest acc.(.00001)
C3            dXW                 = Step on the FS of eq. sp. BEM nodes
C3            NF                  = Actual number of BEM FS nodes
C3            N4=NF+4             = Working number of FS nodes
C3            ISAVE               = (0/1) existing stored grid (no/yes)
C3            ISF                 = File number for storing grid analys.
C3            IO6                 = File for results printing
C3            WK1,WK2,WK3         = Working variables
C3                                                                      RAD00170
C4  RET. ARG. XW,ZW,PHIW,PHINW    = Coord. and potent. and d/dn for s.w RAD00180
C6  INT. CALL SETMAT,DFSBC,SOLUNI,ERRORS                                RAD00190
C9  Sept 89   S. GRILLI, Civil Engng., Univ. of Delaware (Copyright)    RAD00200
C9
CLTANAKA  SUB. WHICH COMPUTES EXACT SOLITARY WAVES SHAPE AND POTENTIAL  RAD00210
C-----------------------------------------------------------------------RAD00220
C                                                                       RAD00230
      PARAMETER (ITRFR2=70, ITRH=10)                                    RAD00240
C                                                                       RAD00250
      IMPLICIT REAL*8(A-H,O-Z)                                          RAD00260
c     IMPLICIT REAL*8(a-h,o-z)                                          RAD00260
C                                                                       RAD00270
      CHARACTER *8 TEXTE
C      
      COMMON/TANAK/ fmax,alpha,TOLH,Hd,dXW,h,m,ISAVE,ISF,IO6,NW

C                                                                       RAD00290
      DIMENSION XW(NF),ZW(NF),PHIW(NF),PHINW(NF),SCOE(5),ATAB(0:50)     RAD00360
      DIMENSION ao(NN,NN),be(NN,NN),WK1(NN,NN),WK2(NF),WK3(NF),S(10)
      DIMENSION x(NF),y(NF),dx(-3:N4),dy(-3:N4),fai(NF),q(NF),th(NF)
      DIMENSION IWK(NN)
C                                                                       RAD00380
      DATA ZERO/0.D0/,ONE/1.D0/,TWO/2.D0/,TEXTE/'TANAKA01'/,HALF/.5D0/, RAD00390
     .     THREE/3.D0/
C                                                                       RAD00400
C.....look up table for amplitudes: qc=0. to 1. by 0.02
C
      DATA ATAB /0.833197D0,0.832926D0,0.831829D0,0.830159D0,0.827981D0,
     .0.825356D0,0.822278D0,0.818686D0,0.814479D0,0.809538D0,0.803735D0,
     .0.796952D0,0.789083D0,0.780045D0,0.769777D0,0.758244D0,0.745435D0,
     .0.731359D0,0.716050D0,0.699555D0,0.681936D0,0.663267D0,0.643628D0,
     .0.623107D0,0.601793D0,0.579778D0,0.557152D0,0.534004D0,0.510421D0,
     .0.486484D0,0.462273D0,0.437861D0,0.413316D0,0.388704D0,0.364083D0,
     .0.339507D0,0.315025D0,0.290682D0,0.266518D0,0.242568D0,0.218864D0,
     .0.195434D0,0.172301D0,0.149486D0,0.127004D0,0.104869D0,0.083090D0,
     .0.061682D0,0.040666D0,0.020083D0,0.000000D0/
C
C.....General constant data
C
      DATA SCON/7257600.D0/,SCOE/2497.D0,-28939.D0,162680.D0,
     .     -641776.D0,4134338.D0/,P02/.02D0/,P10/1.D-10/
C
C.....Input existing mesh data or generate mesh and Cauchy's matrix
C     and its explicit inverse (Crout's method)
C
      IF(ISAVE.EQ.1) THEN
         READ(ISF) fmax,alpha,h,n,m
         READ(ISF) ((ao(i,j),i=1,n),j=1,n)
      ELSE
C
C        Create grid data
C
         CALL SETMAT(ao,be,WK1,WK2,IWK,NN,NX,n)
C        -----------
      END IF
C
C.....Set up 10 pt lagrangian interpolation formula
C
      DO l=1,5
         S(l)=SCOE(l)/SCON * h
         S(11-l)=S(l)
      END DO
C
C.... Set up initial guess for qc as fct. of Hd from table ATAB
C     Keep the previous table element (qc0,h0) for later re-correction
C
      DO i=1,50
C
         ampp=ATAB(i)
C
         IF (Hd.gt.ampp) THEN
C
            qc=P02*(i-(Hd-ampp)/(ATAB(i-1)-ampp))
C
            IF ((i*P02-qc).GT.(qc-(i-1)*P02)) THEN
              qc0=i*P02
              h0=ATAB(i)
            ELSE
              qc0=(i-1)*P02
              h0=ATAB(i-1)
            END IF
            GOTO 10
         END IF
      END DO
 
   10 CONTINUE
C
C.....Iterations on the amplitude Hd :
C     Repeat until err. smaller than TOL or iterations over ITRH
C
      iterH=1
C
   15 CONTINUE
C
C     ...Initial guess for tau=log(q)
C        Potential values according to stretching formula
C
         DO i=1,n
            g=(i-1)*h
            fai(i)=g**m+alpha*g
         END DO
C
C        Crest value of q at first point
C
         q(1)=DLOG(qc)
C
C        Initial guess fct. for tau, placed in vector q : ln(1-(1-qc)/e)
C
         DO i=2,n
            q(i)=DLOG(ONE-(ONE-qc)/DEXP(fai(i)))
         END DO
C
C        Initial Froude number **2
C
         Fr2=ZERO
C
C        Iterations on the nonlinear FSBC :
C        Repeat until err. smaller than 1D-10 or iterations over ITRFR2
C
         iters=1
C
   20    CONTINUE
C
C        ...Compute fr2 corrected and new (q,th). Evaluate error on Fr2
C
            CALL DFSBC(error,Fr2,ao,dy,q,th,S,qc,NN,N4,n)
C           ----------
C
            iters=iters+1
C
            IF((error.GT.P10).AND.(iters.LE.ITRFR2)) GOTO 20
C
C     ...Final amplitude and crest velocity based on final Fr2
C
c23456
         ampc=(ONE-qc*qc)*Fr2*HALF
C
         IF ((DABS(Hd-ampc).LE.TOLH).OR.(iterH.GT.ITRH)) GOTO 25
C
C     ...Corrected qc value for restarting calculations for Hd by
C        re-interpolating in ATAB, for the corrected amplitude
C
         qc=qc0+(Hd-h0)/(ampc-h0)*(qc-qc0)
         iterH=iterH+1

C
      GOTO 15
C
C.....End of iterations
C
   25 CONTINUE
C
C.....Final Froude number and energy
C
      Fr=sqrt(Fr2)
C
C.....Calculation of x(j),y(j) (x(1)=0. and y=0. at infin.), from tau,q
C     NB: Fixed x-steps h are used for the Cauchy's integral colloca-
C     tion, coupled with a stretching for Phi. The actual x,y are then
C     built from the stret. formula and lagrangian interpolation.
C
      DO j=1,n
         g=(j-1)*h
         dfdg=m*g**(m-1) + alpha
         q(j) = DEXP(q(j))
         dx(j)= DCOS(th(j))/q(j)*dfdg
         dy(j)=-DSIN(th(j))/q(j)*dfdg
      END DO
C
C     Steps dx,dy are deduced from dW/dz=q(cos(th)+i sin(th)) and
C     lagrangian smoothing (with 4 extra extremity values of (dx,dy)).
C
C     Initial values of x,y at wave crest. Extreme values of dx,dy
C
      y(1)=(ONE-qc*qc)*Fr2*HALF
      x(1)=ZERO
C
      dy(1)=ZERO
      dy(0)=-dy(2)
      dy(-1)=-dy(3)
      dy(-2)=-dy(4)
      dy(-3)=-dy(5)
      dy(n+1)=ZERO
      dy(n+2)=ZERO
      dy(n+3)=ZERO
      dy(n+4)=ZERO
C
      dx(0)=dx(2)
      dx(-1)=dx(3)
      dx(-2)=dx(4)
      dx(-3)=dx(5)
C
      DO i=n+1,n+4
         g=(i-1)*h
         dx(i)=m*g**(m-1) + alpha
      END DO
C
C     Smoothing of actual (dx,dy). Computation of x,y by summing (dx,dy)
C
      DO j=1,n-1
         delx=ZERO
         dely=ZERO
C
         DO l=1,10
            delx=delx+S(l)*dx(j+l-5)
            dely=dely+S(l)*dy(j+l-5)
         END DO
C
         x(j+1)=x(j)+delx
         y(j+1)=y(j)+dely
      END DO
C
C.....Check of the accuracy of Bernoulli's equation on FS
C
      p0=q(1)*q(1)*HALF+y(1)/Fr2
      err=ZERO
C
      DO j=2,n
         err2=DABS(q(j)*q(j)*HALF+y(j)/Fr2-p0)
C
         IF(err2.GT.err) THEN
            err=err2
         END IF
      END DO
C
C.....Corrected y for exactly satisfying Bernoulli's equation
C
      DO i=1,n
         y(i)=(ONE-q(i)*q(i))*Fr2*HALF
      END DO
C
C.....Calculation of excess mass (em), impulse (rimp) and energy(tv,t,v)
C     The integrals are calculated by trapezoidal formula.
C
      t =ZERO
      em=ZERO
C
      DO j=1,n
         g=(j-1)*h
         dfdg=m*g**(m-1) + alpha
         del=y(j)*dfdg*h
         t=t+del
         del=del*DCOS(th(j))/q(j)
         em=em+del*TWO
C
C        Extremity correction in trapez. formula
C
         IF(j.EQ.1) THEN
            t=t*HALF
            em=em*HALF
         END IF
      END DO
C
      rimp=Fr*em
      v = (Fr2-ONE)*em/THREE
      t =-(t-em*HALF)*Fr2
      tv= t + v
C
C.....Calculation for use in SOLUNI of velocity components
C     (u,v) (in dx,dy) and potential Phi, in dimensional form.
C     The wave propagates right-left: u>0,v<0,f>0, x,y>0 in the
C     moving frame.
C
      DO i=1,n
         dx(i) = q(i)*DCOS(th(i))*Fr
         dy(i) =-q(i)*DSIN(th(i))*Fr
         fai(i)= fai(i)*Fr
      END DO
C
C.....Saving of the check results and initial data on file IO6
C
      write(IO6,2300)
      write(IO6,2350) fmax,n
      write(IO6,2400) alpha,m
      write(IO6,2450)
      write(IO6,2500) ampc,qc
      write(IO6,2550) Fr
      write(IO6,2600) em,rimp
      write(IO6,2650) t,v,tv
      write(IO6,2700) iters,iterH,error
C
C.....Processing by SOLUNI for getting the wave at BEM discretization :
C     i.e. coordinates (XW,ZW) in a constant x-step dXW grid and
C     potential PHIW and its normal gradient PHINW.
C
      NW=n
      N5=NF+5
C
      OPEN(47,FILE='maple')
      write(*,*) 'coucou'
      write(47,*) 'xyfuv:=['
      do i =1,n-1
         write(47,*) '[',x(i),',',y(i),',',
     +        fai(i),',',dx(i),',',dy(i),'],'
      enddo
      write(47,*) '[',x(n),',',y(n),',',
     +     fai(n),',',dx(n),',',dy(n),']];'
      close(47)

      do i=1,n
         write(*,*) i,x(i),y(i),-fai(i)+fr*x(i),fr-dx(i),-dy(i)
      enddo
      CALL SOLUNI(x,y,fai,dx,dy,WK2,WK3,NW,Fr,dXW,XW,ZW,PHIW,PHINW,N5)
C     -----------
C
      RETURN
C
 2300 FORMAT('input data:')
 2350 FORMAT(1x,'fmax=',F7.4,5x,'n=',I4)
 2400 FORMAT(1x,'alpha=',F11.4,', m=',I2)
 2450 FORMAT(//,1x,'output data:')
 2500 FORMAT(/,1x,'amp=',F16.10,5x,'qc=',F16.10)
 2550 FORMAT(1x,'celerity=',F16.10)
 2600 FORMAT(/,1x,'mass=',F16.10,5x,'impulse=',F16.10)
 2650 FORMAT(1x,'k.e.=',F16.10,5x,'p.e.=',F16.10,5x,'energy=',F16.10)
 2700 FORMAT(/,1x,'iters=',I3,5x,'iterH=',I3,5x,'error=',F15.10)
C
      END
C----------------------------------------------------------------------
      SUBROUTINE DFSBC(error,Fr2,ao,si,q,th,S,qc,NN,N4,n)
C----------------------------------------------------------------------
C
      IMPLICIT REAL*8 (A-H,O-Z)
c     IMPLICIT REAL*8 (a-h,o-z)                                          
C
      COMMON/TANAK/ fmax,alpha,TOLH,Hd,dXW,h,m,ISAVE,ISF,IO6,NW

      DIMENSION ao(NN,n),si(-3:N4),q(n),th(n),S(10)
C
      DATA ZERO/0.D0/,THREE/3.D0/,ONE/1.D0/
C
C.....Save old Froude number
C
      oFr2=Fr2
C
C.....Compute teta and Bernoulli integ. si on the free surface (q=tau !)
C
      DO i=1,n
         thi=ZERO
C
C        Application of Cauchy's integ. th. = transf. fct. teta=f(tau)
C
         DO j=1,n
            thi=thi+ao(i,j)*q(j)
         END DO
C
         g=(i-1)*h
         si(i)=DSIN(thi)*THREE*(m*g**(m-1)+alpha)
         th(i)=thi
      END DO
C
      si(1)  = ZERO
      si(0)  =-si(2)
      si(-1) =-si(3)
      si(-2) =-si(4)
      si(-3) =-si(5)
      si(n+1)= ZERO
      si(n+2)= ZERO
      si(n+3)= ZERO
      si(n+4)= ZERO
C
C.....Final Bernoulli integ. computed in q by lagrangian integration
C
      q(1)=ZERO
C
      DO j=1,n-1
C
         del=q(j)
         DO l=1,10
            del=del+S(l)*si(j+l-5)
         END DO
         q(j+1)=del
      END DO
C
C.....Computes new Froude number based on extremity Bernoulli integral
C
      qc3=qc**3
      Fr2=q(n)/(ONE-qc3)
C
C.....Corrected tau (in q) on FS based on Bernoulli. Error on Fr2
C
      DO j=1,n
         q(j)=DLOG(qc3+q(j)/Fr2)/THREE
      END DO
C
      error=DABS(Fr2-oFr2)
C
      RETURN
C
      END
C---------------------------------------------------------------------
      SUBROUTINE SETMAT(ao,be,vwm,vw,IWK,NN,NX,n)     
C---------------------------------------------------------------------
C
      IMPLICIT REAL*8 (A-H,O-Z)
c     IMPLICIT REAL*8 (a-h,o-z)                                          
C
      COMMON/TANAK/ fmax,alpha,TOLH,Hd,dXW,h,m,ISAVE,ISF,IO6,NW

      DIMENSION ao(NN,NN),be(NN,NN),vw(NN),vwm(NN,NN),SCOE(5),S(10),
     .          IWK(NN)
C
      DATA rpi/3.1830988618379067153D-1/,ZERO/0.D0/,ONE/1.D0/,
     .     HALF/.5D0/,TWO/2.D0/,FOUR/4.D0/
      DATA SCOE/35.D0,-405.D0,2268.D0,-8820.D0,39690.D0/,SCON/65536.D0/
C
C.....Initialization
C
      DO l=1,5
         S(l)    = SCOE(l)/SCON
         S(11-l) = S(l)
      END DO
C
C     Free surface constant comput. step, and Phi max
C
      gmax=fmax**(ONE/m)
      h=gmax/n
C
C.....Fill up ao (tempor.) with lagrangian interpol. polynomial
C
      DO i=1,n
         DO j=1,n
            ao(i,j) = ZERO
         END DO
      END DO
C
      DO i=1,n
         DO l=1,10
            IF ((i+l-5.LE.n).AND.(i+l-5.GT.0)) THEN
               ao(i,i+l-5)=S(l)
            END IF
         END DO
      END DO
C
      DO i=1,4
         DO l=2,6-i
            ao(i,l)=ao(i,l)+S(7-i-l)
         END DO
      END DO
C
C.....Cauchy's matrix be set up by integration with trapezoidal rule
C
      DO i=1,n
         DO j=1,n
            gr=(j-1)*h + h*HALF
            g0=(i-1)*h
            f0=g0**m + alpha*g0
            fr=gr**m + alpha*gr
            gr=h*(m*gr**(m-1)+alpha)/(fr-f0)*rpi
            gl=gr*(fr-f0)/(-fr-f0)
C
            be(i,j)=gr+gl
C
         END DO
      END DO
C
C.....Matrix multiplication : be=be*ao
C
      DO i=1,n
         DO j=1,n
            vw(j)=ZERO
C
            DO k=1,n
               vw(j)=vw(j)+be(i,k)*ao(k,j)
            END DO
         END DO
C
         DO j=1,n
            be(i,j)=vw(j)
         END DO
      END DO
C
C.....Integration by trapezoidal rule for ao,be
C
      DO i=1,n
         DO j=1,n
            gj=(j-1)*h
            gi=(i-1)*h
            fj=gj**m + alpha*gj
            fi=gi**m + alpha*gi
            dfdg=m*gj**(m-1) + alpha
            fjmfi=fj-fi
            fjpfi=fj+fi
            bm=dfdg/(fjmfi*fjmfi+FOUR)*TWO*rpi*h
            bp=dfdg/(fjpfi*fjpfi+FOUR)*TWO*rpi*h
            cm=-bm*(fjmfi)*HALF
            cp= bp*(fjpfi)*HALF
            dao=bm-bp
            dbe=cm+cp
C
            IF(j.EQ.1) THEN
               dao=dao*HALF
               dbe=dbe*HALF
            END IF
C
            ao(i,j)=dao
            be(i,j)=be(i,j)+dbe
C
         END DO
C
C        Accounts for Cauchy's free term
C
         ao(i,i)=ao(i,i)+ONE
C
      END DO
C
C.....Matrix ao inversion by NAG routine : ao**-1 (result is in vwm())
C
      DO i = 1, n
         DO j = 1, n
             vwm(i,j) = ZERO
         END DO  
         vwm(i,i) = ONE
      END DO
C
      nc = n
      NMA = NN
      CALL LUDCTA(ao,nc,NMA,IWK,dval)
C     -----------

      DO j = 1, n
         CALL LUBKTA(ao,nc,NMA,IWK,vwm(1,j))
C        -----------
      END DO
C
C.....Matrix multiplication : ao=vwm*be
C
      DO i=1,n
         DO j=1,n
C
            ao(i,j)=ZERO
C
            DO k=1,n
               ao(i,j)=ao(i,j)+vwm(i,k)*be(k,j)
            END DO
         END DO
      END DO
C
C.....Save results for later use
C
      WRITE(ISF) fmax,alpha,h,n,m
      WRITE(ISF) ((ao(i,j),i=1,n),j=1,n)
C
      ISAVE=1
C
      RETURN
C
      END
C----------------------------------------------------------------------
      SUBROUTINE SOLUNI(z,w,p,  u, q, v, s,  n, Fr,dxu,x, y, f,   fn,N5)
c      CALL      SOLUNI(x,y,fai,dx,dy,WK2,WK3,NW,Fr,dXW,XW,ZW,PHIW,PHINW,N5)
C----------------------------------------------------------------------
C      converts 'steady data d' files from solits to 'chyin data d'
C
C----------------------------------------------------------------------
C
      IMPLICIT REAL *8 (A-H,O-Z)
c      IMPLICIT REAL *8 (a-h,o-z)                                          
C
      CHARACTER *8 TEXTE
C
      DIMENSION x(-4:N5), y(-4:N5), f(-4:N5), fn(-4:N5), u(-4:N5),
     .          v(-4:N5), z(-4:N5), w(-4:N5), p(-4:N5) , q(-4:N5),
     .          s(-4:N5), a(0:10)
C
      DATA ZERO/0.D0/,TEN/10.D0/,HALF/0.5D0/,TEXTE/'SOLUNI01'/,
     .     TWO/2.D0/,TOLxj/1.D-12/
C
C.....Initialisation
C
C     Moves the working field 5 and 1 positions forward (for lagr.)
C
      DO i=n,1,-1
         z(i)=z(i-5)
         w(i)=w(i-5)
         p(i)=p(i-5)
         u(i)=u(i-1)
         q(i)=q(i-1)
      END DO

C     Horizontal velocity at inf.+-=Fr>0 for later norm. (prp. l-r)
C
      uinf=u(n)
C
C     Normalise first point (peak) to be at  x=0  where  f=0:
C
      z1=z(1)
      f1=p(1)
C
      DO i=1,n
         z(i)=z(i)-z1
         p(i)=p(i)-f1
      END DO
C
C     Initialize for lagrangian interpolation (5 initial <0 values)
C
      DO i=1,5
         z(1-i)=-z(1+i)
         w(1-i)= w(1+i)
         p(1-i)=-p(1+i)
         u(1-i)= u(1+i)
         q(1-i)=-q(1+i)
      END DO
C
C     Renormalise initial values (i=1) to the exact ones
C
      x(1) = ZERO
      y(1) = w(1)
      f(1) = ZERO
      s(1) = u(1)
      v(1) = ZERO
C
C.....General loop on x(j)=xx (new fixed step x) by interpolation on FS
C
      j=1
      nx=1
C
   20 CONTINUE
C
         j=j+1
         IF (2*(j+1).GT.N5) THEN
C	    WRITE(*,*) j
C	    CALL ERRORS(TEXTE)
C           -----------
            GOTO 50
         END IF
C	 
C     ...New x(j)=xx of BEM discretization
C
         xx=x(j-1)+dxu
         x(j)=xx
C
C     ...Loop on original values x(i)=z(i) for locating x(j)
C
         i=nx-1
C
   30    i=i+1
         IF (i.GT.(n-5)) GOTO 50
            nx=i
            xm=HALF*(z(i+1)+z(i))
            IF (xm.LT.xx) GOTO 30
C
C        ...Fits 11-point polynomial to z(i) : pol(xj)=S(a(i)*xj**i)
C           with i=0 to 10 and xj in [-1,1].
C
            CALL POL11I(z,nx,a)
C           -----------
C
C        ...Loop on xj :
C           Estimate of the initial xj for running Newton's interat.
C           for finding xj such that pol11(xj)=xx or F(xj)=pol(xj)-xx
C           Stop after 8 iterations
C
            ixj=0
            xj=TWO*(xx-z(nx))/(z(nx+1)-z(nx-1))
C
   40       CONTINUE
C
               xjp=xj
               Fxj  = a(10)*xj
               dFxj = Fxj*TEN
C
               DO i=9,2,-1
                  Fxj  = (Fxj + a(i))*xj
                  dFxj = (dFxj + a(i)*DFLOAT(i))*xj
               END DO
C
               Fxj  = (Fxj + a(1))*xj + a(0) - xx
               dFxj = dFxj + a(1)
C
               xj = xjp - Fxj/dFxj
               ixj= ixj + 1
C
            if ((ixj.LT.8).AND.(DABS(xj-xjp).GT.TOLxj)) GOTO 40
C
C        ...Computes the other fields at xj of the new x(j) on FS
C
            y(j)=POLY11(xj,w,nx)
C                ------
            f(j)=POLY11(xj,p,nx)
C                ------
            v(j)=POLY11(xj,q,nx)
C                ------
            s(j)=POLY11(xj,u,nx)
C                ------
         GOTO 20
C
   50 CONTINUE
C
C.....Symmetrization of the grid, imposition of zero velocity at end
C
C     Number of points on half the free surface
C
      n=j-1
C
      n1=n+1
      nl=n-1
C
C     Generates total solitary wave propagating towards x<0 (use
C     z,w,p,u,q as working variables  x,y,f,u>0 v<0)
C
      DO i=1,n
         z(n1-i)=-x(i)
         z(nl+i)= x(i)
         w(n1-i)= y(i)
         w(nl+i)= y(i)
         p(n1-i)=-f(i)
         p(nl+i)= f(i)
         u(n1-i)= s(i)
         u(nl+i)= s(i)
         q(n1-i)=-v(i)
         q(nl+i)= v(i)
      END DO
C
C.....Final FS grid dimension n ignores last points of 0 velocity
C
      n=2*n-1
C
C.....Computation of normal velocity fn(i) in the steady frame at fr
C     Correction of the velocities to account for the steady frame
C
      DO i=1,n
         vm=DSQRT(u(i)*u(i)+q(i)*q(i))
         sint=q(i)/vm
         fn(i)= Fr*sint
      END DO
C
C.....Store working variables into actual values. Correct for moving
C     frame by subtracting uinf of u, and uinf*x of f. The final wave
C     propagates from right (x>0) to left (x<0).
C
      DO i=1,n
         x(i)=z(i)
         y(i)=w(i)
         v(i)=q(i)
         u(i)=u(i)-uinf
         f(i)=p(i)-uinf*x(i)
      END DO
C
C.....Move results 5 backward for transmission
C
      DO i=1,n
         x (i-5)=x (i)
         y (i-5)=y (i)
         f( i-5)=f (i)
         fn(i-5)=fn(i)
         u (i-5)=u (i)
         v (i-5)=v (i)
      END DO
C
      RETURN
C
      END
C---------------------------------------------------------------------
      FUNCTION POLY11(x,f,n)
C---------------------------------------------------------------------
C
      IMPLICIT REAL* 8 (A-H,O-Z)
c     IMPLICIT REAL*8(a-h,o-z)                                          
C
      DIMENSION f(-4:*),a(0:10)
C
      CALL POL11I(f,n,a)
C     -----------
C
      POLY11 = a(10)*x
C
      DO i=9,1,-1
         POLY11 = (POLY11 + a(i))*x
      END DO
C
      POLY11=POLY11 + a(0)
C
      RETURN
C
      END
C---------------------------------------------------------------------
      SUBROUTINE POL11I(z,n,a)
C---------------------------------------------------------------------
C
C     Fitting an 11-point polynomial centred on the n-th value of z
C
      IMPLICIT REAL* 8 (A-H,O-Z)
c     IMPLICIT REAL*8(a-h,o-z)                                          RAD00260
C
      DIMENSION z(-4:*),a(0:10),d(5),C0(5),C1(5,5),C2(5,5),D1(5),D2(5)
C
      DATA C0/-73766.D0,192654.D0,-12276.D0,462.D0,-252.D0/,
     .     D1/1260.D0,181440.D0,34560.D0,120960.D0,725760.D0/,
     .     D2/50400.D0,362880.D0,172800.D0,120960.D0,3628800.D0/,
     .     ZERO/0.D0/
      DATA C1/1050.D0,-300.D0,75.D0,-12.5D0,1.D0,
     .        -70098.D0,52428.D0,-14607.D0,2522.D0,-205.D0,
     .        1938.D0,-1872.D0,783.D0,-152.D0,13.D0,
     .        -378.D0,408.D0,-207.D0,52.D0,-5.D0,
     .        42.D0,-48.D0,27.D0,-8.D0,1.D0/
      DATA C2/42000.D0,-6000.D0,1000.D0,-125.D0,8.D0,
     .        -140196.D0,52428.D0,-9738.D0,1261.D0,-82.D0,
     .        9690.D0,-4680.D0,1305.D0,-190.D0,13.D0,
     .        -378.D0,204.D0,-69.D0,13.D0,-1.D0,
     .        210.D0,-120.D0,45.D0,-10.D0,1.D0/
C
      a(0)=z(n)
C
      DO j=1,5
         d(j) = z(n+j) - z(n-j)
      END DO
C
      DO i=1,5
         a(2*i-1) = ZERO
         DO j=1,5
            a(2*i-1) = a(2*i-1) + C1(j,i)*d(j)
         END DO
         a(2*i-1) = a(2*i-1)/D1(i)
      END DO
C
      DO j=1,5
         d(j) = z(n+j) + z(n-j)
      END DO
C
      DO i=1,5
         a(2*i) = C0(i) * a(0)
         DO j=1,5
            a(2*i) = a(2*i) + C2(j,i)*d(j)
         END DO
         a(2*i) = a(2*i)/D2(i)
      END DO
C
      RETURN
 
      END
C---------------------------------------------------------------------       SUBROUTINE ludcmp(a,n,np,indx,d)
      SUBROUTINE ludcTA(a,n,np,indx,d)
C---------------------------------------------------------------------       SUBROUTINE ludcmp(a,n,np,indx,d)
C
      PARAMETER (NMAX=500)
C
      IMPLICIT REAL*8 (A-H,O-Z)
c     IMPLICIT REAL*8 (a-h,o-z)                                          
C
      DIMENSION a(np,np),vv(NMAX),indx(n)
C
      DATA ZERO/0.D0/,ONE/1.D0/,TINY/1.0D-25/
C
      d = ONE
C
      DO i=1,n
        aamax = ZERO
        DO j=1,n
          if (DABS(a(i,j)).GT.aamax) aamax = DABS(a(i,j))
        END DO
        if (aamax.EQ.ZERO) pause 'singular matrix in ludcmp'
        vv(i) = ONE/aamax
      END DO
C
      DO j=1,n
        DO i=1,j-1
          sum = a(i,j)
          DO k=1,i-1
            sum = sum - a(i,k)*a(k,j)
          END DO
          a(i,j) = sum
        END DO
C
        aamax = ZERO
        DO i=j,n
          sum = a(i,j)
          DO k=1,j-1
            sum = sum - a(i,k)*a(k,j)
          END DO
C
          a(i,j) = sum
          dum = vv(i)*DABS(sum)
          IF (dum.GE.aamax) THEN
            imax  = i
            aamax = dum
          END IF
        END DO

        IF (j.NE.imax) THEN
          DO k=1,n
            dum       = a(imax,k)
            a(imax,k) = a(j,k)
            a(j,k)    = dum
          END DO
          d        = -d
          vv(imax) = vv(j)
        END IF

        indx(j) = imax

        if(a(j,j).LT.TINY) a(j,j) = TINY

        IF (j.NE.n) THEN
          dum = ONE/a(j,j)
          DO i=j+1,n
            a(i,j) = a(i,j)*dum
          END DO
        END IF
      END DO

      RETURN
      END
C  (C) Copr. 1986-92 Numerical Recipes Software 1B=K'VIka5k.
C---------------------------------------------------------------------       SUBROUTINE ludcmp(a,n,np,indx,d)
      SUBROUTINE lubkTA(a,n,np,indx,b)
C---------------------------------------------------------------------       SUBROUTINE ludcmp(a,n,np,indx,d)
C
      IMPLICIT REAL*8 (A-H,O-Z)
c     IMPLICIT REAL*8 (a-h,o-z)                                          
C
      DIMENSION indx(n),a(np,np),b(n)
C
      DATA ZERO/0.D0/
C      
      ii = 0
      DO i=1,n
        ll    = indx(i)
        sum   = b(ll)
        b(ll) = b(i)
        IF (ii.NE.0) THEN
          DO j=ii,i-1
            sum = sum - a(i,j)*b(j)
          END DO
        ELSE IF (sum.NE.ZERO) THEN
          ii=i
        END IF	
        b(i) = sum
      END DO
      
      DO i=n,1,-1
        sum = b(i)
        DO j=i+1,n
          sum = sum - a(i,j)*b(j)
        END DO
        b(i) = sum/a(i,i)
      END DO
      
      RETURN
      END
C  (C) Copr. 1986-92 Numerical Recipes Software 1B=K'VIka5k.