      SUBROUTINE NEWMODES(XI,XN,IBMOD,IBI,IPP,NEWMDS,IGENMDS,           &
     &    IMODE,VELH,ZDISP,IFLAG,IFILEDLL,FILENDLL,IERROR,IWAMLOG,IER)
!DEC$ ATTRIBUTES DLLEXPORT :: NEWMODES
!-----------------------------------------------------------------------
!
!     COPYRIGHT (C) 2000-2008      WAMIT INCORPORATED
!
!-----------------------------------------------------------------------
!
!     Version : 6.4
!
!-----------------------------------------------------------------------
!
!     Source-code file : NEWMODES
!
!-----------------------------------------------------------------------
!
!     DESCRIPTION :  This routine calls subroutines prepared by
!     users which evaluate the vertical displacement and
!     the normal component of the user specified generalized modes.
!
!     NEWMODES is included in a DLL file to facilitate modifications by
!     users.  For information on the use of the DLL file, and generalized
!     modes, refer to the WAMIT User Manual.
!
!     This standard release version of NEWMODES includes subroutines
!     intended to define the generalized modes for the standard Test Runs
!     TEST08/18, TEST16/16a, TEST17a/17b and TEST23.  Refer to the User Manual
!     (appendix) for further information.
!
!     Two different types of subroutines are illustrated.  The generic
!     subroutine DEFINE is of the same form as in the DEFMOD program.
!     Different versions of DEFINE for other applications can be used
!     either in this subroutine or in DEFMOD, or both, as is done with
!     TEST08 and TEST18.  The other subroutines are more specific, but
!     provide a more logical and efficient scheme for use in NEWMODES.
!
!     NEWMODES is only called if IF (NEWMDS(IB)>0)
!-----------------------------------------------------------------------
!
!     Input Arguments (in alphabetical order):
!
!     Symbolic  Where value  Abbreviated Description (see routine where
!       name    is assigned    value is assigned for full description)
!     ------------------------------------------------------------------
!     IBMOD     RDPOT        Starting mode index for each body minus 1
!     IBI                    Body index of the panel or patch
!     IGENMDS                User-assigned integer>0 (optional)
!     IPP                    Local index of panel or patch:
!                               ILOWHI=0: IPP=local panel index on body
!                               ILOWHI=1: IPP=local patch index
!     IFLAG                  IFLAG= -1 return IMODE only
!     NEWMDS    RDPOT        Number of generalized modes
!     XI        BODYGEOM     Cartesian coordinate in body coordinates
!     XN        BODYGEOM     Normal vector
!
!
!     Output Arguments (in alphabetical order):
!
!     Symbolic    Description
!       name
!     ------------------------------------------------------------------
!     IMODE       Mode indices chosen from the rigid body modes 1,...,6
!                   for each generalized mode to have same symmetry
!                   (used only for reflection about planes of symmetry
!                    of the body)
!     VELH        Generalized normal for the radiation problem
!                    (denoted VEL/VEH for ILOWHI=0/1 in other files)
!     ZDISP       Vertical displacement of the surface
!
!  NOTE: The arrays VELH and ZDISP are initialized to zero
!        before calls to NEWMODES.  Only nonzero elements
!        of these arrays must be assigned.
!
!   MODS:
!      10/04 Subroutines added for general use or guidance:
!		   NAVIS_LIDS, GAP_FS (see comments below for explanations)
!       6/06 IFILEDLL,IERROR,IWAMLOG,IER added for use in subroutines
!                 which open/read files
!      10/06 IPP changed to local patch index in call from BODYGEOM
!       5/07 WAVEMAKER extended for piston, tent-function modes
!       5/08 FILENDLL added for optional input to subroutines
!           (FILENDLL is filename of GDF file for body IBI without `.gdf')
!-----------------------------------------------------------------------
      IMPLICIT NONE
      REAL XI(3),XN(3),VELH(*),ZDISP(*)
      INTEGER IBMOD(*),NEWMDS(*),IMODE(*),IBI,IPP,IGENMDS,IFLAG,IER,    &
     &        IFILEDLL,IERROR,IWAMLOG  
      CHARACTER*16 FILENDLL   
!-----------------------------------------------------------------------
!  U,V are local arrays, used only with DEFINE
!-----------------------------------------------------------------------
      REAL U(6+NEWMDS(IBI)),V(6+NEWMDS(IBI))
      INTEGER I1,I2
      IF(IGENMDS == 16) THEN
!-----------------------------------------------------------------------
!  Special subroutine for free-free bending modes along the horizontal
!  axis of a body, used in TEST16 (elastic modes of rectangular barge)
!  The patch index is not needed since these generalized modes are
!     nonzero on all patches of the body
!-----------------------------------------------------------------------
        CALL FREEBEAM_X                                                 &
     &        (XI,XN,IBMOD,IBI,NEWMDS,IMODE,VELH,ZDISP,IFLAG)
      ELSEIF(IGENMDS == 17) THEN
!-----------------------------------------------------------------------
!  Special subroutine for the lid of a cylinder with moonpool,
!  used in TEST17.  The patch index IPP is used to identify the lid.
!-----------------------------------------------------------------------
        CALL MOONPOOL_FS(XI,XN,IBMOD,IBI,IPP,NEWMDS,                    &
     &                                          IMODE,VELH,ZDISP,IFLAG)
      ELSEIF(IGENMDS == 18) THEN
!-----------------------------------------------------------------------
!  This code can be used with the generic subroutine DEFINE in the same
!  format as DEFINE is used in DEFMOD if ILOWHI=0 and IGENMDS=0
!  The version of DEFINE included in this file is used for TEST18
!-----------------------------------------------------------------------
        CALL DEFINE(XI(1),XI(2),XI(3),U,V,ZDISP,IMODE,NEWMDS(IBI))
        I1=IBMOD(IBI)+7
        I2=IBMOD(IBI)+6+NEWMDS(IBI)
        VELH(I1:I2)=XN(1)*U(I1:I2)+XN(2)*V(I1:I2)+XN(3)*ZDISP(I1:I2)
      ELSE IF(IGENMDS == 19) THEN
!-----------------------------------------------------------------------
!  Subroutine for three moonpools of drillship NAVIS EXPLORER I
!  See comments in subroutine, also the paper ``Integation of Geometry
!    Definition and Wave Analyis Software,'' by C.-H. Lee et al, 
!    Proceedings of 21st OMAE Conference, Oslo, Norway, (2002). 
!-----------------------------------------------------------------------
        CALL NAVIS_LIDS(XI,XN,IBMOD,IBI,IPP,NEWMDS,                     &
     &                                          IMODE,VELH,ZDISP,IFLAG)
      ELSEIF (IGENMDS == 20) THEN
!-----------------------------------------------------------------------
!  Special subroutine for the lid of a gap between two hulls
!  Adapted from CAT_FS routine in V62.  IGENMDS =  gap width
!  The body index IBI=3 for the lid 
!-----------------------------------------------------------------------
        CALL GAP_FS(XI,XN,IBMOD,IBI,IPP,NEWMDS,IGENMDS,                 &
     &                                          IMODE,VELH,ZDISP,IFLAG)
      ELSEIF ((IGENMDS==21).OR.(IGENMDS==211).OR.(IGENMDS==212)) THEN
!-----------------------------------------------------------------------
!  Subroutine for a bank of segmented wavemakers with hinges (21), 
!  piston modes (211) or tent-funtion modes (212) 
!  All modes are symmetric about x=0 and y=0 
!  IER is included since this subroutine opens/reads a special data file
!-----------------------------------------------------------------------
        CALL WAVEMAKER(XI,XN,IBMOD,IBI,IPP,NEWMDS,IGENMDS,              &
     &     IMODE,VELH,ZDISP,IFLAG,IFILEDLL,FILENDLL,IERROR,IWAMLOG,IER)
      ELSEIF (IGENMDS==22) THEN
!-----------------------------------------------------------------------
!  Subroutine for a segmented vessel with vertical hinge modes
!  IER is included since this subroutine opens/reads a special data file
!-----------------------------------------------------------------------
        CALL HINGE_MODES(XI,XN,IBMOD,IBI,IPP,NEWMDS,                    &
     &     IMODE,VELH,ZDISP,IFLAG,IFILEDLL,FILENDLL,IERROR,IWAMLOG,IER)
      ELSEIF (IGENMDS == 23) THEN
!-----------------------------------------------------------------------
!  Special subroutine for the lid of a gap between two hulls
!  Adopeted from GAP_FS but Fourier modes are used in both longitudinal
!  and lateral direction  
!  The body index IBI=3 for the lid 
!-----------------------------------------------------------------------
        CALL GAP_FS_FOURIER(XI,XN,IBMOD,IBI,IPP,NEWMDS,IGENMDS,         &
     &                                          IMODE,VELH,ZDISP,IFLAG)
      ELSEIF (IGENMDS == 24) THEN
!-----------------------------------------------------------------------
!   Hinged 2D wavemaker with two planes of symmetry
!-----------------------------------------------------------------------
        CALL WAVEMKR_SYM2D(XI,XN,IBMOD,IBI,IPP,NEWMDS,                  &
     &     IMODE,VELH,ZDISP,IFLAG,IFILEDLL,FILENDLL,IERROR,IWAMLOG,IER)
      ELSEIF (IGENMDS == 25) THEN
!-----------------------------------------------------------------------
!   Bank of hinged 3D wavemakers with two planes of symmetry  
!-----------------------------------------------------------------------
        CALL WAVEMKR_SYM3D(XI,XN,IBMOD,IBI,IPP,NEWMDS,                  &
     &     IMODE,VELH,ZDISP,IFLAG,IFILEDLL,FILENDLL,IERROR,IWAMLOG,IER)
      ENDIF
      RETURN
      END

      SUBROUTINE FREEBEAM_X                                             &
     &        (XI,XN,IBMOD,IBI,NEWMDS,IMODE,VELH,ZDISP,IFLAG)
!-----------------------------------------------------------------------
!
!     COPYRIGHT (C) 2000      WAMIT INCORPORATED
!
!-----------------------------------------------------------------------
!
!     Version : 6.0
!
!-----------------------------------------------------------------------
!
!     Source-code file : NEWMODES
!
!-----------------------------------------------------------------------
!
!     DESCRIPTION :  The vertical displacement and the normal component
!     of free-free beam modes along the x-axis (-XL < X < XL).  XL is
!     assigned in a PARAMETER statement.
!
!-----------------------------------------------------------------------
!
!     Input Arguments (in alphabetical order):
!
!     Symbolic  Where value  Abbreviated Description (see routine where
!       name    is assigned    value is assigned for full description)
!     ------------------------------------------------------------------
!     IBMOD     GENMDS       Starting mode index for each body minus 1
!     IBI                    Body index of the patch
!     IFLAG                  IFLAG= -1 return IMODE only
!                            IFLAG= 1 return Normal only
!     NEWMDS    RDPOT        Number of generalized modes
!     XI        BODYGEOM     Cartesian coordinate in body coordinates
!     XN        BODYGEOM     Normal vector
!
!
!     Output Arguments (in alphabetical order):
!
!     Symbolic    Description
!       name
!     ------------------------------------------------------------------
!     IMODE       Mode indices chosen from the rigid body modes 1,...,6
!                   for each generalized mode to have same symmetry
!     VELH        Generalized normal for the radiation problem
!     ZDISP       Vertical displacement of the surface
!
!-----------------------------------------------------------------------
      IMPLICIT NONE
      REAL XI(*),XN(*),VELH(*),ZDISP(*)
      INTEGER IBMOD(*),NEWMDS(*),IMODE(*),IBI,IFLAG
!-----------------------------------------------------------------------
!   Eigenvalues of the first 10 bending modes of a free-free beam
!-----------------------------------------------------------------------
      REAL, PARAMETER :: WD(1:10) =                                     &
     &(/ 2.36502,   3.92660,   5.49780,   7.06858,   8.63938,           &
     &  10.21018,  11.78097,  13.35177,  14.92257,  16.49336 /)
!-----------------------------------------------------------------------
!   non-dimensionalize x-coordinate relative to half-length of barge
!-----------------------------------------------------------------------
      REAL, PARAMETER :: XL = 80.0
      REAL XI1
      INTEGER I,IB6
      XI1=2.*XI(1)/XL
      IF(IFLAG.EQ.-1) THEN
!-----------------------------------------------------------------------
!     Assign mode symmetry indices for symmetric/antisymmetric modes
!     These are the same as heave/pitch (modes 3/1), respectively
!-----------------------------------------------------------------------
         DO I=1,NEWMDS(IBI),2
            IMODE(IBMOD(IBI)+I+6)=3
         ENDDO
         DO I=2,NEWMDS(IBI),2
            IMODE(IBMOD(IBI)+I+6)=1
         ENDDO
      ELSE
!-----------------------------------------------------------------------
!     Symmetric modes
!-----------------------------------------------------------------------
         DO I=1,NEWMDS(IBI),2
            IB6=IBMOD(IBI)+I+6
            ZDISP(IB6)=0.5*(COS(WD(I)*XI1)/COS(WD(I))+                  &
     &                               COSH(WD(I)*XI1)/COSH(WD(I)))
            VELH(IB6)=XN(3)*ZDISP(IB6)
         ENDDO
!-----------------------------------------------------------------------
!     Antisymmetric modes
!-----------------------------------------------------------------------
         DO I=2,NEWMDS(IBI),2
            IB6=IBMOD(IBI)+I+6
            ZDISP(IB6)=0.5*(SIN(WD(I)*XI1)/SIN(WD(I))+                  &
     &                            SINH(WD(I)*XI1)/SINH(WD(I)))
            VELH(IB6)=XN(3)*ZDISP(IB6)
         ENDDO
      ENDIF
      RETURN
      END

      SUBROUTINE DEFINE(X,Y,Z,U,V,W,ISYM,NEWMDS)
!-----------------------------------------------------------------------
!
!     COPYRIGHT (C) 1998 MASSACHUSETTS INSTITUTE OF TECHNOLOGY
!     COPYRIGHT (C) 2000 WAMIT, Inc.
!
!-----------------------------------------------------------------------
!
!     Version : 6.0
!
!-----------------------------------------------------------------------
!
!     Source-code files : NEWMODES and DEFMOD
!
!    DEFMOD is a separate utility intended for use when ILOWHI=0 and
!    IGENMDS=0  (see WAMIT User Manual, Chapter 8)
!-----------------------------------------------------------------------
!
!     DESCRIPTION : Evaluate new mode vectors U,V,W for specified
!                   position (X,Y,Z)
!                   (VEL = normal velocity in mode J = U*nx+V*ny+W*nz)
!
!                   This example of DEFINE defines four Jacobi
!                   bending modes for a vertical column
!                   (for use with Test Run 08/18)
!
!-----------------------------------------------------------------------
!
!     THIS SUBROUTINE MAY BE MODIFIED BY THE USER FOR USER-DEFINED
!             MODE SHAPES
!
!-----------------------------------------------------------------------
!
!     Input Arguments (in alphabetical order):
!
!     Symbolic  Where value  Description
!       name    is assigned
!     ------------------------------------------------------------------
!     X,Y,Z     DEFMOD      Coordinates of panel centroid
!                           (nondimensionalized by ULEN)
!
!     Output Arguments (in alphabetical order):
!
!     Symbolic    Description
!       name
!     ------------------------------------------------------------------
!     ISYM      Symmetry index of new modes
!     NEWMDS    Number of new modes
!     U(J)      X-component of new modes (J=7,8,...,6+NEWMDS)
!     V(J)      Y-component of new modes (J=7,8,...,6+NEWMDS)
!     W(J)      Z-component of new modes (J=7,8,...,6+NEWMDS)
!-----------------------------------------------------------------------
      IMPLICIT NONE
      REAL X,Y,Z,U(*),V(*),W(*)
      INTEGER ISYM(*)
!-----------------------------------------------------------------------
!   User-defined code follows:
!
!   The following example defines NEWMDS new modes corresponding to
!   bending modes of the vertical beam  (Shifted Jacobi polynomials).
!   Q is the normalized vertical coordinate for a depth of 200 meters,
!   (Q=0 at bottom, 1 at free surface).
!-----------------------------------------------------------------------
      REAL Q,DENA,A,DENBC,B,C
      INTEGER NEWMDS,J,N
      NEWMDS=4
      Q=1.+Z/200.
!-----------------------------------------------------------------------
!   First loop over all new modes initializes symmetry index and V,W
!   Modes have same symmetries as surge
!   Velocities V,W are zero for all modes
!-----------------------------------------------------------------------
      DO J=7,NEWMDS+6
        ISYM(J)=1
        V(J)=0.0E0
        W(J)=0.0E0
      END DO
!-----------------------------------------------------------------------
!   Evaluate first two polynomials explicitly, others by recursion
!-----------------------------------------------------------------------
      U(7) =Q*Q
      U(8) =U(7)*(6.*Q-5.)
      DO J=9,NEWMDS+6
        N=J-6
        DENA=(N-1)*(N+3)
        A=2*(N+1)*(N+N+1)/DENA
        DENBC=N*(N-1)*(N+3)
        B=-(N+N+1)*(N*N+N+4)/DENBC
        C=-(N-2)*(N+1)*(N+2)/DENBC
        U(J)=(A*Q+B)*U(J-1)+C*U(J-2)
      END DO
      RETURN
      END

      SUBROUTINE MOONPOOL_FS                                            &
     &        (XI,XN,IBMOD,IBI,IPP,NEWMDS,IMODE,VELH,ZDISP,IFLAG)
!-----------------------------------------------------------------------
!
!     COPYRIGHT (C) 2000 WAMIT INC.
!
!-----------------------------------------------------------------------
!
!     Version : 6.0
!
!-----------------------------------------------------------------------
!
!     Source-code file : NEWMODES
!
!-----------------------------------------------------------------------
!
!     DESCRIPTION :  The vertical displacement and the normal component
!     of the pumping (heave) and sloshing (pitch) modes for a moonpool
!     lid defined by patch number 4, used in TEST17a,b
!
!     The geometry is a circular cylinder with moonpool generated by
!     subroutine CYLMP in GEOMXACT (IGDEF=-7).  There are 3 patches on
!     the outer, bottom, inner surfaces of the cylinder and patch 4 is
!     included to define the lid when NPATCH=4 is input.  For incident
!     waves from oblique directions a transverse sloshing mode should be
!     added.  For more accurate representation of these modes when the
!     moonpool radius is large (relative to the wavelength),
!     Bessel function modes can be used.
!
!-----------------------------------------------------------------------
!
!     Input Arguments (in alphabetical order):
!
!     Symbolic  Where value  Abbreviated Description (see routine where
!       name    is assigned    value is assigned for full description)
!     ------------------------------------------------------------------
!     IBMOD     GENMDS       Starting mode index for each body minus 1
!     IBI                    Body index
!     IPP                    Patch index
!     IFLAG                  IFLAG= -1 return IMODE only
!     NEWMDS    RDPOT        Number of generalized modes
!     XI        BODYGEOM     Cartesian coordinate in body coordinates
!     XN        BODYGEOM     Normal vector
!
!
!     Output Arguments (in alphabetical order):
!
!     Symbolic    Description
!       name
!     ------------------------------------------------------------------
!     IMODE       Mode indices chosen from the rigid body modes 1,...,6
!                   for each generalized mode to have same symmetry
!     VELH        Generalized normal for the radiation problem
!     ZDISP       Vertical displacement of the surface
!
!-----------------------------------------------------------------------
      IMPLICIT NONE
      REAL XI(3),XN(3),VELH(*),ZDISP(*)
      INTEGER IBMOD(*),NEWMDS(*),IMODE(*),IBI,IPP,IFLAG
      IF(IFLAG.EQ.-1) THEN
!-----------------------------------------------------------------------
!   Assign mode symmetry indices
!-----------------------------------------------------------------------
         IMODE(7)=3
         IMODE(8)=1
      ELSEIF (IPP == 4) THEN
!-----------------------------------------------------------------------
!     Assign two modes, heave and pitch (pumping mode and approximation
!       for 1st longitudinal sloshing mode)
!-----------------------------------------------------------------------
         ZDISP(7)=1.0
         ZDISP(8)=XI(1)
         VELH(7:8)=ZDISP(7:8)
      ENDIF
 99   RETURN
      END

      SUBROUTINE NAVIS_LIDS                                             &
     &        (XI,XN,IBMOD,IBI,IPP,NEWMDS,IMODE,VELH,ZDISP,IFLAG)
!-----------------------------------------------------------------------
!
!     COPYRIGHT (C) 2001 WAMIT INC.
!
!-----------------------------------------------------------------------
!
!     Version : 6.1
!
!-----------------------------------------------------------------------
!
!     Source-code file : NEWMODES
!
!-----------------------------------------------------------------------
!
!     DESCRIPTION :  The vertical displacement and the normal component
!     of the pumping (heave) modes for 3 moonpools on Navis Explorer
!     defined by patch numbers 56.57.58
!
!-----------------------------------------------------------------------
!
!     Input Arguments (in alphabetical order):
!
!     Symbolic  Where value  Abbreviated Description (see routine where
!       name    is assigned    value is assigned for full description)
!     ------------------------------------------------------------------
!     IBMOD     GENMDS       Starting mode index for each body minus 1
!     IBI                    Body index
!     IPP                    Patch index
!     IFLAG                  IFLAG= -1 return IMODE only
!     NEWMDS    RDPOT        Number of generalized modes
!     XI        BODYGEOM     Cartesian coordinate in body coordinates
!     XN        BODYGEOM     Normal vector
!
!
!     Output Arguments (in alphabetical order):
!
!     Symbolic    Description
!       name
!     ------------------------------------------------------------------
!     IMODE       Mode indices chosen from the rigid body modes 1,...,6
!                   for each generalized mode to have same symmetry
!     VELH        Generalized normal for the radiation problem
!     ZDISP       Vertical displacement of the surface
!
!-----------------------------------------------------------------------
      IMPLICIT NONE
      REAL XI(3),XN(3),VELH(*),ZDISP(*)
      INTEGER IBMOD(*),NEWMDS(*),IMODE(*),IBI,IPP,IFLAG
      IF(IFLAG.EQ.-1) THEN
!-----------------------------------------------------------------------
!   Assign mode symmetry indices
!-----------------------------------------------------------------------
         IMODE(7:9)=3
      ELSEIF (IPP == 56) THEN
         VELH(9)=1.0
      ELSEIF (IPP == 57) THEN
         VELH(8)=1.0
      ELSEIF (IPP == 58) THEN
         VELH(7)=1.0
      ENDIF
      ZDISP(7:9)=VELH(7:9)
 99   RETURN
      END

      SUBROUTINE GAP_FS                                                 &
     &      (XI,XN,IBMOD,IBI,IPP,NEWMDS,IGENMDS,IMODE,VELH,ZDISP,IFLAG)
!-----------------------------------------------------------------------
!
!     COPYRIGHT (C) 2004 WAMIT INC.
!
!-----------------------------------------------------------------------
!
!     Version : 6.2
!
!-----------------------------------------------------------------------
!
!     Source-code file : NEWMODES
!
!-----------------------------------------------------------------------
!
!     DESCRIPTION :  Adapted from CAT_FS  
!  Modes for lid on free surface between hulls of two independent hulls
!  Generalized modes represent vertical motions of the gap free surface:
!  Longitudinal Chebyshev modes T_n(x')= (1, cos t, cos 2t, ...)
!        cos t = x',   x' normalized between -1,1 in gap
!  Lateral modes (1, sin(y'), cos(2y')   y' normalized between +/- pi/2
!  Domain of the gap:  X1<X<X2, -HGAP < Y < HGAP
!  IBI=3 is the body index of the lid, modes 19,20,...,(18+NEWMDS(3))
!   NEWMDS(IBI) is arbitrary but should be Nfourier*Nchebyshev
!
!   NFOURIER upto 3 can be specified. Larger Nfourier, additional fourier
!   modes must be added in this suroutine, sin(3y')..etc.
!
!   Mode sequence is Ncheb * 1, Ncheb * sin(y'), Ncheb * cos(2y')
!  
!     (X1=-80, X2=+80) and HGAP=4 are assigned in subroutine
!
!   User must adjust paramters X1,X2,HGAP,Nfourier below 
!   X1/X2 lower/upper limit in x in body coordinate
!   HGAP  half of the gap
!   Nfouier number of fourier modes in lateral direction
!   (Nchebyshev is determine by newmds/Nfourier)
!-----------------------------------------------------------------------
!
!     Input Arguments (in alphabetical order):
!
!     Symbolic  Where value  Abbreviated Description (see routine where
!       name    is assigned    value is assigned for full description)
!     ------------------------------------------------------------------
!     IBMOD     GENMDS       Starting mode index for each body minus 1
!     IBI                    Body index
!     IPP                    Patch index
!     IFLAG                  IFLAG= -1 return IMODE only
!     NEWMDS    RDPOT        Number of generalized modes
!     XI        BODYGEOM     Cartesian coordinate in body coordinates
!     XN        BODYGEOM     Normal vector
!
!
!     Output Arguments (in alphabetical order):
!
!     Symbolic    Description
!       name
!     ------------------------------------------------------------------
!     IMODE       Mode indices chosen from the rigid body modes 1,...,6
!                   for each generalized mode to have same symmetry
!     VELH        Generalized normal for the radiation problem
!     ZDISP       Vertical displacement of the surface
!
!-----------------------------------------------------------------------
      IMPLICIT NONE
      REAL XI(3),XN(3),VELH(*),ZDISP(*)
      INTEGER IBMOD(*),NEWMDS(*),IMODE(*),IBI,IPP,IFLAG
      INTEGER J,NCHEB,IGENMDS
      INTEGER, SAVE :: J1,J2,J3,J4,J5,J6
      REAL T,TH,FY
      REAL, PARAMETER :: PI2=1.570796
      REAL, SAVE :: FACY
      INTEGER, PARAMETER :: NFOURIER=2
      REAL, PARAMETER :: X1=-80,X2=80,HGAP=4.
      IF(IFLAG==-1) THEN  
         FACY=PI2/HGAP
         J1=19
         NCHEB=NEWMDS(IBI)/NFOURIER
         J2=J1-1+NCHEB
!-----------------------------------------------------------------------
!   Assign mode symmetry indices (arbitrary, since ISX=ISY=0)
!-----------------------------------------------------------------------
         IMODE(J1:J2)=3
         IF(NFOURIER >=2) THEN
           J3=J2+1
           J4=J2+NCHEB
           IMODE(J3:J4)=3
         ENDIF
         IF(NFOURIER >=3) THEN
           J5=J4+1
           J6=J4+NCHEB
           IMODE(J5:J6)=3
         ENDIF
!-----------------------------------------------------------------------
!     Assign Chebyshev modes using recursion relation
!-----------------------------------------------------------------------
      ELSEIF (IBI == 3) THEN
         T=(2.*XI(1)-X1-X2)/(X2-X1)
!-----------------------------------------------------------------------
!     Assign Chebyshev modes using recursion relation
!     Modes (J1,J2) are for f(y)=constant
!-----------------------------------------------------------------------
         ZDISP(J1)=1.0
         IF (J2>J1) ZDISP(J1+1)=T
         DO J=J1+2,J2
                ZDISP(J)=2.*T*ZDISP(J-1)-ZDISP(J-2)
         ENDDO
         VELH(J1:J2)=ZDISP(J1:J2)
!-----------------------------------------------------------------------
!     Transverse sine modes
!-----------------------------------------------------------------------
         IF(NFOURIER >=2) THEN
           TH=XI(2)*FACY
           FY=SIN(TH)
           ZDISP(J3:J4)=ZDISP(J1:J2)*FY
           VELH(J3:J4)=ZDISP(J3:J4)
         ENDIF
         IF(NFOURIER >=3) THEN
           FY=COS(2.*TH)
           ZDISP(J5:J6)=ZDISP(J1:J2)*FY
           VELH(J5:J6)=ZDISP(J5:J6)
         ENDIF
      ENDIF
 99   RETURN
      END

       SUBROUTINE WAVEMKR_SYM2D(XI,XN,IBMOD,IBI,IPP,NEWMDS,             &
     &     IMODE,VELH,ZDISP,IFLAG,IFILEDLL,FILENDLL,IERROR,IWAMLOG,IER)
!-----------------------------------------------------------------------
!
!     COPYRIGHT (C) 2006 WAMIT INC.
!
!-----------------------------------------------------------------------
!
!     Version : 6.4
!
!-----------------------------------------------------------------------
!
!     Source-code file : NEWMODES
!
!-----------------------------------------------------------------------
!
!     DESCRIPTION :  Flap wavemaker with hinge at depth ZHINGE on patch 1
!                    This subroutine is for two planes of symmetry, 
!                    symmetric and antisymmetric modes for usein rectangular
!                    wave tank to generate/absorb waves
!                    NEWMDS=2 (symmetric and antisymmetric modes about x=0
!                            no antisymmetric modes about y=0)
!                    Modes are normalized with unit amplitude at Z=0
!     NB:  1)  additional patches
!              can be used for sides of tank and ends below wavemaker
!          2) only one patch is permitted on each wavemaker segment
!              Low-order method can be used, but only with one panel
!              per wavemaker.  To remove these restrictions, modify the
!              subroutine so that the index IPP has the correct relation
!              to each segment so that there is one generalized mode for
!              each segment
!
!   Error return (IER=1) if WAVEMAKER_DEPTH.DAT cannot be opened and read 
!
!  N.B. This subroutine has access to three external files, with unit
!        numbers IFILEDLL, IWAMLOG, IERROR, which can be used as follows:
!    
!      IFILEDLL:  input data file 
!      IWAMLOG:   WAMITLOG.TXT (log file for WAMIT run)
!      IERROR:    ERRORP.LOG or ERRORF.LOG (error files for Poten/Force)
!
!-----------------------------------------------------------------------
!
!     Input Arguments (in alphabetical order):
!
!     Symbolic  Where value  Abbreviated Description (see routine where
!       name    is assigned    value is assigned for full description)
!     ------------------------------------------------------------------
!     IBMOD     GENMDS       Starting mode index for each body minus 1
!     IBI                    Body index
!     IPP                    Patch index
!     IFLAG                  IFLAG= -1 return IMODE only
!     NEWMDS    RDPOT        Number of generalized modes
!     XI        BODYGEOM     Cartesian coordinate in body coordinates
!     XN        BODYGEOM     Normal vector
!
!
!     Output Arguments (in alphabetical order):
!
!     Symbolic    Description
!       name
!     ------------------------------------------------------------------
!     IMODE       Mode indices chosen from the rigid body modes 1,...,6
!                   for each generalized mode to have same symmetry
!     VELH        Generalized normal for the radiation problem
!     ZDISP       Vertical displacement of the surface
!
!  MODS
!    6/08 optional input IPPMAX, set IPPMAX=1 to only assign modes on 
!          patch=1  for use with oval tank study  (TANKOV)
!         If NEWMDS=1 only one generalized mode assigned
!-----------------------------------------------------------------------
      IMPLICIT NONE
      REAL XI(3),XN(3),VELH(*),ZDISP(*)
      INTEGER IBMOD(*),NEWMDS(*),IMODE(*),IBI,IPP,IFLAG,IER,IFILEDLL,   &
     &        IERROR,IWAMLOG        
      REAL, SAVE :: ZHINGE
      INTEGER , SAVE :: IPPMAX
      CHARACTER*80 HEADDAT
      CHARACTER*16 FILENDLL
      CHARACTER*14 :: FILEN='_wmkrhinge.dat'
      IF(IFLAG.EQ.-1) THEN
!-----------------------------------------------------------------------
!   Derive filename of GDF file and add extension for data file
!   Open input file with data for hinge depth, copy file to wamitlog.txt
!     Read header line and ZHINGE
!-----------------------------------------------------------------------
        OPEN (UNIT=IFILEDLL,FILE=TRIM(FILENDLL)//FILEN,                 &
     &        STATUS='OLD',ERR=89)
        READ (IFILEDLL,900,ERR=89) HEADDAT
        READ (IFILEDLL,*,ERR=89) ZHINGE
        IPPMAX=3
        READ (IFILEDLL,*,END=10,ERR=89) IPPMAX
  10    REWIND (IFILEDLL)
!-----------------------------------------------------------------------
!   Copy inputs to log file
!-----------------------------------------------------------------------
        WRITE (IWAMLOG,910) FILEN
        DO  
          READ (IFILEDLL,900,END=20,ERR=20) HEADDAT
          WRITE (IWAMLOG,900) TRIM(HEADDAT)
        END DO
   20   CLOSE (IFILEDLL)
!-----------------------------------------------------------------------
!   Assign generalized-mode symmetry indices (same as surge/heave)
!-----------------------------------------------------------------------
        IMODE(IBMOD(IBI)+7)=1
        IF (NEWMDS(IBI)==2) IMODE(IBMOD(IBI)+8)=3
      ELSEIF (IPP==1) THEN
!-----------------------------------------------------------------------
!     Assign hinge mode, normalized =1 at free surface. 
!     Patch index and generalized mode index are the same 
!-----------------------------------------------------------------------
        VELH(IBMOD(IBI)+7)=1.0-XI(3)/ZHINGE
        IF (NEWMDS(IBI)==2) VELH(IBMOD(IBI)+8)=VELH(IBMOD(IBI)+7)
      ELSEIF ((IPP==3).AND.(IPP.LE.IPPMAX)) THEN
!-----------------------------------------------------------------------
!     Assign hinge mode, normalized =1 at free surface. 
!     Patch index and generalized mode index are the same 
!-----------------------------------------------------------------------
        VELH(IBMOD(IBI)+7)=-(1.0-XI(3)/ZHINGE)
        IF (NEWMDS(IBI)==2) VELH(IBMOD(IBI)+8)=-VELH(IBMOD(IBI)+7)
      ENDIF
      GOTO 99
!-----------------------------------------------------------------------
!   Write error statement to monitor and error file, then stop
!     Patch index and generalized mode index are the same 
!-----------------------------------------------------------------------
  89  WRITE (*,900) ' Error in input file '//FILEN
      WRITE (IERROR,900) ' Error in input file '//FILEN
      IER=1
  99  RETURN
 900  FORMAT(A)
 910  FORMAT (/,' Input file ',A)
      END

      SUBROUTINE WAVEMKR_SYM3D(XI,XN,IBMOD,IBI,IPP,NEWMDS,              &
     &     IMODE,VELH,ZDISP,IFLAG,IFILEDLL,FILENDLL,IERROR,IWAMLOG,IER)
!-----------------------------------------------------------------------
!
!     COPYRIGHT (C) 2006 WAMIT INC.
!
!-----------------------------------------------------------------------
!
!     Version : 6.4
!
!-----------------------------------------------------------------------
!
!     Source-code file : NEWMODES
!
!-----------------------------------------------------------------------
!
!     DESCRIPTION :  Flap wavemakers with hinge at depth ZHINGE  
!                    This subroutine is for two planes of symmetry, all
!                    symmetric and antisymmetric modes for use in rectangular
!                    or circular wave tank to generate/absorb waves
!                    NEWMDS=4*number of wavemakers in one quadrant
!                    
!                    Modes are normalized with unit amplitude at Z=0
!     NB:      only one patch is permitted on each wavemaker segment
!              Low-order method can be used, but only with one panel
!              per wavemaker.  To remove these restrictions, modify the
!              subroutine so that the index IPP has the correct relation
!              to each segment so that there is one generalized mode for
!              each segment
!
!   Error return (IER=1) if _wmkrhinge.dat cannot be opened and read 
!
!  N.B. This subroutine has access to three external files, with unit
!        numbers IFILEDLL, IWAMLOG, IERROR, which can be used as follows:
!    
!      IFILEDLL:  input data file  
!      IWAMLOG:   WAMITLOG.TXT (log file for WAMIT run)
!      IERROR:    ERRORP.LOG or ERRORF.LOG (error files for Poten/Force)
!
!-----------------------------------------------------------------------
!
!     Input Arguments (in alphabetical order):
!
!     Symbolic  Where value  Abbreviated Description (see routine where
!       name    is assigned    value is assigned for full description)
!     ------------------------------------------------------------------
!     IBMOD     GENMDS       Starting mode index for each body minus 1
!     IBI                    Body index
!     IPP                    Patch index
!     IFLAG                  IFLAG= -1 return IMODE only
!     NEWMDS    RDPOT        Number of generalized modes
!     XI        BODYGEOM     Cartesian coordinate in body coordinates
!     XN        BODYGEOM     Normal vector
!
!
!     Output Arguments (in alphabetical order):
!
!     Symbolic    Description
!       name
!     ------------------------------------------------------------------
!     IMODE       Mode indices chosen from the rigid body modes 1,...,6
!                   for each generalized mode to have same symmetry
!     VELH        Generalized normal for the radiation problem
!     ZDISP       Vertical displacement of the surface
!
!  MODS:
!    12/07 extended to permit NW wavemakers on one quadrant of tank
!          patches IPP=1,2,...,NW used for each wavemaker
!          NEWMDS=4*NW  for symmetries ISYM=1,2,3,6
!     5/08 filename changed to gdf filename + wmkrhinge.dat
!-----------------------------------------------------------------------
      IMPLICIT NONE
      REAL XI(3),XN(3),VELH(*),ZDISP(*)
      INTEGER IBMOD(*),NEWMDS(*),IMODE(*),IBI,IPP,IFLAG,IER,IFILEDLL,   &
     &        IERROR,IWAMLOG         
      REAL, SAVE :: ZHINGE
      INTEGER, SAVE :: NW,NMDS
      INTEGER  J,J1
      CHARACTER*80 HEADDAT
      CHARACTER*16 FILENDLL
      CHARACTER*14 :: FILEN='_wmkrhinge.dat'
      IF(IFLAG.EQ.-1) THEN
!-----------------------------------------------------------------------
!   Derive filename of GDF file and add extension for data file
!   Open input file with data for hinge depth, copy file to wamitlog.txt
!     Read header line and ZHINGE
!-----------------------------------------------------------------------
        OPEN (UNIT=IFILEDLL,FILE=TRIM(FILENDLL)//FILEN,                 &
     &        STATUS='OLD',ERR=89)
        READ (IFILEDLL,900,ERR=89) HEADDAT
        READ (IFILEDLL,*,ERR=89) ZHINGE
        REWIND (IFILEDLL)
!-----------------------------------------------------------------------
!   Copy inputs to log file
!-----------------------------------------------------------------------
        WRITE (IWAMLOG,910) FILEN
        DO  
          READ (IFILEDLL,900,END=10,ERR=10) HEADDAT
          WRITE (IWAMLOG,900) TRIM(HEADDAT)
        END DO
   10   CLOSE (IFILEDLL)
        NMDS=NEWMDS(IBI)
        NW=NMDS/4
        IF (4*NW<NMDS) GOTO 88
!-----------------------------------------------------------------------
!   Assign mode symmetry indices (same as surge/sway/heave/yaw)
!-----------------------------------------------------------------------
        DO J=1,NMDS,4
          IMODE(IBMOD(IBI)+J+6)=1
          IMODE(IBMOD(IBI)+J+7)=2
          IMODE(IBMOD(IBI)+J+8)=3
          IMODE(IBMOD(IBI)+J+9)=6
        ENDDO
      ELSE
        J1=IBMOD(IBI)+4*(IPP-1)+7
!-----------------------------------------------------------------------
!     Assign hinge mode, normalized =1 at free surface. 
!     Patch index and generalized mode index are the same 
!-----------------------------------------------------------------------
        VELH(J1:J1+3)=1.0-XI(3)/ZHINGE
      ENDIF
      GOTO 99
!-----------------------------------------------------------------------
!   Write error statement to monitor and error file, then stop
!     Patch index and generalized mode index are the same 
!-----------------------------------------------------------------------
  88  WRITE (*,900) ' NEWMDS < 4*number of wavemakers '//FILEN
      GOTO 90
!-----------------------------------------------------------------------
!   Write error statement to monitor and error file, then stop
!     Patch index and generalized mode index are the same 
!-----------------------------------------------------------------------
  89  WRITE (*,900) ' Error in input file '//FILEN
      WRITE (IERROR,900) ' Error in input file '//FILEN
  90  IER=1
  99  RETURN
 900  FORMAT(A)
 910  FORMAT (/,' Input file ',A)
      END

       SUBROUTINE WAVEMAKER(XI,XN,IBMOD,IBI,IPP,NEWMDS,IGENMDS,         &
     &     IMODE,VELH,ZDISP,IFLAG,IFILEDLL,FILENDLL,IERROR,IWAMLOG,IER)
!-----------------------------------------------------------------------
!
!     COPYRIGHT (C) 2005-2007 WAMIT INC.
!
!-----------------------------------------------------------------------
!
!     Version : 6.4
!
!-----------------------------------------------------------------------
!
!     Source-code file : NEWMODES
!
!-----------------------------------------------------------------------
!
!     DESCRIPTION :  Segmented wavemaker with either flap, piston, or 
!                           tent-function modes (depending on IGENMDS):
!
!          IGENMDS=21 hinged wavemakers (flap), hinged about Z=ZHINGE, 
!                      which must be specified in the file WAVEMAKER_HINGE.DAT
!                      Modes are normalized with unit amplitude at Z=0
!                     (same as V6.3 version of this subroutine, except
!                     filename is changed from WAVEMAKER_DEPTH.DAT in V6.3) 
! 
!                    The number of individual wavemakers is equal to 
!                       NEWMDS(IBI) and these correspond to patches
!                       1:NEWMDS(IBI).  
!
!          IGENMDS=211 piston wavemakers  file WAVEMAKER_PISTON.DAT is not required,
!                    but may be used to specify the optional parameter NPPW
!                    Modes are normalized with unit amplitude 
!
!          IGENMDS=212 tent-function modes, file WAVEMAKER_TENT.DAT is required
!                    to specify Y1=first side of first wavemaker segment and
!                    DELY=width of each segment.  Tent-function modes have unit
!                    amplitude at the joints between adjacent segments, zero at
!                    all other joints.  Modes are independent of depth on the
!                    patches.  Optional parameter NPPW may be specified.
!
!     Flap and piston wavemakers can be in the planes x=0 and/or y=0
!         (or parallel to these planes, if they protrude from the walls)
!     Tent-function wavemakers must be in or parrallel to the plane x=0
!        (mode shapes are piecewise linear functions of y)
!
!                    Modes are defined with symmetry about x=0 and y=0
!                    to simulate adjacent walls of a rectangular tank
!                    
!     NB:  1) patches must be assigned as above.  additional patches
!              can be used for sides and bottom of wavemaker, if it
!              protrudes from the wall X=0, provided these patches are
!              assigned with index > NEWMDS (or NEWMDS-1 for tent-functions)
!          2) only one patch is permitted on each wavemaker segment
!              Low-order method can be used, but only with one panel
!              per wavemaker.  To remove these restrictions, modify the
!              subroutine so that the index IPP has the correct relation
!              to each segment so that there is one generalized mode for
!              each segment
!
!   MOD 1/07:  
!  THE ABOVE RESTRICTIONS ARE PARTIALLY REMOVED BY INPUTTING NPPW>1
!     The parameter NPPW is optional and can be included on 
!     line 3 of the input file WAVEMAKER.DAT.  NPPW= number of panels
!     or patches per each generalized mode.  In the low-order method, set
!     NPPW equal to the number of panels on each wavemaker (must be the same)
!     Can also use NPPW to apply the same generalized mode to more than one
!     wavemaker, moving with the same amplitude and phase.  NPPW should
!     correspond to the number of panels per wavemaker in the GDF file.
!
!   Error return (IER=1) if input data file cannot be opened and read 
!
!  N.B. This subroutine has access to three external files, with unit
!        numbers IFILEDLL, IWAMLOG, IERROR, which can be used as follows:
!    
!      IFILEDLL:  input data file  
!      IWAMLOG:   WAMITLOG.TXT (log file for WAMIT run)
!      IERROR:    ERRORP.LOG or ERRORF.LOG (error files for Poten/Force)
!
!-----------------------------------------------------------------------
!
!     Input Arguments (in alphabetical order):
!
!     Symbolic  Where value  Abbreviated Description (see routine where
!       name    is assigned    value is assigned for full description)
!     ------------------------------------------------------------------
!     IBMOD     GENMDS       Starting mode index for each body minus 1
!     IBI                    Body index
!     IGENMDS   RDCONFIG     Index specifying the mode shapes:
!                            IGENMDS=21 hinged wavemakers (flap)
!                            IGENMDS=211 piston wavemakers
!                            IGENMDS=212 tent function modes
!     IPP                    Patch index
!     IFLAG                  IFLAG= -1 return IMODE only
!     NEWMDS    RDPOT        Number of generalized modes
!     XI        BODYGEOM     Cartesian coordinate in body coordinates
!     XN        BODYGEOM     Normal vector
!
!
!     Output Arguments (in alphabetical order):
!
!     Symbolic    Description
!       name
!     ------------------------------------------------------------------
!     IMODE       Mode indices chosen from the rigid body modes 1,...,6
!                   for each generalized mode to have same symmetry
!     VELH        Generalized normal for the radiation problem
!     ZDISP       Vertical displacement of the surface
! MODS:
!   1/07  optional parameter NPPW to specify panels per wavemaker
!   5/07  IGENMDS added, options for piston or tent-function modes
!   7/07  different filenames for .dat files for each mode type
!   5/08  gdf filename used in filename of .dat file
!-----------------------------------------------------------------------
      IMPLICIT NONE
      REAL XI(3),XN(3),VELH(*),ZDISP(*)
      INTEGER IBMOD(*),NEWMDS(*),IMODE(*),IBI,IPP,IFLAG,IER,IFILEDLL,   &
     &        IERROR,IWAMLOG,IMD,IGENMDS,NSEG  
      CHARACTER*16 FILENDLL           
      REAL, SAVE :: ZHINGE,Y1,DELY
      INTEGER, SAVE :: NPPW
      REAL YMIN,YMAX
      CHARACTER*80 HEADDAT
      CHARACTER*15 :: FILEN(3)=(/ '_wmkrhinge.dat',                     &
     &                  '_wmkrpiston.dat','_wmkrtent.dat' /)
      IF(IFLAG.EQ.-1) THEN
!-----------------------------------------------------------------------
!   Open input file, read data for wavemaker, copy file to wamitlog.txt
!     IF (IGENMDS==21) Read header line and ZHINGE 
!     IF (IGENMDS==211) Read header line (file is optional)
!     IF (IGENMDS==212) Read header line and Y1,DELY
!     Write complete input filename to wamitlog.txt file
!     Read optional parameter NPPW, if not present then NPPW=1
!-----------------------------------------------------------------------
        NPPW=1
        IF (IGENMDS==21) THEN	 
          OPEN (UNIT=IFILEDLL,FILE=TRIM(FILENDLL)//FILEN(1),            &
     &           STATUS='OLD',ERR=20)
          READ (IFILEDLL,900,ERR=89) HEADDAT
          READ (IFILEDLL,*,ERR=89) ZHINGE
          WRITE (IWAMLOG,910) TRIM(FILENDLL)//FILEN(1)
        ELSEIF (IGENMDS==211) THEN
          OPEN (UNIT=IFILEDLL,FILE=TRIM(FILENDLL)//FILEN(2),            &
     &           STATUS='OLD',ERR=20)
          READ (IFILEDLL,900,ERR=89) HEADDAT
          WRITE (IWAMLOG,910) TRIM(FILENDLL)//FILEN(2)
        ELSEIF (IGENMDS==212) THEN
          OPEN (UNIT=IFILEDLL,FILE=TRIM(FILENDLL)//FILEN(3),            &
     &           STATUS='OLD',ERR=20)
          READ (IFILEDLL,900,ERR=89) HEADDAT
          READ (IFILEDLL,*,ERR=89) Y1,DELY
          WRITE (IWAMLOG,910) TRIM(FILENDLL)//FILEN(3)
        ENDIF
        READ (IFILEDLL,*,END=5,ERR=5) NPPW
  5     REWIND (IFILEDLL)
!-----------------------------------------------------------------------
!   Copy inputs to log file
!-----------------------------------------------------------------------
        DO  
          READ (IFILEDLL,900,END=10,ERR=10) HEADDAT
          WRITE (IWAMLOG,900) TRIM(HEADDAT)
        END DO
  10    CLOSE (IFILEDLL)
        GOTO 30
  20    IF (IGENMDS.NE.211) GOTO 89
!-----------------------------------------------------------------------
!   Assign generalized-mode symmetry indices (same as heave)
!-----------------------------------------------------------------------
  30    IMODE(IBMOD(IBI)+7:IBMOD(IBI)+6+NEWMDS(IBI))=3
      ELSEIF ((IPP.LE.NPPW*NEWMDS(IBI)).AND.(IGENMDS<212)) THEN
!-----------------------------------------------------------------------
!  IGENMDS=21 or 211: Assign hinge or piston mode, normalized=1 at XI(3)=0
!     Patch index and generalized mode index are the same if NPPW=1
!     NSEG= patch index or index of wavemaker element
!-----------------------------------------------------------------------
        NSEG=(IPP-1)/NPPW+1
        IMD=IBMOD(IBI)+6+NSEG
        IF (IGENMDS==21) THEN	 
           VELH(IMD)=1.0-XI(3)/ZHINGE
        ELSE
           VELH(IMD)=1.0 
        ENDIF
      ELSEIF ((IPP.LE.NPPW*(1+NEWMDS(IBI))).AND.(IGENMDS==212)) THEN
!-----------------------------------------------------------------------
!  IGENMDS=212: Assign tent-function modes.  
!     NSEG= patch index or index of wavemaker element, which spans the
!     interval  YMIN<Y<YMAX, where YMIN=(Y1+(NSEG-1)*DELY), YMAX=Y1+NSEG*DELY
!     1st segment has only one mode, equal to (XI(2)-YMIN)/DELY
!     Last segment has only one mode, equal to (YMAX-XI(2))/DELY
!     All other segments have both modes
!-----------------------------------------------------------------------
        NSEG=(IPP-1)/NPPW+1
        YMAX=Y1+NSEG*DELY
        YMIN=YMAX-DELY
        IMD=IBMOD(IBI)+6+NSEG
        IF (NSEG.LE.NEWMDS(IBI)) VELH(IMD)=(XI(2)-YMIN)/DELY
        IF (NSEG>1) VELH(IMD-1)=(YMAX-XI(2))/DELY	  	          
      ENDIF
      GOTO 99
!-----------------------------------------------------------------------
!   Write error statement to monitor and error file, then stop
!     Patch index and generalized mode index are the same 
!-----------------------------------------------------------------------
  89  WRITE (*,900) ' Error in input file '//FILEN
      WRITE (IERROR,900) ' Error in input file '//FILEN
      IER=1
  99  RETURN
 900  FORMAT(A)
 910  FORMAT (/,' Input file ',A)
      END

      SUBROUTINE HINGE_MODES(XI,XN,IBMOD,IBI,IPP,NEWMDS,                &
     &     IMODE,VELH,ZDISP,IFLAG,IFILEDLL,FILENDLL,IERROR,IWAMLOG,IER)                                            &
!-----------------------------------------------------------------------
!
!     COPYRIGHT (C) 2006 WAMIT INC.
!
!-----------------------------------------------------------------------
!
!     Version : 6.4
!
!-----------------------------------------------------------------------
!
!     Source-code file : NEWMODES
!
!-----------------------------------------------------------------------
!
!     DESCRIPTION :  Vertical modes of a hinged structure, consisting of 
!                    NSEG rigid segments connected by NSEG-1 hinges
!                    Hinges are in the plane z=0, and parallel to the y-axis
!
!          Options:  No symmetry (ISX=0) or symmetric about x=0 (ISX=1), 
!                    ISX must correspond to ISX in the body GDF file
!                    The nodes (hinges and ends) are at x=XH(n) 
!                    Array XH(n) is input from special data file 
!     Format of special input file `XHINGE.DAT':
!     header (comment line)
!     ISX, NSEG
!     XH (array)  (x coordinates of ends of segments)
!
!   If ISX=0, XH(0:NSEG) has index running from 0 (left end) to NSEG
!        (right end)
!   If ISX=1, XH(1:NXH) has index running from 1 (first hinge at x>0) to 
!        NXH=(NSEG+1)/2  (right end)
!   (XH must be monotonic increasing, i.e. XH(N)>xH(N-1) for all elements)
!
!   Error return (IER=1) if XHINGE.DAT cannot be opened and read 
!      or if NSEG is too large relative to NEWMDS  (NSEG > NEWMDS+1)
!
!  N.B. This subroutine has access to three external files, with unit
!        numbers IFILEDLL, IWAMLOG, IERROR, which can be used as follows:
!    
!      IFILEDLL:  input file XHINGE.DAT  
!      IWAMLOG:   WAMITLOG.TXT (log file for WAMIT run)
!      IERROR:    ERRORP.LOG or ERRORF.LOG (error files for Poten/Force)
!-----------------------------------------------------------------------
!
!     Input Arguments (in alphabetical order):
!
!     Symbolic  Where value  Abbreviated Description (see routine where
!       name    is assigned    value is assigned for full description)
!     ------------------------------------------------------------------
!     IBMOD     GENMDS       Starting mode index for each body minus 1
!     IBI                    Body index
!     IPP                    Patch index
!     IFLAG                  IFLAG= -1 return IMODE only
!     NEWMDS    RDPOT        Number of generalized modes
!     XI        BODYGEOM     Cartesian coordinate in body coordinates
!     XN        BODYGEOM     Normal vector
!
!     Output Arguments (in alphabetical order):
!
!     Symbolic    Description
!       name
!     ------------------------------------------------------------------
!     IMODE       Mode indices chosen from the rigid body modes 1,...,6
!                   for each generalized mode to have same symmetry
!     VELH        Generalized normal for the radiation problem
!     ZDISP       Vertical displacement of the surface
!
! MODS:
!   5/08  gdf filename used in filename of .dat file
!-----------------------------------------------------------------------
      IMPLICIT NONE
      REAL XI(3),XN(3),VELH(*),ZDISP(*)
      INTEGER IBMOD(*),NEWMDS(*),IMODE(*),IBI,IPP,IFLAG,IER,IFILEDLL,   &
     &        IERROR,IWAMLOG         
      REAL, DIMENSION(:), ALLOCATABLE, SAVE :: XH
      INTEGER, SAVE :: NSEG, NXH, ISX, NSEGEO 
      INTEGER N,J1,J2
!-----------------------------------------------------------------------
!   Unit numbers for external files.  See note in header above.
!-----------------------------------------------------------------------
      CHARACTER*80 HEADDAT
      CHARACTER*16 FILENDLL
      CHARACTER*11 :: FILEN='_xhinge.dat'
      IF(IFLAG.EQ.-1) THEN
!-----------------------------------------------------------------------
!   Open input file with data for hinge coordinates 
!   Read data and copy file to wamitlog.txt
!   The indexing of the array XH depends on the symmetry index ISX and
!      number of segments NSEG (NSEG=total including reflection about x=0)
!   The two cases:
!      ISX=0: XH(0:NSEG) no plane of symmetry, XH includes ends and hinges
!      ISX=1: XH(1:NXH)  NXH=(NSEG+1)/2 is the number of hinges+end for x>0
!   When ISX=1 and NSEG is even,the hinge at x=0 is assumed, not input
!   NSEGEO=(0,1) for NSEG (even,odd)
!   Since XH is an array it must be allocated with the correct dimensions
!   Since this subroutine may be initialized again in FORCE, it is
!   necessary to use the IF (.NOT.ALLOCATED ...  statement
!-----------------------------------------------------------------------
        OPEN (UNIT=IFILEDLL,FILE=TRIM(FILENDLL)//FILEN,                 &
     &        STATUS='OLD',ERR=89)
        READ (IFILEDLL,900,ERR=89) HEADDAT
        READ (IFILEDLL,*,ERR=89) ISX,NSEG
        IF (NSEG>NEWMDS(IBI)+1) GOTO 88
        IF (ISX==0) THEN
           NXH=NSEG
           IF (.NOT.ALLOCATED(XH)) ALLOCATE (XH(0:NXH))
           READ (IFILEDLL,*,ERR=89) (XH(N),N=0,NXH)
        ELSE
           NXH=(NSEG+1)/2 
           NSEGEO=MOD(NSEG,2)
           IF (.NOT.ALLOCATED(XH)) ALLOCATE (XH(NXH))
           READ (IFILEDLL,*,ERR=89) (XH(N),N=1,NXH)
        ENDIF
        REWIND (IFILEDLL)
!-----------------------------------------------------------------------
!   Copy inputs to log file
!-----------------------------------------------------------------------
        WRITE (IWAMLOG,910) TRIM(FILENDLL)//FILEN
        DO  
          READ (IFILEDLL,900,END=10,ERR=10) HEADDAT
          WRITE (IWAMLOG,900) TRIM(HEADDAT)
        END DO
   10   CLOSE (IFILEDLL)
!-----------------------------------------------------------------------
!   Assign generalized-mode symmetry indices 
!     odd modes (J=7,9,11,... ) are symmetric like heave  
!     even modes (J=8,10,12,...) are anti-symmetric like surge)  
!-----------------------------------------------------------------------
        IMODE(7:6+NEWMDS(IBI):2)=3
        IMODE(8:6+NEWMDS(IBI):2)=1
      ELSE 
!-----------------------------------------------------------------------
!     Assign generalized modes 
!     First locate the segment, then assign the modes for it
!     J1/J2 are the first/last modes that are active on the segment
!     N=NXH after completing DO loop, if x is on last segment
!-----------------------------------------------------------------------
        DO N=1,NXH-1
          IF (XI(1).LE.XH(N)) EXIT
        ENDDO
!-----------------------------------------------------------------------
!   ISX=0  (no symmetry)
!-----------------------------------------------------------------------
        IF (ISX==0) THEN
          IF (N==1) THEN
            J1=7
            J2=7 	
            ZDISP(J2)=(XI(1)-XH(0))/(XH(N)-XH(N-1))
          ELSEIF (N<NSEG) THEN	   	  
	      J1=N+5
            J2=N+6	    	
            ZDISP(J1)=(XH(N)-XI(1))/(XH(N)-XH(N-1))
            ZDISP(J2)=(XI(1)-XH(N-1))/(XH(N)-XH(N-1))
          ELSE
	      J1=N+5
            J2=J1
            ZDISP(J1)=(XH(N)-XI(1))/(XH(N)-XH(N-1))  	            
          ENDIF      
!-----------------------------------------------------------------------
!  ISX=1, NSEG even
!-----------------------------------------------------------------------
        ELSEIF (NSEGEO==0) THEN
          IF (N==1) THEN
            J1=7
            J2=7 	
	      ZDISP(7)=(XH(1)-XI(1))/XH(1)
            IF (NSEG>2) THEN
              ZDISP(8:9)=XI(1)/XH(1)
              J2=9
            ENDIF
          ELSEIF (N<NXH) THEN	   	  
	      J1=2*N+4
            J2=J1+3	    	
	      ZDISP(J1:J1+1)=(XH(N)-XI(1))/(XH(N)-XH(N-1))
	      ZDISP(J2-1:J2)=(XI(1)-XH(N-1))/(XH(N)-XH(N-1))
          ELSE
	      J1=2*N+4
            J2=J1+1
	      ZDISP(J1:J2)=(XH(N)-XI(1))/(XH(N)-XH(N-1))  	            
          ENDIF    
!-----------------------------------------------------------------------
!  ISX=1, NSEG odd
!-----------------------------------------------------------------------
        ELSE 
          IF (N==1) THEN
            J1=7
            J2=8 
            ZDISP(7)=1.0	
            ZDISP(8)=XI(1)/XH(1)
          ELSEIF (N<NXH) THEN	   	  
	      J1=2*N+3
            J2=J1+3	    	
	      ZDISP(J1:J1+1)=(XH(N)-XI(1))/(XH(N)-XH(N-1))
	      ZDISP(J2-1:J2)=(XI(1)-XH(N-1))/(XH(N)-XH(N-1))
          ELSE
	      J1=2*N+3
            J2=J1+1
	      ZDISP(J1:J2)=(XH(N)-XI(1))/(XH(N)-XH(N-1))  	            
          ENDIF    
        ENDIF
        VELH(J1:J2)=XN(3)*ZDISP(J1:J2)
      ENDIF
      GOTO 99
!-----------------------------------------------------------------------
!   Write error statements to monitor and error file
!-----------------------------------------------------------------------
  88  WRITE (*,900) ' NSEG > NEWMDS+1 '
      WRITE (IERROR,900) ' NSEG > NEWMDS+1 '
      GOTO 98
  89  WRITE (*,900) ' Error in input file '//FILEN
      WRITE (IERROR,900) ' Error in input file '//FILEN
  98  IER=1
  99  RETURN
 900  FORMAT(A)
 910  FORMAT (/,' Input file ',A)
      END

      SUBROUTINE GAP_FS_FOURIER                                         &
     &      (XI,XN,IBMOD,IBI,IPP,NEWMDS,IGENMDS,IMODE,VELH,ZDISP,IFLAG)
!-----------------------------------------------------------------------
!
!     COPYRIGHT (C) 2004 WAMIT INC.
!
!-----------------------------------------------------------------------
!
!     Version : 6.2
!
!-----------------------------------------------------------------------
!
!     Source-code file : NEWMODES
!
!-----------------------------------------------------------------------
!
!     DESCRIPTION :  Adapted from GAP_FS  
!  Modes for lid on free surface between hulls of two independent hulls
!  Generalized modes represent vertical motions of the gap free surface:
!
!  Longitudinal Fourier modes X_n(x')=[1, sin(pi*x'), cos(pi*x'), .., 
!    sin(n*pi*x'), cos(n*pi*x').)
!    x'=[x-(xbow+xstern)/2]/[(xbow-xstearn)/2], xbow > xstearn
!    x' is normalzed value between -1 and 1
!    
!  Lateral modes: i) general or ii) standing wave in lateral direction
!
!   i) Y_m[1, sin(pi*y'), cos(pi*y'), .., 
!      sin(m*pi*y'), cos(m*pi*y').)
!
!   ii) (Y_m(y')=[1, sin(pi/2*y'), cos(pi*y'), sin(m_odd*pi/2*y') 
!                 cos(m_even*(pi/2)*y')
!   y'=y/HGAP (HGAP is half of the gap)
!   y' is normalized value between -1 and 1
!
!  
!  In the program
!    T=x'
!    XI(1)=x
!    X2=xbow
!    X1=xstearn
!
!    P=y'
!    XI(2)=y
!  
!  User must specify
!
!  X1: x coordinates of the bow 
!  X2: x coordinates of the stearn
!  HGAP: half gap (the body coordinate is at the mid point of th gap)
!  NLNG: Number of modes in longtudinal direction
!  MLTR: Number of modes in lateral direction
!  (NEWMDS(IBI)=NLOGNI*NLATERAL)
!  ISTANDING: =0 full Fourier series, =1 subset for standing waves 
!
!  Mode sequence is X_0*Y_0,X_1*Y0,...X_NLNG-1*Y0,,,,X_NLNG*Y_MLTR
! 
!-----------------------------------------------------------------------
!
!     Input Arguments (in alphabetical order):
!
!     Symbolic  Where value  Abbreviated Description (see routine where
!       name    is assigned    value is assigned for full description)
!     ------------------------------------------------------------------
!     IBMOD     GENMDS       Starting mode index for each body minus 1
!     IBI                    Body index
!     IPP                    Patch index
!     IFLAG                  IFLAG= -1 return IMODE only
!     NEWMDS    RDPOT        Number of generalized modes
!     XI        BODYGEOM     Cartesian coordinate in body coordinates
!     XN        BODYGEOM     Normal vector
!
!
!     Output Arguments (in alphabetical order):
!
!     Symbolic    Description
!       name
!     ------------------------------------------------------------------
!     IMODE       Mode indices chosen from the rigid body modes 1,...,6
!                   for each generalized mode to have same symmetry
!     VELH        Generalized normal for the radiation problem
!     ZDISP       Vertical displacement of the surface
!
!-----------------------------------------------------------------------
      IMPLICIT NONE
      REAL XI(3),XN(3),VELH(*),ZDISP(*)
      INTEGER IBMOD(*),NEWMDS(*),IMODE(*),IBI,IPP,IFLAG
      INTEGER J,IGENMDS
      INTEGER, SAVE :: J1,J2,J3,J4,J5,J6
      INTEGER NTHETA,NTHETAMOD,JL,JL2
      REAL T,FY,TNTHETA,P,PNTHETA
      REAL, PARAMETER :: PI2=1.570796,PI=3.141592654
      INTEGER, PARAMETER :: NLNG=16,MLTR=2,ISTANDING=1
      REAL, PARAMETER :: X1=-80,X2=80,HGAP=4.
      IF(IFLAG==-1) THEN
         J1=19
         J2=J1-1+NLNG
!-----------------------------------------------------------------------
!   Assign mode symmetry indices (arbitrary, since ISX=ISY=0)
!-----------------------------------------------------------------------
         IMODE(J1:J1-1+NEWMDS(IBI))=3
      ELSEIF (IBI == 3) THEN
!-----------------------------------------------------------------------
!     Assign longitudinal modes 
!     Modes (J1,J2) are for Y_0=constant
!-----------------------------------------------------------------------
         T=(2.*XI(1)-(X1+X2))/(X2-X1)
         ZDISP(J1)=1.0
         NTHETA=0
         DO J=J1+1,J2,2
           NTHETA=NTHETA+1
           TNTHETA=NTHETA*T
           NTHETAMOD=TNTHETA/2
           TNTHETA=TNTHETA-NTHETAMOD*2
           ZDISP(J)=SIN(TNTHETA*PI)
           IF(J+1.LE.J2) ZDISP(J+1)=COS(TNTHETA*PI)
         ENDDO           
         VELH(J1:J2)=ZDISP(J1:J2)
!-----------------------------------------------------------------------
!     Transverse modes
!-----------------------------------------------------------------------
         P=XI(2)/HGAP
         IF(ISTANDING==1) THEN
           NTHETA=0
           DO J=1,MLTR-1
             NTHETA=NTHETA+1
             PNTHETA=NTHETA*P
             NTHETAMOD=NTHETA/4
             PNTHETA=PNTHETA-NTHETAMOD*4
             IF(MOD(J,2)==1) THEN
               FY=SIN(PI2*PNTHETA)
             ELSE
               FY=COS(PI2*PNTHETA)
             ENDIF
             DO JL=J1,J2
               JL2=JL+J*NLNG
               ZDISP(JL2)=ZDISP(JL)*FY
               VELH(JL2)=ZDISP(JL2)
             ENDDO
           ENDDO
         ELSE
           NTHETA=0
           DO J=1,MLTR-1
             NTHETA=NTHETA+1
             PNTHETA=NTHETA*P
             NTHETAMOD=NTHETA/2
             PNTHETA=PNTHETA-NTHETAMOD*2
             IF(MOD(J,2)==1) THEN
               FY=SIN(PI*PNTHETA)
             ELSE
               FY=COS(PI*PNTHETA)
             ENDIF
             DO JL=J1,J2
               JL2=JL+J*NLNG
               ZDISP(JL2)=ZDISP(JL)*FY
               VELH(JL2)=ZDISP(JL2)
             ENDDO
           ENDDO
         ENDIF
      ENDIF
 99   RETURN
      END
