      SUBROUTINE GEOMXACT(U,V,IBI,IPI,IS,NPATCH,IRR,X,XU,XV,IGDEF,      &
     &                    IGDFSCRB,IER)
!DEC$ ATTRIBUTES DLLEXPORT :: GEOMXACT
!-----------------------------------------------------------------------
!
!     COPYRIGHT (C) 2000-2008  WAMIT INCORPORATED
!
!-----------------------------------------------------------------------
!
!     Version : 6.4
!
!-----------------------------------------------------------------------
!
!     Source-code file :  GEOMXACT.F
!
!-----------------------------------------------------------------------
!
!     DESCRIPTION :  
!
!     Define the body geometry exactly by analytic relations
!     This subroutine is intended to be extended by users for specific
!     applications where the geometry of the body can be defined
!     explicitly.  This option in WAMIT is specified by
!     setting the parameter IGDEF > 2 or IGDEF<0.  IGDEF<0 is used
!     for the subroutines provided by WAMIT Inc, and IGDEF>2 are reserved
!     for users.
!
!     This Version of GEOMXACT is distributed with the standard release
!     of WAMIT Version 6.4.  It includes 36 subroutines which are described
!     in the User Manual, Section 6.8.
!     See the code below and subroutine headers for further information.
!
!     Inputs:
!
!     IBI    body number
!     IGDEF  index specifying which subroutine/body is requested
!     IGDFSCRB unit number of file with GDF data, used only if IPI=0)
!     IPI    patch number (local value for body IBI)
!     IS(1:2) Symmetry indices of body IBI in GDF file
!     NPATCH Number of patches on body IBI in GDF file
!     U,V    parametric coordinates on the patch
!
!     Outputs:
!
!     X      Cartesian coordinates of the point on body surface
!     XU     Partial derivatives of X with respect to U
!     XV     Partial derivatives of X with respect to V
!
!     X,XU,XV are initialized to zero in the calling routine
!     Only nonzero values are required to be assigned here
!
!     MODS:
!     8/00     Dec directive to export GEOMXACT as DLL
!              Add IGDEF and IGDFSCRB to the argument lists
!              Dec directive temporarily blocked, old arguments restored
!              All subroutines revised or new
!     3/01     Error return added if IGDEF.NE.1 at end of IF block
!     9/05     Extended to be used with CSF 
!              Since CSF is called one, ICDEF(IBI) is passed via IGDEF
!              when GEOMXACT is called for CSF. 
!              Any of IGDEF(1:NBODY) should not be same as 
!              ICDEF in a run. If same geometry is used for GDF and CSF
!              make duplicate of the subroutine, one with different name
!     6/06     Calls from RDCSF,CSFGEOM changed to include IBI,ICDEF array
!              ICDEF and IBI<0 IF BLOCK removed
!     5/08     NPATCH, IS	 added to calling arguments
!     7/08     IRR(IBI) added to calling arguments.  IRR is not passed
!              for initial calls to initialize data (IPI=0).
!-----------------------------------------------------------------------
      IMPLICIT NONE
      INTEGER IGDEF(*),IER,IBI,IPI,IGDFSCRB,NPATCH,IS(2),IRR,JS(2)
      REAL U,V,X(3),XU(3),XV(3)
      IER=0
!-----------------------------------------------------------------------
!     Body geometries for standard release version (IGDEF= -1 to -32)
!-----------------------------------------------------------------------      
      IF (IGDEF(IBI).EQ.-1) THEN
         CALL CIRCCYL(U,V,IBI,IPI,X,XU,XV,IGDFSCRB,IER)
      ELSEIF (IGDEF(IBI).EQ.-2) THEN
         CALL ELLIPCYL(U,V,IBI,IPI,X,XU,XV,IGDFSCRB,IER)
      ELSEIF (IGDEF(IBI).EQ.-3) THEN
         CALL SPHERE(U,V,IBI,IPI,X,XU,XV,IGDFSCRB,IER)
      ELSEIF (IGDEF(IBI).EQ.-4) THEN
         CALL ELLIPSOID(U,V,IBI,IPI,X,XU,XV,IGDFSCRB,IER)
      ELSEIF (IGDEF(IBI).EQ.-5) THEN
         CALL BARGE(U,V,IBI,IPI,X,XU,XV,IGDFSCRB,IER)
      ELSEIF (IGDEF(IBI).EQ.-6) THEN
         CALL BARGEMP(U,V,IBI,IPI,X,XU,XV,IGDFSCRB,IER)
      ELSEIF (IGDEF(IBI).EQ.-7) THEN
         CALL CYLMP(U,V,IBI,IPI,X,XU,XV,IGDFSCRB,IER)
      ELSEIF (IGDEF(IBI).EQ.-8) THEN
         CALL TORUS(U,V,IBI,IPI,X,XU,XV,IGDFSCRB,IER)
      ELSEIF (IGDEF(IBI).EQ.-9) THEN
         CALL TLP(U,V,IBI,IPI,X,XU,XV,IGDFSCRB,IER)
      ELSEIF (IGDEF(IBI).EQ.-10) THEN
         CALL SEMISUB(U,V,IBI,IPI,X,XU,XV,IGDFSCRB,IER)
      ELSEIF (IGDEF(IBI).EQ.-11) THEN
         CALL FPSO(U,V,IBI,IPI,X,XU,XV,IGDFSCRB,IER)
      ELSEIF (IGDEF(IBI).EQ.-12) THEN
         CALL SPAR(U,V,IBI,IPI,X,XU,XV,IGDFSCRB,IER)
      ELSEIF (IGDEF(IBI).EQ.-13) THEN
         CALL AUV(U,V,IBI,IPI,X,XU,XV,IGDFSCRB,IER)
      ELSEIF (IGDEF(IBI).EQ.-14) THEN
         CALL SPAR2(U,V,IBI,IPI,X,XU,XV,IGDFSCRB,IER)
      ELSEIF (IGDEF(IBI).EQ.-15) THEN
         CALL SPHERXYZ(U,V,IBI,IPI,IS,X,XU,XV,IGDFSCRB,IER)
      ELSEIF (IGDEF(IBI).EQ.-16) THEN
         CALL FPSO2(U,V,IBI,IPI,X,XU,XV,IGDFSCRB,IER)
      ELSEIF (IGDEF(IBI).EQ. -17) THEN
         CALL FPSO12(U,V,IBI,IPI,X,XU,XV,IGDFSCRB,IER)
      ELSEIF (IGDEF(IBI).EQ.-18) THEN
         CALL TORUS_ELLIP(U,V,IBI,IPI,X,XU,XV,IGDFSCRB,IER)
      ELSEIF (IGDEF(IBI).EQ. -19) THEN
         CALL TORUS2(U,V,IBI,IPI,X,XU,XV,IGDFSCRB,IER)
      ELSEIF (IGDEF(IBI).EQ.-20) THEN
         CALL CIRCCYLH(U,V,IBI,IPI,X,XU,XV,IGDFSCRB,IER)
      ELSEIF (IGDEF(IBI).EQ.-21) THEN
         CALL FPSOINT(U,V,IBI,IPI,IRR,X,XU,XV,IGDFSCRB,IER)
      ELSEIF (IGDEF(IBI).EQ.-22) THEN
         CALL CIRCCYL_ARRAY(U,V,IBI,IPI,X,XU,XV,IGDFSCRB,IER)
      ELSEIF (IGDEF(IBI).EQ.-23) THEN
         CALL ELLIPINT(U,V,IBI,IPI,IRR,X,XU,XV,IGDFSCRB,IER)
      ELSEIF (IGDEF(IBI).EQ.-24) THEN
         CALL GAPLID(U,V,IBI,IPI,X,XU,XV,IGDFSCRB,IER)
      ELSEIF (IGDEF(IBI).EQ.-25) THEN
         CALL CYLFIN(U,V,IBI,IPI,X,XU,XV,IGDFSCRB,IER)
      ELSEIF (IGDEF(IBI).EQ.-26) THEN
         CALL CYLFIN4(U,V,IBI,IPI,IS,NPATCH,X,XU,XV,IGDFSCRB,IER)
      ELSEIF (IGDEF(IBI).EQ.-27) THEN
         CALL SKEW_SPHERE(U,V,IBI,IPI,X,XU,XV,IGDFSCRB,IER)
      ELSEIF (IGDEF(IBI).EQ.-28) THEN
         CALL CIRCCYL_NOSYM(U,V,IBI,IPI,X,XU,XV,IGDFSCRB,IER)
      ELSEIF (IGDEF(IBI).EQ.-29) THEN
         CALL ELLIPSOID_NOSYM_TANK(U,V,IBI,IPI,IRR,X,XU,XV,IGDFSCRB,IER)
      ELSEIF (IGDEF(IBI).EQ.-30) THEN
         CALL BARGE_INT(U,V,IBI,IPI,X,XU,XV,IGDFSCRB,IER)
      ELSEIF (IGDEF(IBI).EQ.-31) THEN
         CALL BARGENUC(U,V,IBI,IPI,X,XU,XV,IGDFSCRB,IER)
      ELSEIF (IGDEF(IBI).EQ.-32) THEN
	  CALL CCYLHSP(U,V,IBI,IPI,X,XU,XV,IGDFSCRB,IER)
!-----------------------------------------------------------------------
!  Control surface geometries for standard release version:
!    Circular cylinder with annular free surface: 2 planes of symmetry
!    Circular cylinder with annular free surface: no planes of symmetry
!    Ellipsoid with annular free surface: 2 planes of symmetry
!    Ellipsoid with annular free surface: no planes of symmetry
!-----------------------------------------------------------------------
      ELSEIF(IGDEF(IBI)==-1001) THEN
         CALL CCYL_CS(U,V,IBI,IPI,X,XU,XV,IGDFSCRB,IER) ! circular cylinder 
      ELSEIF(IGDEF(IBI)==-1002) THEN
         CALL CCYL_CS_NOSYSM(U,V,IBI,IPI,X,XU,XV,IGDFSCRB,IER) ! cir cylinder 
      ELSEIF(IGDEF(IBI)==-1003) THEN
         CALL ELLIPSOID_CS(U,V,IBI,IPI,X,XU,XV,IGDFSCRB,IER) ! ellipsoid 
      ELSEIF(IGDEF(IBI)==-1004) THEN
         CALL ELLIPSOID_CS_NOSYM(U,V,IBI,IPI,X,XU,XV,IGDFSCRB,IER) ! ellipsoid 
      ELSEIF (IGDEF(IBI).NE.1) THEN
         IER=1
      ENDIF
      IF (IER.NE.0) IER=1
      RETURN
      END

      SUBROUTINE CIRCCYL(U,V,IBI,IPI,X,XU,XV,IGDFSCRB,IER)
!-----------------------------------------------------------------------
!
!     COPYRIGHT (C) 2000      WAMIT INCORPORATED
!
!-----------------------------------------------------------------------
!
!     Version : 6.0
!
!-----------------------------------------------------------------------
!
!     Source-code file :  GEOMXACT.F
!
!-----------------------------------------------------------------------
!
!     DESCRIPTION :  (IGDEF=-1)
!
!     Subroutine defines the first quadrant of a circular cylinder
!     RADIUS and DRAFT are read from GDF input file initially
!     INONUMAP             read from GDF input file initially
!
!     Patch 1: side
!     Patch 2: bottom
!     Patch 3: interior free surface (optional)
!
!     Options:
!       INONUMAP=0 use uniform mapping on patches
!       INONUMAP=1 use nonuniform mapping on patches (see comments below)
!
!       Use NPATCH=1 for bottom-mounted cylinder (draft=depth)
!                    or for disc on free surface (draft=0.0)
!       Use NPATCH=2 for conventional case
!       Use NPATCH=3 for irregular-frequency removal (IRR=1)
!
!       6/02  Sign corrected for NPATCH=3 and INONUMAP=1
!-----------------------------------------------------------------------
      IMPLICIT NONE
      REAL X(3),XU(3),XV(3),U,V
      INTEGER IPI,IBI,IGDFSCRB,IER
!-----------------------------------------------------------------------
!  user-assigned definitions -- inputs from GDF must have SAVE attribute
!  HALFDRAFT, HALFRAD are evaluated/saved at first call
!-----------------------------------------------------------------------
      REAL, SAVE :: RADIUS,DRAFT,HALFDRAFT,HALFRAD,QRAD,QDRAFT,THREEQDR
      INTEGER, SAVE :: INONUMAP
      REAL THETA,CTHETA,STHETA,SR,ONEPMV,XV3FAC
!-----------------------------------------------------------------------
!   PIO4=pi/4
!-----------------------------------------------------------------------
      REAL, PARAMETER :: PIO4=0.785398163397E0
!-----------------------------------------------------------------------
!   Tolerance for zero draft input
!-----------------------------------------------------------------------
      REAL, PARAMETER :: TOLD=1.0E-8
!-----------------------------------------------------------------------
!   On initial call read data from GDF and evaluate parameters to SAVE
!   On all subsequent calls evaluate azimuthal angle and cos/sine
!-----------------------------------------------------------------------
      IF(IPI.EQ.0) THEN
        READ (IGDFSCRB,*,IOSTAT=IER) RADIUS,DRAFT
        READ (IGDFSCRB,*,IOSTAT=IER) INONUMAP
        IF (IER.NE.0) GOTO 99
        HALFRAD=0.5E0*RADIUS
        QRAD=0.25E0*RADIUS
        HALFDRAFT=0.5E0*DRAFT
        QDRAFT=0.25E0*DRAFT
        THREEQDR=0.75E0*DRAFT
!-----------------------------------------------------------------------
!   For all patches, first evaluate the angle THETA in quadrant one:
!      U=-1 <--> THETA=0,   U=+1 <--> THETA=PI/2
!-----------------------------------------------------------------------
      ELSE
         THETA=PIO4*(U+1.E0)
         CTHETA=COS(THETA)
         STHETA=SIN(THETA)
      ENDIF
!-----------------------------------------------------------------------
!  Patch 1: side        V=-1 at free surface, V=+1 at corner
!    If INONUMAP=0 X(3) is a linear function of V
!    If INONUMAP=1 X(3) is a cubic in V with derivatives=0 at V=+1,-1
!-----------------------------------------------------------------------
      IF((IPI.EQ.1).AND.(DRAFT.GE.TOLD)) THEN
         X(1)=RADIUS*CTHETA
         X(2)=RADIUS*STHETA
         XU(1)=-PIO4*X(2)
         XU(2)= PIO4*X(1)
         IF (INONUMAP==0) THEN
            X(3)=-HALFDRAFT*(V+1.E0)
            XV(3)=-HALFDRAFT
         ELSE
            ONEPMV=V+1.E0
            X(3)= QDRAFT*(V-2.E0)*ONEPMV*ONEPMV
            XV(3)=THREEQDR*ONEPMV*(V-1.E0)
         ENDIF
!-----------------------------------------------------------------------
!  Patch 2: bottom      V=+1 on axis, -1 at corner
!    If INONUMAP=0 SR (radius) is a linear function of V
!    If INONUMAP=1 SR is quadratic in V with derivative=0 at V=-1
!-----------------------------------------------------------------------
      ELSEIF((IPI.EQ.2).OR.(DRAFT<TOLD)) THEN
         IF (INONUMAP==0) THEN
            SR=HALFRAD*(1.E0-V)
            XV3FAC=-HALFRAD
         ELSE
            ONEPMV=1.E0+V
            SR=QRAD*(4.E0-ONEPMV*ONEPMV)
            XV3FAC=-HALFRAD*ONEPMV
         ENDIF
         X(1)=SR*CTHETA
         X(2)=SR*STHETA
         X(3)=-DRAFT
         XU(1)=-PIO4*X(2)
         XU(2)= PIO4*X(1)
         XV(1)=XV3FAC*CTHETA
         XV(2)=XV3FAC*STHETA
!-----------------------------------------------------------------------
!  Patch 3: interior free surface   V=-1 on axis,  V=+1 on side
!-----------------------------------------------------------------------
      ELSEIF(IPI.EQ.3) THEN
         IF (INONUMAP==0) THEN
            SR=HALFRAD*(1.E0+V)
            XV3FAC=HALFRAD
         ELSE
            ONEPMV=1.E0-V
            SR=QRAD*(4.E0-ONEPMV*ONEPMV)
            XV3FAC=HALFRAD*ONEPMV
         ENDIF
         X(1)=SR*CTHETA
         X(2)=SR*STHETA
         XU(1)=-PIO4*X(2)
         XU(2)= PIO4*X(1)
         XV(1)=XV3FAC*CTHETA
         XV(2)=XV3FAC*STHETA
      ENDIF
 99   RETURN
      END

      SUBROUTINE ELLIPCYL(U,V,IBI,IPI,X,XU,XV,IGDFSCRB,IER)
!-----------------------------------------------------------------------
!
!     COPYRIGHT (C) 2000      WAMIT INCORPORATED
!
!-----------------------------------------------------------------------
!
!     Version : 6.0
!
!-----------------------------------------------------------------------
!
!     Source-code file :  GEOMXACT.F
!
!-----------------------------------------------------------------------
!
!     DESCRIPTION :  (IGDEF=-2)
!
!     Subroutine defines the first quadrant of an elliptical cylinder
!     with semi-axes A,B in the X,Y directions.
!     Dimensions A,B,DRAFT are read from GDF input file
!
!     Extended from subroutine CIRCCYL
!
!     Patch 1: side
!     Patch 2: bottom
!     Patch 3: interior free surface (optional)
!
!     Options:
!        Use NPATCH=1 for bottom-mounted cylinder (draft=depth)
!        Use NPATCH=2 for conventional case
!        Use NPATCH=3 for irregular-frequency removal (IRR=1)
!
!-----------------------------------------------------------------------
      IMPLICIT NONE
      REAL X(3),XU(3),XV(3),U,V
      INTEGER IPI,IBI,IGDFSCRB,IER
      REAL, SAVE :: A, B, DRAFT, HALFA, HALFB, HALFDRAFT
      REAL THETA,CTHETA,STHETA,SR,SRA,SRB
      REAL, PARAMETER :: PIO4=0.785398163397E0
      IF(IPI.EQ.0) THEN
         READ (IGDFSCRB,*,IOSTAT=IER) A,B,DRAFT
         IF (IER.NE.0) GOTO 99
         HALFDRAFT=0.5E0*DRAFT
         HALFA=0.5E0*A
         HALFB=0.5E0*B
      ELSE
         THETA=PIO4*(U+1.E0)
         CTHETA=COS(THETA)
         STHETA=SIN(THETA)
      ENDIF
!-----------------------------------------------------------------------
!  Patch 1: side
!-----------------------------------------------------------------------
      IF(IPI.EQ.1) THEN
         X(1)=A*CTHETA
         X(2)=B*STHETA
         X(3)=-HALFDRAFT*(V+1.E0)
         XU(1)=-PIO4*A*STHETA
         XU(2)= PIO4*B*CTHETA
         XV(3)=-HALFDRAFT
!-----------------------------------------------------------------------
!  Patch 2: bottom
!-----------------------------------------------------------------------
      ELSEIF(IPI.EQ.2) THEN
         SR=1.E0-V
         SRA=SR*HALFA
         SRB=SR*HALFB
         X(1)=SRA*CTHETA
         X(2)=SRB*STHETA
         X(3)=-DRAFT
         XU(1)=-PIO4*SRA*STHETA
         XU(2)= PIO4*SRB*CTHETA
         XV(1)=-HALFA*CTHETA
         XV(2)=-HALFB*STHETA
!-----------------------------------------------------------------------
!  Patch 3: interior free surface
!-----------------------------------------------------------------------
      ELSEIF(IPI.EQ.3) THEN
         SR=V+1.E0
         SRA=SR*HALFA
         SRB=SR*HALFB
         X(1)=SRA*CTHETA
         X(2)=SRB*STHETA
         XU(1)=-PIO4*SRA*STHETA
         XU(2)= PIO4*SRB*CTHETA
         XV(1)=HALFA*CTHETA
         XV(2)=HALFB*STHETA
      ENDIF
 99   RETURN
      END

      SUBROUTINE SPHERE(U,V,IBI,IPI,X,XU,XV,IGDFSCRB,IER)
!-----------------------------------------------------------------------
!
!     COPYRIGHT (C) 2000-2008      WAMIT INCORPORATED
!
!-----------------------------------------------------------------------
!
!     Version : 6.4
!
!-----------------------------------------------------------------------
!
!     Source-code file :  GEOMXACT.F
!
!-----------------------------------------------------------------------
!
!     DESCRIPTION :  (IGDEF=-3)
!
!     Subroutine defines one quadrant of a floating hemisphere
!     Center in free surface
!     Radius input from GDF file
!
!     Patch 1: body surface
!     Patch 2: interior free surface (optional)
!
!     Options:
!        Use NPATCH=1 for conventional case
!        Use NPATCH=2 for irregular-frequency removal (IRR=1)
!        INONUMAP=0 (default) uniform mapping in azimuthal direction
!        INONUMAP=1 nonuniform mapping near waterline, both patches
!
!   MODS
!     10/08 INUONUMAP option added
!-----------------------------------------------------------------------
      IMPLICIT NONE
      REAL X(3),XU(3),XV(3),U,V
      INTEGER IPI,IBI,IGDFSCRB,IER
      INTEGER, SAVE :: INONUMAP      
      REAL, SAVE :: RADIUS,HALFRAD
      REAL THETA,CTHETA,STHETA,PSI,RADCPSI,RADSPSI,PRADCPSI,SR,VP1,ONEMV
!-----------------------------------------------------------------------
!   PI2=pi/2, PI4=pi/4, PI8=pi/8
!-----------------------------------------------------------------------
      REAL, PARAMETER :: PI2=1.5707963E0, PI4=0.7853981634E0,           &
     &                   PI8=0.39269908E0
      IF(IPI.EQ.0) THEN
         READ (IGDFSCRB,*,IOSTAT=IER) RADIUS
         IF (IER.NE.0) GOTO 99
         READ (IGDFSCRB,*,END=2) INONUMAP
         GOTO 4
  2      INONUMAP=0
  4      HALFRAD=0.5E0*RADIUS
      ELSE
         THETA=PI4*(U+1.E0)
         CTHETA=COS(THETA)
         STHETA=SIN(THETA)
      ENDIF
!-----------------------------------------------------------------------
!  Patch 1: body surface  V=-1: waterline (psi=pi/2)  V=+1: bottom axis (psi=pi)
!-----------------------------------------------------------------------
      IF(IPI.EQ.1) THEN
         IF (INONUMAP==0) THEN
            PSI=PI4*(V+3.E0)
         ELSE 
            VP1=V+1.E0
            PSI=PI2+PI8*VP1*VP1         
         ENDIF   
         RADCPSI=RADIUS*COS(PSI)
         RADSPSI=RADIUS*SIN(PSI)
         X(1)=RADSPSI*CTHETA
         X(2)=RADSPSI*STHETA
         X(3)=RADCPSI
         XU(1)=-PI4*X(2)
         XU(2)= PI4*X(1)
         IF (INONUMAP==0) THEN         
            PRADCPSI=PI4*RADCPSI
            XV(1)=PRADCPSI*CTHETA
            XV(2)=PRADCPSI*STHETA
            XV(3)=-PI4*RADSPSI
         ELSE 
            PRADCPSI=PI4*RADCPSI*VP1
            XV(1)=PRADCPSI*CTHETA
            XV(2)=PRADCPSI*STHETA
            XV(3)=-PI4*RADSPSI*VP1        
         ENDIF            
!-----------------------------------------------------------------------
!  Patch 2: interior free surface: V=-1: axis,  V=+1: waterline
!-----------------------------------------------------------------------
      ELSEIF (IPI.EQ.2) THEN
         IF (INONUMAP==0) THEN 
            SR=HALFRAD*(V+1.)
         ELSE
            ONEMV=1.-V
            SR=RADIUS*(1.-.25*ONEMV*ONEMV)
         ENDIF      
         X(1)=SR*CTHETA
         X(2)=SR*STHETA
         XU(1)=-PI4*X(2)
         XU(2)= PI4*X(1)
         IF (INONUMAP==0) THEN
            XV(1)=HALFRAD*CTHETA
            XV(2)=HALFRAD*STHETA
         ELSE
            XV(1)=HALFRAD*CTHETA*ONEMV
            XV(2)=HALFRAD*STHETA*ONEMV
         ENDIF   
      ENDIF
 99   RETURN
      END

      SUBROUTINE ELLIPSOID(U,V,IBI,IPI,X,XU,XV,IGDFSCRB,IER)
!-----------------------------------------------------------------------
!
!     COPYRIGHT (C) 2000      WAMIT INCORPORATED
!
!-----------------------------------------------------------------------
!
!     Version : 6.0
!
!-----------------------------------------------------------------------
!
!     Source-code file :  GEOMXACT.F
!
!-----------------------------------------------------------------------
!
!     DESCRIPTION :  (IGDEF=-4)
!
!
!     Subroutine defines the ellipsoid (IPI=1) by one patch
!     Center in free surface
!     Semi-axes A,B,C input from GDF file
!
!     Patch 1: body surface
!     Patch 2: interior free surface (optional)
!
!     Options:
!        Use NPATCH=1 for conventional case
!        Use NPATCH=2 for irregular-frequency removal (IRR=1)
!
!-----------------------------------------------------------------------
      IMPLICIT NONE
      REAL X(3),XU(3),XV(3),U,V
      INTEGER IPI,IBI,IGDFSCRB,IER
!-----------------------------------------------------------------------
!  user-assigned definitions -- inputs from GDF must have SAVE attribute
!  PI/4*(A,B,C) are evaluated once to reduce CPU time
!-----------------------------------------------------------------------
      REAL, SAVE :: A,B,C,PIO4A,PIO4B,PIO4C,HALFA,HALFB
      REAL THETA,CTHETA,STHETA,PSI,CPSI,SPSI,SR,SRA,SRB
      REAL, PARAMETER :: PIO4=0.785398163397E0
      IF (IPI == 0) THEN
        READ (IGDFSCRB,*,IOSTAT=IER) A,B,C
        IF (IER.NE.0) GOTO 99
        PIO4A=PIO4*A
        PIO4B=PIO4*B
        PIO4C=PIO4*C
        HALFA=0.5E0*A
        HALFB=0.5E0*B
      ELSE
        THETA=PIO4*(U+1.E0)
        CTHETA=COS(THETA)
        STHETA=SIN(THETA)
      ENDIF
!-----------------------------------------------------------------------
!  Patch 1: body surface
!-----------------------------------------------------------------------
      IF (IPI == 1) THEN
         PSI=PIO4*(V+3.E0)
         CPSI=COS(PSI)
         SPSI=SIN(PSI)
         X(1)=A*SPSI*CTHETA
         X(2)=B*SPSI*STHETA
         X(3)=C*CPSI
         XU(1)=-PIO4A*SPSI*STHETA
         XU(2)= PIO4B*SPSI*CTHETA
         XV(1)=PIO4A*CPSI*CTHETA
         XV(2)=PIO4B*CPSI*STHETA
         XV(3)=-PIO4C*SPSI
!-----------------------------------------------------------------------
!  Patch 2: interior free surface
!-----------------------------------------------------------------------
      ELSEIF (IPI == 2) THEN
         SR=V+1.E0
         SRA=SR*HALFA
         SRB=SR*HALFB
         X(1)=SRA*CTHETA
         X(2)=SRB*STHETA
         XU(1)=-PIO4*SRA*STHETA
         XU(2)= PIO4*SRB*CTHETA
         XV(1)=HALFA*CTHETA
         XV(2)=HALFB*STHETA
      ENDIF
 99   RETURN
      END

      SUBROUTINE BARGE(U,V,IBI,IPI,X,XU,XV,IGDFSCRB,IER)
!-----------------------------------------------------------------------
!
!     COPYRIGHT (C) 2000-03      WAMIT INCORPORATED
!
!-----------------------------------------------------------------------
!
!     Version : 6.2
!
!-----------------------------------------------------------------------
!
!     Source-code file :  GEOMXACT.F
!
!-----------------------------------------------------------------------
!
!     DESCRIPTION :  (IGDEF=-5)
!
!     Subroutine defines the first quadrant of a rectangular barge
!
!     Inputs from GDF file:
!        HALFLEN = .5*Length
!        HALFBEAM= .5*Beam
!        DRAFT   = Draft
!
!    NORMAL CASE (DRAFT>=TOL):
!      Patch 1: end
!      Patch 2: side
!      Patch 3: bottom
!      Patch 4: internal free surface (optional, use if IRR=1)
!
!     Options:
!        Use NPATCH=2 for bottom-mounted caisson (draft=depth)
!        Use NPATCH=3 for conventional case
!        Use NPATCH=4 for irregular-frequency removal (IRR=1)
!
!    SPECIAL CASE (DRAFT<TOL):
!      Patch 1: bottom
!
!       INONUMAP=0 use uniform mapping on patches (default)
!       INONUMAP=1,2 use nonuniform mapping on patches (see comments below)
!                  only installed for transverse mapping on side/bottom
!       4/03  option added for INONUMAP=1,2
!       6/03  option added for DRAFT=0.0
!-----------------------------------------------------------------------
      IMPLICIT NONE
      REAL X(3),XU(3),XV(3),U,V
      INTEGER IPI,IBI,IGDFSCRB,IER
      REAL, SAVE :: HALFLEN,HALFBEAM,DRAFT,QDRAFT,THREEQDR,             &
     &              HALFDRAFT,QUARTLEN,QUARTBEAM,EIGHTHBEAM
      INTEGER, SAVE :: INONUMAP
      REAL ONEPMV,XV3FAC
      REAL, PARAMETER :: TOL=1.E-6
      IF(IPI.EQ.0) THEN
         READ (IGDFSCRB,*,IOSTAT=IER) HALFLEN,HALFBEAM,DRAFT
         IF (IER.NE.0) GOTO 99
         READ (IGDFSCRB,*,END=2) INONUMAP
         GOTO 4
  2      INONUMAP=0
  4      HALFDRAFT=0.5E0*DRAFT
         QDRAFT=0.25E0*DRAFT
         THREEQDR=0.75E0*DRAFT
         QUARTLEN=0.5E0*HALFLEN
         QUARTBEAM=0.5E0*HALFBEAM
         EIGHTHBEAM=0.25E0*HALFBEAM
!-----------------------------------------------------------------------
!   IPI=1 or 3: BOTTOM
!    If INONUMAP=0 X(2) is a linear function of V
!    If INONUMAP=1 X(2) is quadratic in V with derivative=0 at V=+1
!-----------------------------------------------------------------------
      ELSEIF((DRAFT<TOL).OR.(IPI.EQ.3)) THEN
         X(1)=(U+1.E0)*QUARTLEN
         X(3)=-DRAFT
         XU(1)=QUARTLEN
         IF (INONUMAP==0) THEN
            X(2)=(V+1.E0)*QUARTBEAM
            XV(2)=QUARTBEAM
         ELSE
            ONEPMV=1.E0-V
            X(2)=EIGHTHBEAM*(4.E0-ONEPMV*ONEPMV)
            XV(2)=QUARTBEAM*ONEPMV
         ENDIF
!-----------------------------------------------------------------------
!   IPI=1: END
!-----------------------------------------------------------------------
      ELSEIF(IPI.EQ.1) THEN
         X(1)=HALFLEN
         X(2)=(U+1.E0)*QUARTBEAM
         X(3)=-HALFDRAFT*(V+1.E0)
         XU(2)=QUARTBEAM
         XV(3)=-HALFDRAFT
!-----------------------------------------------------------------------
!   IPI=2: SIDE
!    If INONUMAP=0 X(3) is a linear function of V
!    If INONUMAP=1 X(3) is quadratic in V with derivative=0 at V=+1
!    If INONUMAP=2 X(3) is a cubic in V with derivatives=0 at V=+1,-1
!-----------------------------------------------------------------------
      ELSEIF(IPI.EQ.2) THEN
         X(1)=(1.E0-U)*QUARTLEN
         X(2)=HALFBEAM
         XU(1)=-QUARTLEN
         IF (INONUMAP==0) THEN
            X(3)=-HALFDRAFT*(V+1.E0)
            XV(3)=-HALFDRAFT
         ELSEIF (INONUMAP==1) THEN
            ONEPMV=1.E0-V
            X(3)=-QDRAFT*(4.E0-ONEPMV*ONEPMV)
            XV(3)=-HALFDRAFT*ONEPMV
         ELSE
            ONEPMV=V+1.E0
            X(3)= QDRAFT*(V-2.E0)*ONEPMV*ONEPMV
            XV(3)=THREEQDR*ONEPMV*(V-1.E0)
         ENDIF
!-----------------------------------------------------------------------
!   IPI=4: INTERIOR FREE SURFACE
!-----------------------------------------------------------------------
      ELSEIF(IPI.EQ.4) THEN
         X(1)=(U+1.E0)*QUARTLEN
         X(2)=(1.E0-V)*QUARTBEAM
         XU(1)=QUARTLEN
         XV(2)=-QUARTBEAM
      ENDIF
 99   RETURN
      END

      SUBROUTINE BARGEMP(U,V,IBI,IPI,X,XU,XV,IGDFSCRB,IER)
!-----------------------------------------------------------------------
!
!     COPYRIGHT (C) 2000 WAMIT Inc
!
!-----------------------------------------------------------------------
!
!     Source-code file :  GEOMXACT.F
!
!-----------------------------------------------------------------------
!
!     VERSION : 6.0
!
!-----------------------------------------------------------------------
!
!     DESCRIPTION :  (IGDEF=-6)
!
!     Subroutine defines the geometry of one quadrant of
!     a rectangular barge with a rectangular moonpool at its center.
!     All relevant dimensions are read from the GDF file during the
!     initial call with the patch index IPI=0.  The variables storing
!     these dimensions are defined with the attribute SAVE so that they
!     are available for subsequent calls to the subroutine.
!
!     The following parameters are input from GDF file
!     to prescribe the geometry:
!
!     HALFLEN : Half-Length (x-coordinate of bow)
!     HALFBEAM: Half-Beam   (y-coordinate of side)
!     DRAFT: Draft
!     XMP: x-coordinate of the forward face of moon pool
!     YMP: transverse coordinate of the moon pool side (half width)
!
!     The following patches are defined:
!        IPI=1 bow end X=HALFLEN=L/2
!        IPI=2 side Y=HALFBEAM=B/2
!        IPI=3 bottom surface outboard of moon pool (YMP<Y<B/2)
!        IPI=4 outbd surface of moon pool 1
!        IPI=5 fwd surface of moon pool 1
!        IPI=6 bottom surface fwd of moon pool 1
!        IPI=7 free surface of moon pool (optional to damp resonance)
!
!   If NPATCH=6 the moonpool free surface is a physical free surface
!   If NPATCH=7 generalized modes must be used to describe the moonpool
!               (this option permits the input of external damping
!                to suppress moonpool reonances)
!
!-----------------------------------------------------------------------
      IMPLICIT NONE
      REAL X(3),XU(3),XV(3),U,V
      INTEGER IPI,IBI,IGDFSCRB,IER
      REAL, SAVE :: HALFLEN,HALFBEAM,DRAFT,XMP,YMP,                     &
     &              HALFDRAFT,QUARTLEN,QUARTBEAM,HALFXMP,HALFYMP
!-----------------------------------------------------------------------
!  Read data from GDF file and evaluate SAVE parameters
!-----------------------------------------------------------------------
      IF(IPI.EQ.0) THEN
        READ (IGDFSCRB,*,IOSTAT=IER) HALFLEN,HALFBEAM,DRAFT
        READ (IGDFSCRB,*,IOSTAT=IER) XMP,YMP
        IF (IER.NE.0) GOTO 99
        HALFDRAFT=0.5E0*DRAFT
        QUARTLEN=0.5E0*HALFLEN
        QUARTBEAM=0.5E0*HALFBEAM
        HALFXMP=0.5E0*XMP
        HALFYMP=0.5E0*YMP
      ELSEIF(IPI.EQ.1) THEN
!-----------------------------------------------------------------------
!                                        IPI=1 bow end
!-----------------------------------------------------------------------
         XU(2)=QUARTBEAM
         XV(3)=-HALFDRAFT
         X(1)=HALFLEN
         X(2)=(U+1.E0)*XU(2)
         X(3)=(V+1.E0)*XV(3)
      ELSEIF(IPI.EQ.2) THEN
!-----------------------------------------------------------------------
!                                        IPI=2 side Y=B/2
!-----------------------------------------------------------------------
         XU(1)=-QUARTLEN
         XV(3)=-HALFDRAFT
         X(1)=(U-1.E0)*XU(1)
         X(2)=HALFBEAM
         X(3)=(V+1.E0)*XV(3)
      ELSEIF(IPI.EQ.3) THEN
!-----------------------------------------------------------------------
!                IPI=3 bottom surface outboard of moon pool
!-----------------------------------------------------------------------
         XU(1)=QUARTLEN
         XV(2)=QUARTBEAM-HALFYMP
         X(1)=(U+1.E0)*XU(1)
         X(2)=(V+1.E0)*XV(2)+YMP
         X(3)=-DRAFT
      ELSEIF(IPI.EQ.4) THEN
!-----------------------------------------------------------------------
!                IPI=4 side of moon pool
!-----------------------------------------------------------------------
         XV(3)=-HALFDRAFT
         X(3)=(V+1.E0)*XV(3)
         XU(1)=HALFXMP
         X(1)=(U+1.E0)*XU(1)
         X(2)=YMP
      ELSEIF(IPI.EQ.5) THEN
!-----------------------------------------------------------------------
!                fwd surface of moon pool
!-----------------------------------------------------------------------
         XV(3)=-HALFDRAFT
         X(3)=(V+1.E0)*XV(3)
         XU(2)=-HALFYMP
         X(1)=XMP
         X(2)=(U-1.E0)*XU(2)
      ELSEIF(IPI.EQ.6) THEN
!-----------------------------------------------------------------------
!                bottom surface fwd of moon pool
!-----------------------------------------------------------------------
         XU(1)=QUARTLEN-HALFXMP
         XV(2)=HALFYMP
         X(1)=(U+1.E0)*XU(1)+XMP
         X(2)=(V+1.E0)*XV(2)
         X(3)=-DRAFT
      ELSEIF(IPI.EQ.7) THEN
!-----------------------------------------------------------------------
!                 free surface of moon pool (optional)
!-----------------------------------------------------------------------
         XU(1)=HALFXMP
         XV(2)=HALFYMP
         X(1)=(U+1.E0)*XU(1)
         X(2)=(V+1.E0)*XV(2)
      ENDIF
 99   RETURN
      END

      SUBROUTINE CYLMP(U,V,IBI,IPI,X,XU,XV,IGDFSCRB,IER)
!-----------------------------------------------------------------------
!
!     COPYRIGHT (C) 2000-7      WAMIT INCORPORATED
!
!-----------------------------------------------------------------------
!
!     Version : 6.4
!
!-----------------------------------------------------------------------
!
!     Source-code file :  GEOMXACT.F
!
!-----------------------------------------------------------------------
!
!     DESCRIPTION :  (IGDEF=-7)
!
!     Subroutine defines the first quadrant of a circular cylinder
!     of radius RADIUS with a concentric moonpoll of radius RADMP < RADIUS
!     RADIUS, DRAFT, RADMP are read from GDF input file
!
!     Patch 1: side
!     Patch 2: bottom
!     Patch 3: inside of moonpool (vertical surface)
!     Patch 4: moonpool free surface (optional)
!
!     Options:
!        Use NPATCH=3 for conventional case
!        Use NPATCH=4 with generalized modes to damp moonpool resonances
!  1/07  INONUMAP is optional input:
!        INONUMAP=0 (default) uniform mapping
!        INONUMAP=1 quadratic on sides, cubic on bottom
!        INONUMAP=2 cubic on sides and bottom
!-----------------------------------------------------------------------
      IMPLICIT NONE
      REAL X(3),XU(3),XV(3),U,V
      INTEGER IPI,IBI,IGDFSCRB,IER
      INTEGER, SAVE :: INONUMAP
      REAL, SAVE :: RADIUS,DRAFT,RADMP,HALFDRAFT,HALFRMP,RADAVG,RADDIF, &
     &              QDR,THREEQDR,HRADDIF,THRADDIF
      REAL THETA,CTHETA,STHETA,SR,THREESRV,ONEPMV
      REAL, PARAMETER :: PIO4=0.785398163397E0
      IF(IPI.EQ.0) THEN
        READ (IGDFSCRB,*,IOSTAT=IER) RADIUS,DRAFT,RADMP
        IF (IER.NE.0) GOTO 99
        INONUMAP=0
        READ (IGDFSCRB,*,END=10,ERR=10) INONUMAP
  10    HALFDRAFT=0.5E0*DRAFT
        QDR=0.25E0*DRAFT
        THREEQDR=0.75E0*DRAFT
        RADAVG=0.5E0*(RADIUS+RADMP)
        RADDIF=0.5E0*(RADIUS-RADMP)
        HRADDIF=0.5E0*RADDIF
        THRADDIF=1.5E0*RADDIF
        HALFRMP=0.5E0*RADMP
      ELSE
         THETA=PIO4*(U+1.E0)
         CTHETA=COS(THETA)
         STHETA=SIN(THETA)
      ENDIF
!-----------------------------------------------------------------------
!  Patch 1: side  V=-1 at free surface, V=+1 at corner
!    If INONUMAP=0 X(3) is a linear function of V
!    If INONUMAP=1 X(3) is quadratic in V with derivative=0 at V=+1
!    If INONUMAP=2 X(3) is a cubic in V with derivatives=0 at V=+1,-1
!-----------------------------------------------------------------------
      IF(IPI.EQ.1) THEN
         X(1)=RADIUS*CTHETA
         X(2)=RADIUS*STHETA
         XU(1)=-PIO4*X(2)
         XU(2)= PIO4*X(1)
         ONEPMV=1.E0+V
         IF(INONUMAP==0) THEN
            X(3)=-HALFDRAFT*ONEPMV
            XV(3)=-HALFDRAFT
         ELSEIF(INONUMAP==1) THEN
            X(3)= QDR*(V-3.E0)*ONEPMV
            XV(3)=HALFDRAFT*(V-1.E0)
         ELSE
            X(3)= QDR*(V-2.E0)*ONEPMV*ONEPMV
            XV(3)=THREEQDR*ONEPMV*(V-1.E0)
         ENDIF
!-----------------------------------------------------------------------
!  Patch 2: bottom
!-----------------------------------------------------------------------
      ELSEIF(IPI.EQ.2) THEN
         IF(INONUMAP==0) THEN
            SR=RADAVG-RADDIF*V
            X(1)=SR*CTHETA
            X(2)=SR*STHETA
            XV(1)=-RADDIF*CTHETA
            XV(2)=-RADDIF*STHETA
         ELSE
            ONEPMV=(1.E0+V)*(1.E0-V)
            THREESRV=-THRADDIF*ONEPMV
            SR=RADAVG-HRADDIF*V*(2.E0+ONEPMV)
            X(1)=SR*CTHETA
            X(2)=SR*STHETA
            XV(1)=THREESRV*CTHETA
            XV(2)=THREESRV*STHETA
         ENDIF
         X(3)=-DRAFT
         XU(1)=-PIO4*X(2)
         XU(2)= PIO4*X(1)

!-----------------------------------------------------------------------
!  Patch 3: side of moonpool  V=-1 at corner, V=+1 at free surface
!-----------------------------------------------------------------------
      ELSEIF(IPI.EQ.3) THEN
         X(1)=RADMP*CTHETA
         X(2)=RADMP*STHETA
         XU(1)=-PIO4*X(2)
         XU(2)= PIO4*X(1)
         ONEPMV=1.E0-V
         IF(INONUMAP==0) THEN
            X(3)=-HALFDRAFT*ONEPMV
            XV(3)=HALFDRAFT
         ELSEIF(INONUMAP==1) THEN
            X(3)=-QDR*(V+3.E0)*ONEPMV
            XV(3)=HALFDRAFT*(V+1.E0)
         ELSE
            X(3)=-QDR*(V+2.E0)*ONEPMV*ONEPMV
            XV(3)=-THREEQDR*ONEPMV*(V+1.E0)
         ENDIF
!-----------------------------------------------------------------------
!  Patch 4: moonpool free surface (optional)
!-----------------------------------------------------------------------
      ELSEIF(IPI.EQ.4) THEN
         SR=HALFRMP*(1.E0-V)
         X(1)=SR*CTHETA
         X(2)=SR*STHETA
         XU(1)=-PIO4*X(2)
         XU(2)= PIO4*X(1)
         XV(1)=-HALFRMP*CTHETA
         XV(2)=-HALFRMP*STHETA
      ENDIF
 99   RETURN
      END

      SUBROUTINE TORUS(U,V,IBI,IPI,X,XU,XV,IGDFSCRB,IER)
!-----------------------------------------------------------------------
!
!     COPYRIGHT (C) 2000      WAMIT INCORPORATED
!
!-----------------------------------------------------------------------
!
!     Version : 6.0
!
!-----------------------------------------------------------------------
!
!     Source-code file :  GEOMXACT.F
!
!-----------------------------------------------------------------------
!
!     DESCRIPTION :  (IGDEF=-8)
!
!     Subroutine defines the torus (IPI=1) by one patch
!     Option to define the moonpool free surface by IPI=2 if NPATCH=2
!     Vertical position is arbitrary (floating or submerged), but
!        generating circle should not be tangent to free surface
!
!     Input from GDF file:
!         RCIRC     radius of generating circles (sections of torus)
!         RAXIS     radius of axis (center of generating circles)
!         ZAXIS     vertical coordinate of axis (negative if submerged)
!
!     Restrictions:
!         RAXIS > RCIRC   (open hole in center of torus)
!         |ZAXIS| < RCIRC (floating torus) or ZAXIS < -RCIRC (submerged)
!-----------------------------------------------------------------------
      IMPLICIT NONE
      REAL X(3),XU(3),XV(3),U,V
      INTEGER IPI,IBI,IGDFSCRB,IER
      REAL, SAVE :: RAXIS,RCIRC,ZAXIS,PSI1,HALFRMP
      REAL THETA,CTHETA,STHETA,PSI,RCOSPSI,RSINPSI,SR,AZ,HBWL
!-----------------------------------------------------------------------
!   PIO4=pi/4
!-----------------------------------------------------------------------
      REAL, PARAMETER :: PIO4=0.78539816E0, PI=3.14159263E0
      IF(IPI.EQ.0) THEN
        READ (IGDFSCRB,*,IOSTAT=IER) RCIRC,RAXIS,ZAXIS
        IF (IER.NE.0) GOTO 99
!-----------------------------------------------------------------------
!  First if block is for floating case, second submerged
!  HBWL = half beam of section at waterline
!  PSI1 = angle to waterline from below axis (if axis in Z=0, PSI1=pi/2,
!                                             if submerged, PSI1=pi)
!~ HALFRMP = half * radius of moon pool
!-----------------------------------------------------------------------
        AZ=ABS(ZAXIS)
        IF (AZ < RCIRC) THEN
          HBWL=SQRT(RCIRC*RCIRC-AZ*AZ)
          PSI1=ATAN2(HBWL,ZAXIS)
          HALFRMP=0.5E0*(RAXIS-HBWL)
        ELSE
          PSI1=PI
        ENDIF
      ELSEIF(IPI.EQ.1) THEN
!-----------------------------------------------------------------------
!  THETA = polar angle about z-axis, 0<THETA<PI/4 for the quadrant:
!     U=-1, THETA=0, Y=0 plane;  U=+1, THETA=pi/2, X=0 plane
!  PSI = angle around semi-circle below free surface, -PSI1<PSI< PSI1
!     V=-1, PSI=-PSI1  is the outer waterline (if floating)
!     V= 0, PSI=0      is the point of maximum draft
!     V=+1, PSI=+PSI1  is the inner waterline (if floating)
!-----------------------------------------------------------------------
         THETA=PIO4*(U+1.E0)
         CTHETA=COS(THETA)
         STHETA=SIN(THETA)
         PSI=PSI1*V
         RCOSPSI=RCIRC*COS(PSI)
         RSINPSI=RCIRC*SIN(PSI)
         SR=RAXIS-RSINPSI
         X(1)=SR*CTHETA
         X(2)=SR*STHETA
         X(3)=-RCOSPSI+ZAXIS
         XU(1)=-PIO4*X(2)
         XU(2)= PIO4*X(1)
         XV(1)=-PSI1*RCOSPSI*CTHETA
         XV(2)=-PSI1*RCOSPSI*STHETA
         XV(3)= PSI1*RSINPSI
!-----------------------------------------------------------------------
!  Patch 2: moonpool free surface
!-----------------------------------------------------------------------
      ELSEIF(IPI.EQ.2) THEN
         SR=HALFRMP*(1.E0-V)
         THETA=PIO4*(U+1.E0)
         CTHETA=COS(THETA)
         STHETA=SIN(THETA)
         X(1)=SR*CTHETA
         X(2)=SR*STHETA
         XU(1)=-PIO4*X(2)
         XU(2)= PIO4*X(1)
         XV(1)=-HALFRMP*CTHETA
         XV(2)=-HALFRMP*STHETA
      ENDIF
 99   RETURN
      END

      SUBROUTINE TLP(U,V,IBI,IPI,X,XU,XV,IGDFSCRB,IER)
!-----------------------------------------------------------------------
!
!     COPYRIGHT (C) 2000      WAMIT INCORPORATED
!
!-----------------------------------------------------------------------
!
!     Version : 6.0
!
!-----------------------------------------------------------------------
!
!     Source-code file :  GEOMXACT.F
!
!-----------------------------------------------------------------------
!
!     DESCRIPTION :  (IGDEF=-9)
!
!     Subroutine defines the first quadrant of a TLP with equal spacing
!       between the four columns.
!       Columns are circular cylinders, pontoons are rectangular
!
!     RADIUS: Radius of the cylinder
!     DRAFT: Draft of TLP
!     HSPACE: .5*Space between the centers of adjacent columns
!     WIDTH: Pontoon width
!     HEIGHT: Pontoon height
!
!  note READ statements below read RADIUS,DRAFT,HSPACE on 1st line
!                                  WIDTH, HEIGHT       on 2nd line
!
!  RESTRICTIONS ON THE GEOMETRY:
!
!  The bottom of the pontoons and bottom of the columns must be in the same
!      horizontal plane (Z=-DRAFT) and the height of the pontoons must be
!      less than the draft (pontoons must be completely submerged)
!
!  The width of the pontoons must be less than or equal to SQRT(2)*RADIUS,
!     so that the pontoons do not intersect outside the column radius.
!     In the case where the pontoons intersect precisely on the column
!     (WIDTH = SQRT(2)*RADIUS) there are a total of NPATCH=11 patches, as
!     defined below.  In the case where there is a space between the pontoons
!     on the column, NPATCH=12.
!
!  Mod 9/03 Add interior free surface as patch number 13
!
!-----------------------------------------------------------------------
      IMPLICIT NONE
      REAL X(3),XU(3),XV(3),U,V
      INTEGER IPI,IBI,IGDFSCRB,IER
      REAL, PARAMETER :: PI=3.141592654E0
      REAL, PARAMETER :: HALF=0.5E0
!-----------------------------------------------------------------------
!  Indices for (X,Y) coordinates on two pontoons
!-----------------------------------------------------------------------
      INTEGER, PARAMETER :: IXY(1:8)=(/1,1,1,1,2,2,2,2/)
      INTEGER, PARAMETER :: IYX(1:8)=(/2,2,2,2,1,1,1,1/)
      REAL  AP,T,COST,SINT,TH,COSTH,SINTH,R,RCOSTH,RSINTH,VP1,RS
      REAL, SAVE :: XFAC(1:8),XUFAC(1:8),X3(1:6)
      REAL, SAVE :: AU(10:12),AV(10:12),BU(10:12),BV(10:12)
      REAL, SAVE :: RADIUS,DRAFT,HSPACE,WIDTH,HEIGHT,HRADANG,           &
     &              HWIDTH,ANG,XL,HALFHT,HALFXL,HALFRAD,HALFHSP
!-----------------------------------------------------------------------
!  Read dimensions from GDF file, pre-evaluate geometric parameters
!-----------------------------------------------------------------------
      IF(IPI.EQ.0) THEN
        READ (IGDFSCRB,*,IOSTAT=IER) RADIUS, DRAFT, HSPACE
        READ (IGDFSCRB,*,IOSTAT=IER) WIDTH, HEIGHT
        IF (IER.NE.0) GOTO 99
        HWIDTH=HALF*WIDTH
        ANG=ATAN2(HWIDTH,SQRT(RADIUS*RADIUS-HWIDTH*HWIDTH))
        XL=HSPACE - RADIUS*COS(ANG)
        HALFHT=HALF*HEIGHT
        HALFXL=HALF*XL
        HALFRAD=HALF*RADIUS
        HALFHSP=HALF*HSPACE
        HRADANG=HALFRAD*ANG
!-----------------------------------------------------------------------
!  Factors for top and bottom of pontoons
!-----------------------------------------------------------------------
        XFAC(1)= RADIUS
        XFAC(2)=-RADIUS
        XFAC(5)=-RADIUS
        XFAC(6)= RADIUS
        XUFAC(1)=RADIUS*ANG
        XUFAC(2)=-XUFAC(1)
        XUFAC(5)=-XUFAC(1)
        XUFAC(6)= XUFAC(1)
        X3(1)=-DRAFT+HEIGHT
        X3(2)=-DRAFT
        X3(5)=-DRAFT+HEIGHT
        X3(6)=-DRAFT
!-----------------------------------------------------------------------
!  Factors for sides of pontoons
!-----------------------------------------------------------------------
        XFAC(3)= HWIDTH
        XFAC(4)=-HWIDTH
        XFAC(7)= HWIDTH
        XFAC(8)=-HWIDTH
        XUFAC(3)=-HALFHT
        XUFAC(4)= HALFHT
        XUFAC(7)= HALFHT
        XUFAC(8)=-HALFHT
!-----------------------------------------------------------------------
!  Factors for side of column
!-----------------------------------------------------------------------
        AU(10)=PI
        AU(11)=0.75E0*PI-ANG
        AU(12)=0.25E0*PI-ANG
        BU(10)=PI
        BU(11)=0.25E0*PI
        BU(12)=1.25E0*PI
        AV(10)=-HALF*(DRAFT-HEIGHT)
        AV(11)=-HALFHT
        AV(12)=-HALFHT
        BV(10)=AV(10)
        BV(11)=-DRAFT+HALFHT
        BV(12)=BV(11)
!-----------------------------------------------------------------------
!   Patch 1: top of pontoon parallel to x-axis
!   Patch 2: bottom of pontoon parallel to x-axis
!   Patch 5: top of pontoon parallel to y-axis
!   Patch 6: bottom of pontoon parallel to y-axis
!-----------------------------------------------------------------------
      ELSEIF ((IPI==1).OR.(IPI==2).OR.(IPI==5).OR.(IPI==6)) THEN
         T=ANG*U
         COST=COS(T)
         SINT=SIN(T)
         AP=-HALFRAD*COST+HALFHSP
         VP1=V+1.E0
         X(IXY(IPI))=AP*VP1
         X(IYX(IPI))=HSPACE+XFAC(IPI)*SINT
         XU(IXY(IPI))=HRADANG*SINT*VP1
         XU(IYX(IPI))=XUFAC(IPI)*COST
         XV(IXY(IPI))=AP
         X(3)=X3(IPI)
!-----------------------------------------------------------------------
!   Patch 3: outside of pontoon parallel to x-axis
!   Patch 4: inside of pontoon parallel to x-axis
!   Patch 7: outside of pontoon parallel to y-axis
!   Patch 8: inside of pontoon parallel to y-axis
!-----------------------------------------------------------------------
      ELSEIF ((IPI==3).OR.(IPI==4).OR.(IPI==7).OR.(IPI==8)) THEN
         X(IXY(IPI))=HALFXL*(V+1.)
         X(IYX(IPI))=HSPACE+XFAC(IPI)
         X(3)=XUFAC(IPI)*U+BV(11)
         XU(3)=XUFAC(IPI)
         XV(IXY(IPI))=HALFXL
!-----------------------------------------------------------------------
!   Patch 9:  bottom of column
!-----------------------------------------------------------------------
      ELSEIF(IPI.EQ.9) THEN
         TH=PI*(U+1.)
         R=-HALFRAD*(V-1.)
         COSTH=COS(TH)
         SINTH=SIN(TH)
         RCOSTH=R*COSTH
         RSINTH=R*SINTH
         X(1)=RCOSTH+HSPACE
         X(2)=RSINTH+HSPACE
         X(3)=-DRAFT
         XU(1)=-RSINTH*PI
         XU(2)=RCOSTH*PI
         XV(1)=-HALFRAD*COSTH
         XV(2)=-HALFRAD*SINTH
!-----------------------------------------------------------------------
!   Patch 10:  side of column above pontoons
!   Patch 11:  side of column outside pontoons
!   Patch 12:  side of column between pontoons
!-----------------------------------------------------------------------
      ELSEIF ((IPI==10).OR.(IPI==11).OR.(IPI==12)) THEN
         TH=AU(IPI)*U+BU(IPI)
         RCOSTH=RADIUS*COS(TH)
         RSINTH=RADIUS*SIN(TH)
         X(1)=RCOSTH+HSPACE
         X(2)=RSINTH+HSPACE
         X(3)=AV(IPI)*V+BV(IPI)
         XU(1)=-AU(IPI)*RSINTH
         XU(2)=AU(IPI)*RCOSTH
         XV(3)=AV(IPI)
      ELSE
         TH=PI*(U+1.)
         RS=HALFRAD*(V+1.)
         COSTH=COS(TH)
         SINTH=SIN(TH)
         X(1)= RS*COSTH+HSPACE
         X(2)= RS*SINTH+HSPACE
         XU(1)=-PI*RS*SINTH
         XU(2)= PI*RS*COSTH
         XV(1)= HALFRAD*COSTH
         XV(2)= HALFRAD*SINTH
      ENDIF
 99   RETURN
      END

      SUBROUTINE SEMISUB(U,V,IBI,IPI,X,XU,XV,IGDFSCRB,IER)
!-----------------------------------------------------------------------
!
!     COPYRIGHT (C) 2000      WAMIT INCORPORATED
!
!-----------------------------------------------------------------------
!
!     Version : 6.0
!
!-----------------------------------------------------------------------
!
!     Source-code file :  GEOMXACT.F
!
!-----------------------------------------------------------------------
!
!     DESCRIPTION :  (IGDEF=-10)
!
!     Subroutine describes one quadrant of a semi-sub with two parallel
!     pontoons and equally-spaced cylindrical columns. The pontoon
!     is rectangular with semi-circular ends. All columns have the same
!     radius.  The complete semi-sub is symmetric about x=y=0.  The
!     subroutine describes only the first quadrant (x>0, y>0).
!
!   The geometric input parameters:
!     XL    total length of pontoons (including halfbody x<0)
!     Y1,Y2 y-coordinates of two sides of pontoon  (Y2>Y1>0)
!     Z1,Z2 z-coordinates of bottom,top of pontoon (0>=Z2>Z1)
!     NCOL  number of columns (identical, symmetric about centerline)
!           on one pontoon
!     DCOL  distance between column axes (all equally spaced)
!     RCOL  column radius (all assumed equal)
!
!  note READ statements below read XL,Y1,Y2,Z1,Z2      on 1st line
!                                  DCOL, RCOL, NCOL    on 2nd line
!
!   Additional resctrictions on input parameters:
!     2*RCOL < (Y2-Y1)  (column diameter less than width of pontoon)
!     DCOL > (Y2-Y1)    (column spacing greater than width of pontoon)
!       XL >= (NCOL-1)*DCOL + (Y2-Y1)   (pontoon length must be greater
!        than or equal to the total length occupied by the columns and
!        annular rings described below.
!
!   General description of patches:
!     a) one patch for the bottom of the pontoon
!     b) one patch for the sides and end of the pontoon
!     c)  one patch on each column
!     d) annular ring around the base of each column extending out to
!         the sides of the pontoons
!     e) one patch on pontoon deck between adjacent annular rings
!     f) one patch on pontoon deck at end (if required)
!
!   If Z1=0 and NPATCH=2, the pontoons are floating on the free surface.
!
!   In the normal case where the pontoons are submerged, there are four
!   different cases, depending on the number of columns and length of
!   the pontoons:
!
!   Odd number of columns (NCOL = 1,3,5, ...)
!
!     If    XL = (NCOL-1)*DCOL + (Y2-Y1),   NPATCH=3*((NCOL-1)/2)+4
!
!     If    XL > (NCOL-1)*DCOL + (Y2-Y1),   NPATCH=3*((NCOL-1)/2)+5
!
!      IPI=1       pontoon bottom
!      IPI=2       pontoon sides and end
!      IPI=3       half-surface (x>0) of the column at x=0
!      IPI=4       half-annulus around the base of the column at x=0
!      IPI=5       pontoon deck adjoining the half-annulus
!      IPI=6       surface of 1st column at x>0
!      IPI=7       annulus around 1st column at x>0
!         .
!         .
!      NPATCH      annulus around last column or deck forward of annulus
!
!
!   Even number of columns (NCOL = 2,4,6, ...)
!
!     If    XL = (NCOL-1)*DCOL + (Y2-Y1),   NPATCH=3*(NCOL/2)+2
!
!     If    XL > (NCOL-1)*DCOL + (Y2-Y1),   NPATCH=3*(NCOL/2)+3
!
!      IPI=1       pontoon bottom
!      IPI=2       pontoon sides and end
!      IPI=3       half-surface (x>0) of the deck from x=0 to 1st annulus
!      IPI=4       surface of 1st column at x>0
!      IPI=5       annulus around 1st column at x>0
!      IPI=6       deck forward of 1st annulus
!         .
!         .
!      NPATCH      annulus around last column or deck forward of annulus
!
!   Columns are numbered in increasing order w.r.t. x
!
!---------------------------------------------------------------------
      IMPLICIT NONE
      REAL X(3),XU(3),XV(3),U,V
      INTEGER IPI,IBI,IGDFSCRB,IER
      REAL, SAVE :: Y1,Y2,Z1,Z2,DCOL,RCOL,                              &
     &              HYP,HZP,HYM,HZM,HDCOL,HZ2,XCOL1,HXDECK,             &
     &              ARCSIDE,ARCRAT,XUSIDE,THETAU
      INTEGER, SAVE :: IMODNCOL,NPATCHM
      REAL, PARAMETER :: PIO4=0.78539816E0, PIO2=1.57079633E0,          &
     &      PI=3.1415926E0, ZERO=0.0E0, HALF=0.5E0, ONE=1.E0, TWO=2.E0
      REAL ARCEND,THETA,RCOS,RSIN,XXH,HXXH,COST,SINT,RDIF,RA,XCOL,XL,HXL
      INTEGER IPIL,IMODIPIL,NCOL
      IF (IPI.EQ.0) THEN
         READ(IGDFSCRB,*,IOSTAT=IER) XL,Y1,Y2,Z1,Z2
         READ(IGDFSCRB,*,IOSTAT=IER) DCOL,RCOL,NCOL
         IF (IER.NE.0) GOTO 99
         HXL=HALF*XL
         HYP=HALF*(Y1+Y2)
         HYM=HALF*(Y2-Y1)
         HZP=HALF*(Z1+Z2)
         HZM=HALF*(Z2-Z1)
         HZ2=HALF*Z2
         HDCOL=HALF*DCOL
         HXDECK=HALF*(HXL-(NCOL-1)*HDCOL-HYM)
!-----------------------------------------------------------------------
!  Arc lengths of pontoon side, end, and ratio end/total
!-----------------------------------------------------------------------
         ARCSIDE=HXL-HYM
         ARCEND=PI*HYM
         ARCRAT=ARCEND/(ARCEND+TWO*ARCSIDE)
!-----------------------------------------------------------------------
!  Derivatives of X,theta w.r.t. U on side, end
!-----------------------------------------------------------------------
         XUSIDE=ARCSIDE/(ONE-ARCRAT)
         THETAU=PIO2/ARCRAT
!-----------------------------------------------------------------------
!  IMODNCOL = 0,1 for even,odd number of columns
!  NPATCHM = max value of NPATCH (index of extra patch on deck)
!  XCOL1 = x-coordinate of column before the 1st complete column
!-----------------------------------------------------------------------
         IMODNCOL=MOD(NCOL,2)
         NPATCHM=3*(NCOL/2)+3+2*IMODNCOL
         IF (IMODNCOL == 0) THEN
           XCOL1=-HDCOL
         ELSE
           XCOL1=ZERO
         ENDIF
      ELSEIF (IPI.EQ.1) THEN
!-----------------------------------------------------------------------
!     IPI=1: patch on bottom of pontoon
!-----------------------------------------------------------------------
          THETA=PIO2*V
          RCOS=HYM*COS(THETA)
          RSIN=HYM*SIN(THETA)
          X(2)=RSIN+HYP
          XXH=ARCSIDE+RCOS
          HXXH=HALF*XXH
          X(1)=HXXH*(U+ONE)
          X(3)=Z1
          XU(1)=HXXH
          XV(1)=-RSIN*PIO4*(U+ONE)
          XV(2)=RCOS*PIO2
      ELSEIF (IPI.EQ.2) THEN
!-----------------------------------------------------------------------
!   IPI=2: patch on sides and end of pontoon subdivided in 3 segments:
!     Inboard side of pontoon   (-1 < u < -ARCRAT)
!     End of pontoon            (-ARCRAT < u < ARCRAT)
!     Outboard side of pontoon   (ARCRAT < u < 1)
!-----------------------------------------------------------------------
         IF(U.LE.-ARCRAT) THEN
            X(1)=XUSIDE*(U+ONE)
            X(2)=Y1
            XU(1)=XUSIDE
         ELSEIF(U.LT.ARCRAT) THEN
            THETA=THETAU*U
            RCOS=HYM*COS(THETA)
            RSIN=HYM*SIN(THETA)
            X(1)=RCOS+ARCSIDE
            X(2)=RSIN+HYP
            XU(1)=-THETAU*RSIN
            XU(2)= THETAU*RCOS
         ELSE
            X(1)=XUSIDE*(ONE-U)
            X(2)=Y2
            XU(1)=-XUSIDE
         ENDIF
         X(3)=HZP-V*HZM
         XV(3)=-HZM
      ELSEIF ((IPI==3).AND.(IMODNCOL==1)) THEN
!-----------------------------------------------------------------------
!    Half-column at x=0  (ICOL = odd)
!-----------------------------------------------------------------------
         THETA=PIO2*U
         RCOS=RCOL*COS(THETA)
         RSIN=RCOL*SIN(THETA)
         X(1)=RCOS
         X(2)=RSIN+HYP
         XU(1)=-PIO2*RSIN
         XU(2)= PIO2*RCOS
         X(3)=HZ2*(V+ONE)
         XV(3)=HZ2
      ELSEIF (IPI==3) THEN
!-----------------------------------------------------------------------
!    deck from x=0 to first annulus  (ICOL = even)
!-----------------------------------------------------------------------
         THETA=PIO2*V
         RCOS=HYM*COS(THETA)
         RSIN=HYM*SIN(THETA)
         X(2)=-RSIN+HYP
         XXH=HDCOL-RCOS
         HXXH=HALF*XXH
         X(1)=HXXH*(U+ONE)
         X(3)=Z2
         XU(1)=HXXH
         XV(1)=RSIN*PIO4*(U+ONE)
         XV(2)=-RCOS*PIO2
      ELSEIF ((IPI==4).AND.(IMODNCOL==1)) THEN
!-----------------------------------------------------------------------
!    Half-annulus at x=0  (ICOL = odd)
!-----------------------------------------------------------------------
         THETA=PIO2*U
         COST=COS(THETA)
         SINT=SIN(THETA)
         RDIF=HALF*(HYM-RCOL)
         RA=RDIF*(V+ONE)+RCOL
         RCOS=RA*COST
         RSIN=RA*SINT
         X(1)=RCOS
         X(2)=RSIN+HYP
         X(3)=Z2
         XU(1)=-PIO2*RSIN
         XU(2)= PIO2*RCOS
         XV(1)=RDIF*COST
         XV(2)=RDIF*SINT
      ELSEIF (IPI < NPATCHM) THEN
!-----------------------------------------------------------------------
!    Symmetric patches on columns, annulae, deck between annulae:
!       IMODIPIL = 1: columns centered at XCOL
!       IMODIPIL = 2: annulae centered at XCOL
!       IMODIPIL = 0: deck centered between annulae XCOL-DCOL,XCOL
!-----------------------------------------------------------------------
         IPIL=IPI-2*IMODNCOL
         XCOL=(IPIL/3)*DCOL+XCOL1
         IMODIPIL=MOD(IPIL,3)
         IF (IMODIPIL == 1) THEN
            THETA=PI*(U+ONE)
            RCOS=RCOL*COS(THETA)
            RSIN=RCOL*SIN(THETA)
            X(1)=RCOS+XCOL
            X(2)=RSIN+HYP
            XU(1)=-PI*RSIN
            XU(2)= PI*RCOS
            X(3)= HZ2*(V+ONE)
            XV(3)=HZ2
         ELSEIF (IMODIPIL == 2) THEN
            THETA=PI*(U+ONE)
            COST=COS(THETA)
            SINT=SIN(THETA)
            RDIF=HALF*(HYM-RCOL)
            RA=RDIF*(V+ONE)+RCOL
            RCOS=RA*COST
            RSIN=RA*SINT
            X(1)=RCOS+XCOL
            X(2)=RSIN+HYP
            X(3)=Z2
            XU(1)=-PI*RSIN
            XU(2)= PI*RCOS
            XV(1)=RDIF*COST
            XV(2)=RDIF*SINT
         ELSE
            THETA=PIO2*V+PI
            RCOS=HYM*COS(THETA)
            RSIN=HYM*SIN(THETA)
            X(2)=RSIN+HYP
            XXH=HDCOL+RCOS
            X(1)=XXH*U+XCOL-HDCOL
            X(3)=Z2
            XU(1)=XXH
            XV(1)=-RSIN*PIO2*U
            XV(2)=RCOS*PIO2
         ENDIF
      ELSE
!-----------------------------------------------------------------------
!   Pontoon extends beyond last annulus.  Extra patch for deck at end.
!-----------------------------------------------------------------------
         THETA=PIO2*V
         RCOS=HYM*COS(THETA)
         RSIN=HYM*SIN(THETA)
         X(1)=ARCSIDE+RCOS+HXDECK*(U-ONE)
         X(2)=-RSIN+HYP
         X(3)=Z2
         XU(1)=HXDECK
         XV(1)=-RSIN*PIO2
         XV(2)=-RCOS*PIO2
      ENDIF
 99   RETURN
      END

      SUBROUTINE FPSO(U,V,IBI,IPI,X,XU,XV,IGDFSCRB,IER)
!-----------------------------------------------------------------------
!
!     COPYRIGHT (C) 2000      WAMIT INCORPORATED
!
!-----------------------------------------------------------------------
!
!     Version : 6.0
!
!-----------------------------------------------------------------------
!
!     Source-code file :  GEOMXACT.F
!
!-----------------------------------------------------------------------
!
!     DESCRIPTION :  (IGDEF=-11)
!
!     Subroutine describes one side of an FPSO hull with elliptic bow,
!     rectangular mid-body, prismatic stern.
!
!   The geometric input parameters:
!     XBOW  length of bow (semi-axis of ellipse)
!     XMID  length of midbody
!     XAFT  length of prismatic afterbody
!     HBEAM half beam
!     HTRANSOM half width of transom
!     DRAFT draft
!     DTRANSOM depth of transom
!
!   General description of patches:
!
!    IPI=1 patch for the bottom including bow and midbody
!    IPI=2 patch for bow
!    IPI=3 patch for side
!    IPI=4 patch on transom
!    IPI=5 sloping bottom on prismatic stern
!    IPI=6 sloping side on prismatic stern
!
!   If NPATCH=4 the prismatic patches are omitted and the following
!   dimensions are assigned: XAFT=0.0, HTRANSOM=HBEAM, DTRANSOM=DRAFT
!
!
!-------------------------------------------------------------------------
      IMPLICIT NONE
      REAL X(3),XU(3),XV(3),U,V
      INTEGER IPI,IBI,IGDFSCRB,IER
      REAL, SAVE :: XBOW,XMID,XAFT,HBEAM,HTRANSOM,DRAFT,DTRANSOM,       &
     &              XL,X1,X2,X3,X4,XSIDMID,XUSIDE,HDRAFT,HALFDT,        &
     &              HALFHT,HALFXAFT,XAFTMID,ZBAFTMID,ZBAFTFAC,          &
     &              QHT,QHB,FACUB,YBAFTMID,YSAFTFAC,ZSAFTMID,QDT,QDR,   &
     &              FACUS,THETA,COSTH,SINTH,XXH,HXXH,YSAFTMID
      REAL, PARAMETER :: PI=3.1415926E0, HALF=0.5E0, ONE=1.E0
      REAL, PARAMETER :: PIO8=PI*0.125E0, PIO4=PI*.25E0
      IF (IPI.EQ.0) THEN
         READ(IGDFSCRB,*,IOSTAT=IER) XBOW,XMID,XAFT
         READ(IGDFSCRB,*,IOSTAT=IER) HBEAM,HTRANSOM
         READ(IGDFSCRB,*,IOSTAT=IER) DRAFT,DTRANSOM
         IF (IER.NE.0) GOTO 99
         XL=XBOW+XMID+XAFT
         X4=HALF*XL
         X1=-X4
         X2=X1+XAFT
         X3=X2+XMID
         XSIDMID=HALF*(X2+X3)
         XUSIDE=HALF*(X3-X2)
         HDRAFT=HALF*DRAFT
         HALFDT=HALF*DTRANSOM
         HALFHT=HALF*HTRANSOM
         HALFXAFT=HALF*XAFT
         XAFTMID=HALF*(X1+X2)
         ZBAFTMID=-HALF*(DRAFT+DTRANSOM)
         ZBAFTFAC=-HALF*(DRAFT-DTRANSOM)
         QHT=0.25*HTRANSOM
         QHB=0.25*HBEAM
         FACUB=QHB-QHT
         YBAFTMID=QHB+QHT
         YSAFTMID=HALF*(HBEAM+HTRANSOM)
         YSAFTFAC=HALF*(HBEAM-HTRANSOM)
         ZSAFTMID=0.25E0*(DTRANSOM+DRAFT)
         QDT=0.25*DTRANSOM
         QDR=0.25*DRAFT
         FACUS=QDR-QDT
      ELSEIF (IPI.EQ.1) THEN
!-----------------------------------------------------------------------
!     IPI=1: patch on bottom
!-----------------------------------------------------------------------
          THETA=PIO4*(V+ONE)
          COSTH=COS(THETA)
          SINTH=SIN(THETA)
          X(2)=HBEAM*SINTH
          XXH=XMID+COSTH*XBOW
          HXXH=HALF*XXH
          X(1)=X2+HXXH*(ONE+U)
          X(3)=-DRAFT
          XU(1)=HXXH
          XV(1)=-SINTH*PIO8*XBOW*(U+ONE)
          XV(2)=COSTH*HBEAM*PIO4
      ELSEIF (IPI.EQ.2) THEN
!-----------------------------------------------------------------------
!   IPI=2: patch on bow
!-----------------------------------------------------------------------
          THETA=PIO4*(U+ONE)
          COSTH=COS(THETA)
          SINTH=SIN(THETA)
          X(1)=X3+XBOW*COSTH
          X(2)=HBEAM*SINTH
          XU(1)=-PIO4*XBOW*SINTH
          XU(2)= PIO4*HBEAM*COSTH
          X(3)=-HDRAFT-V*HDRAFT
          XV(3)=-HDRAFT
      ELSEIF (IPI.EQ.3) THEN
!-----------------------------------------------------------------------
!   IPI=3: patch on side
!-----------------------------------------------------------------------
         X(1)=XSIDMID-XUSIDE*U
         X(2)=HBEAM
         XU(1)=-XUSIDE
         X(3)=-HDRAFT-V*HDRAFT
         XV(3)=-HDRAFT
      ELSEIF (IPI==4) THEN
!-----------------------------------------------------------------------
!    Transom
!-----------------------------------------------------------------------
         XV(3)=-HALFDT
         X(3)=(V+1.E0)*XV(3)
         XU(2)=-HALFHT
         X(1)=X1
         X(2)=(U-1.E0)*XU(2)
      ELSEIF (IPI==5) THEN
!-----------------------------------------------------------------------
!    Sloping bottom aft
!-----------------------------------------------------------------------
         X(1)=HALFXAFT*U+XAFTMID
         X(2)=(V+ONE)*(FACUB*U+YBAFTMID)
         X(3)=ZBAFTMID+U*ZBAFTFAC
         XU(1)=HALFXAFT
         XU(2)=FACUB*(V+ONE)
         XV(2)=FACUB*U+YBAFTMID
         XU(3)=ZBAFTFAC
      ELSE
!-----------------------------------------------------------------------
!    Sloping side aft
!-----------------------------------------------------------------------
         X(1)=HALFXAFT*U+XAFTMID
         X(2)=YSAFTMID+U*YSAFTFAC
         X(3)=(V-ONE)*(FACUS*U+ZSAFTMID)
         XU(1)=HALFXAFT
         XU(2)=YSAFTFAC
         XU(3)=(V-ONE)*FACUS
         XV(3)=FACUS*U+ZSAFTMID
      ENDIF
 99   RETURN
      END
      
      SUBROUTINE SPAR(U,V,IBI,IPI,X,XU,XV,IGDFSCRB,IER)
!-----------------------------------------------------------------------
!
!     COPYRIGHT (C) 2002-2003      WAMIT INCORPORATED
!
!-----------------------------------------------------------------------
!
!     Version : 6.2
!
!-----------------------------------------------------------------------
!
!     Source-code file :  GEOMXACT.F
!
!-----------------------------------------------------------------------
!
!     DESCRIPTION :  (IGDEF=-12)
!
!     Subroutine defines a spar with strakes
!
!     Patch 1 : Side
!     Patch 2 : Strake surface
!    (Patch 3): Strake surface on opposite side of Patch 2 (Thickness>0)
!    (Patch 4): Strake tip (Thickness>0)
!    (Patch 5): Strake bottom (Thickness>0)
!       .
!       .
!     Patch NP+1: Bottom (NP=2*Number of strakes when thickness < 0)
!                        (NP=5*Number of strakes when thickness > 0)
!     Patch NP+2: Interior surface (IRR >0. optional)
!     Patch LP+1: Moonpool side wall (optional) (LP=NP+2 when IRR>0)
!                                               (LP=NP+1 when IRR=0)
!     Patch LP+2: Moonpool free surface (NEWMDS > 0,optional)
!
!     GDF Example A-1) Spar with 3 zero-thickness strakes
!
!      TEST21 SPAR with three strakes IGDEF=-12
!      18. 9.80665  ULEN GRAV
!      0  0        ISX  ISY
!      7  -12       NPATCH  IGDEF
!      npatch_dipole = 3
!      ipatch_dipole = 2 4 6
!      5
!      18. 200.              RADIUS DRAFT
!      3.7 0. 1. 3           WIDTH THICKNESS TWIST NSTRAKE
!      0                     IRRFRQ
!      0  0.                 IMOONPOOL, RADIUSMP
!      0                     IMPGEN
!
!     GDF Example A-2) Spar with 3 zero-thickness strakes
!                      with irregular frequency removal option
!
!      TEST21 SPAR with three strakes IGDEF=-12
!      18. 9.80665  ULEN GRAV
!      0  0        ISX  ISY
!      8  -12       NPATCH  IGDEF
!      npatch_dipole = 3
!      ipatch_dipole = 2 4 6
!      5
!      18. 200.              RADIUS DRAFT
!      3.7 0. 1. 3           WIDTH THICKNESS TWIST NSTRAKE
!      0                     IRRFRQ
!      0  0.                 IMOONPOOL, RADIUSMP
!      0                     IMPGEN
!
!     GDF Example A-3) Spar with 3 zero-thickness strakes
!                      with a circular moon pool
!
!      TEST21 SPAR with three strakes IGDEF=-12
!      18. 9.80665  ULEN GRAV
!      0  0        ISX  ISY
!      8  -12       NPATCH  IGDEF
!      npatch_dipole = 3
!      ipatch_dipole = 2 4 6
!      5
!      18. 200.              RADIUS DRAFT
!      3.7 0. 1. 3           WIDTH THICKNESS TWIST NSTRAKE
!      0                     IRRFRQ
!      1  3.                 IMOONPOOL, RADIUSMP
!      0                     IMPGEN
!
!     GDF Example A-4) Spar with 3 zero-thickness strakes
!                      with irregular frequency removal option
!                      with a circular moon pool
!
!      TEST21 SPAR with three strakes IGDEF=-12
!      18. 9.80665  ULEN GRAV
!      0  0        ISX  ISY
!      9  -12       NPATCH  IGDEF
!      npatch_dipole = 3
!      ipatch_dipole = 2 4 6
!      5
!      18. 200.              RADIUS DRAFT
!      3.7 0. 1. 3           WIDTH THICKNESS TWIST NSTRAKE
!      1                     IRRFRQ
!      1  3.                 IMOONPOOL, RADIUSMP
!      0                     IMPGEN
!
!     GDF Example A-5) Spar with 3 zero-thickness strakes
!                      with a circular moon pool
!                      with generalized modes on moon pool free surface
!
!
!      TEST21 SPAR with three strakes IGDEF=-12
!      18. 9.80665  ULEN GRAV
!      0  0        ISX  ISY
!      8  -12       NPATCH  IGDEF
!      npatch_dipole = 3
!      ipatch_dipole = 2 4 6
!      5
!      18. 200.              RADIUS DRAFT
!      3.7 0. 1. 3           WIDTH THICKNESS TWIST NSTRAKE
!      0                     IRRFRQ
!      1  3.                 IMOONPOOL, RADIUSMP
!      1                     IMPGEN
!
!     GDF Example A-6) Spar with 3 zero-thickness strakes
!                      with irregular frequency removal option
!                      with a circular moon pool
!                      with generalized modes on moon pool free surface
!
!      TEST21 SPAR with three strakes IGDEF=-12
!      18. 9.80665  ULEN GRAV
!      0  0        ISX  ISY
!      9  -12       NPATCH  IGDEF
!      npatch_dipole = 3
!      ipatch_dipole = 2 4 6
!      5
!      18. 200.              RADIUS DRAFT
!      3.7 0. 1. 3           WIDTH THICKNESS TWIST NSTRAKE
!      1                     IRRFRQ
!      1  3.                 IMOONPOOL, RADIUSMP
!      1                     IMPGEN
!
!
!
!
!
!
!-----------------------------------------------------------------------
      IMPLICIT NONE
      REAL X(3),XU(3),XV(3),U,V
      INTEGER IPI,IBI,IGDFSCRB,IER
!-----------------------------------------------------------------------
!     R=RADIUS     radius of spar
!     D=DRAFT      draft
!     N=NSTRAKE    number of strakes
!     W=WIDTH      width of the strakes
!     T=THICKNESS  thickness of the strakes
!     TN=TURN      twisting of strakes at the bottom relative to the top
!     IM=IMOONPOOL IM=0 no moonpool, IM=1 with circular moonpool
!     IR=IRRFRQ    IR=1 interior free surface, IR=0 no interior surface
!     IG=IMPGEN    IG=1 moonpool free surface, IG=0 no moonpool surface
!     RM=RADIUSMP  moonpool radius
!
!     RADIUS,DRAFT
!     WIDTH,THICKNESS,TURN,NSTRAKE
!     IRRFRQ
!     IMOONPOOL,RADIUSMP,IMPGEN
!-----------------------------------------------------------------------
      REAL, SAVE :: R,D,W,T,TN,RM
      INTEGER, SAVE :: N,IM,IR,IG
      REAL, SAVE ::  HD,HDN,HT,XS,TA,TTA,TU,R2M1,R1M2,RMC
      INTEGER, SAVE :: NP,NC,LP
      REAL TH,SR,CTH,STH,TH1,TH2,TZ,CTZ,STZ,XL(2),XLU(2),X1,X2
      REAL XLV(2),X2M1,X1M2,T2M1,CTZV,STZV
      INTEGER MODP
!-----------------------------------------------------------------------
!   PIO4=pi/4
!-----------------------------------------------------------------------
      REAL, PARAMETER :: TWOPI=6.283185307,PI=3.141592654,TOL=1.E-6
      REAL, PARAMETER :: ZERO=0.
      IF(IPI.EQ.0) THEN
        READ (IGDFSCRB,*,IOSTAT=IER) R,D
        READ (IGDFSCRB,*,IOSTAT=IER) W,T,TN,N
        READ (IGDFSCRB,*,IOSTAT=IER) IR
        READ (IGDFSCRB,*,IOSTAT=IER) IM,RM
        READ (IGDFSCRB,*,IOSTAT=IER) IG
        IF (IER.NE.0) GOTO 99
        HT=0.5E0*T
        XS=SQRT(R*R-HT*HT)
        HD=0.5E0*D
        HDN=-HD
        TU=TN*TWOPI*0.5
        IF(T>TOL) THEN
          TA=ATAN2(HT,XS)
          NC=5
        ELSE
          TA=ZERO
          NC=2
        ENDIF
        NP=NC*N
        TTA=2.*TA
        IF(IR==0) THEN
          LP=NP+1
        ELSE
          LP=NP+2
        ENDIF
        IF(IM==0) RM=ZERO
        R2M1=(RM-R)*0.5E0
        R1M2=-R2M1
        RMC=-RM*0.5
      ELSE
        IF(IPI.LE.NP) THEN
          MODP=MOD(IPI,NC)
          IF(MODP==0) MODP=NC
          IF(MODP<5) THEN
            TZ=TU*(V+1.)+(TWOPI/N)*((IPI-1)/NC)
            CTZ=COS(TZ)
            STZ=SIN(TZ)
            CTZV=-TU*STZ
            STZV= TU*CTZ
          ELSEIF(MODP==5) THEN
            TZ=2.*TU+(TWOPI/N)*((IPI-1)/NC)
            CTZ=COS(TZ)
            STZ=SIN(TZ)
          ENDIF
        ENDIF
      ENDIF
!-----------------------------------------------------------------------
!     Bottom of the spar
!-----------------------------------------------------------------------
      IF(IPI==NP+1) THEN
         TH=PI*(U+1.)
         SR=R2M1*(V+1.)+R
         CTH=COS(TH)
         STH=SIN(TH)
         X(1)=SR*CTH
         X(2)=SR*STH
         X(3)=-D
         XU(1)=-PI*X(2)
         XU(2)= PI*X(1)
         XV(1)=R2M1*CTH
         XV(2)=R2M1*STH
      ENDIF
!-----------------------------------------------------------------------
!     Interior free surface
!-----------------------------------------------------------------------
      IF(IR==1.AND.IPI==NP+2) THEN
         TH=PI*(U+1.)
         SR=R1M2*(V+1.)+RM
         CTH=COS(TH)
         STH=SIN(TH)
         X(1)=SR*CTH
         X(2)=SR*STH
         X(3)=ZERO
         XU(1)=-PI*X(2)
         XU(2)= PI*X(1)
         XV(1)=R1M2*CTH
         XV(2)=R1M2*STH
      ENDIF
!-----------------------------------------------------------------------
!     Moonpool side wall
!-----------------------------------------------------------------------
      IF(IM/=0) THEN
        IF(IPI==LP+1) THEN
           TH=-PI*(U+1.)
           CTH=COS(TH)
           STH=SIN(TH)
           X(1)=RM*CTH
           X(2)=RM*STH
           X(3)=HDN*(V+1.)
           XU(1)= PI*X(2)
           XU(2)=-PI*X(1)
           XV(3)= HDN
        ENDIF
!-----------------------------------------------------------------------
!     Moonpool free surface
!-----------------------------------------------------------------------
        IF(IG/=0) THEN
          IF(IPI==LP+2) THEN
             TH=PI*(U+1.)
             SR=RMC*(V+1.)+RM
             CTH=COS(TH)
             STH=SIN(TH)
             X(1)=SR*CTH
             X(2)=SR*STH
             XU(1)=-PI*X(2)
             XU(2)= PI*X(1)
             XV(1)=RMC*CTH
             XV(2)=RMC*STH
          ENDIF
        ENDIF
      ENDIF
      IF(IPI.LE.NP) THEN
!-----------------------------------------------------------------------
!     Spar exterior cylindrical side between the strakes
!-----------------------------------------------------------------------
	   IF(MODP==1) THEN
           TH1=TA
           TH2=TWOPI/N-TA
           T2M1=0.5*(TH2-TH1)
           TH=T2M1*(U+1.)+TH1
           CTH=COS(TH)
           STH=SIN(TH)
           XL(1)=R*CTH
           XL(2)=R*STH
           XLU(1)=-T2M1*XL(2)
           XLU(2)= T2M1*XL(1)
!-----------------------------------------------------------------------
!     Strake surface, side 1
!-----------------------------------------------------------------------
         ELSEIF(MODP==2) THEN
           X1=R+W
           X2M1=(XS-X1)*0.5
           XL(1)=X2M1*(U+1.)+X1
           XL(2)=HT
           XLU(1)=X2M1
           XLU(2)=ZERO
!-----------------------------------------------------------------------
!     Strake surface, side 2 when thickness is not zero
!-----------------------------------------------------------------------
         ELSEIF(MODP==3) THEN
           X2=R+W
           X2M1=(X2-XS)*0.5
           XL(1)=X2M1*(U+1.)+XS
           XL(2)=-HT
           XLU(1)=X2M1
           XLU(2)=ZERO
!-----------------------------------------------------------------------
!     Strake tip surface, when thickness is not zero
!-----------------------------------------------------------------------
         ELSEIF(MODP==4) THEN
           XL(1)=R+W
           XL(2)=HT*(U+1)-HT
           XLU(1)=ZERO
           XLU(2)=HT
         ENDIF
!-----------------------------------------------------------------------
!     Rotate due to the helical form of the strakes
!     Repeated forms of the strakes is represented by additional rotation
!     by taking the rotational image of the patches on the first strake
!     the patch between the first and second strake
!-----------------------------------------------------------------------
         IF(MODP.LE.4) THEN
           X(1)=XL(1)*CTZ-XL(2)*STZ
           X(2)=XL(1)*STZ+XL(2)*CTZ
           X(3)=HDN*(V+1.)
           XU(1)=XLU(1)*CTZ-XLU(2)*STZ
           XU(2)=XLU(1)*STZ+XLU(2)*CTZ
           XV(1)=XL(1)*CTZV-XL(2)*STZV
           XV(2)=XL(1)*STZV+XL(2)*CTZV
           XV(3)=HDN
         ENDIF
!-----------------------------------------------------------------------
!     Strake bottom, when thickness is not zero
!     X=0.5*{R*COS[TA(V+1)-TA]-(R+W)}*(U+1.)+R+W
!     Y=     R*SIN[TA(V+1)-TA]
!-----------------------------------------------------------------------
         IF(MODP==5) THEN
           TH=TA*(V+1.)-TA
           CTH=COS(TH)
           STH=SIN(TH)
           X2=R+W
           X1=R*CTH
           X2M1=0.5*(X2-X1)
           XL(1)=X2M1*(U+1.)+X1
	     XL(2)=R*SIN(TH)
           XLU(1)=X2M1
           XLU(2)=ZERO
           XLV(1)=0.5*R*STH*TA*(U+1.)
           XLV(2)=R*CTH*TA
           X(1)=XL(1)*CTZ-XL(2)*STZ
           X(2)=XL(1)*STZ+XL(2)*CTZ
           X(3)=-D
           XU(1)=XLU(1)*CTZ-XLU(2)*STZ
           XU(2)=XLU(1)*STZ+XLU(2)*CTZ
           XV(1)=XLV(1)*CTZ-XLV(2)*STZ
           XV(2)=XLV(1)*STZ+XLV(2)*CTZ
         ENDIF
      ENDIF
 99   RETURN
      END

      SUBROUTINE AUV(U,V,IBI,IPI,X,XU,XV,IGDFSCRB,IER)
!-----------------------------------------------------------------------
!
!     COPYRIGHT (C) 2002      WAMIT INCORPORATED
!
!-----------------------------------------------------------------------
!
!     Version : 6.1
!
!-----------------------------------------------------------------------
!
!     Source-code file :  GEOMXACT.F
!
!-----------------------------------------------------------------------
!
!     DESCRIPTION :  (IGDEF=-13)
!
!     Subroutine defines one quadrant of an axisymmetric body with
!     hemispherical nose, cylindrical midbody, conical tail
!     Body is axisymmetric about z-axis
!     Radius, midbody length, tail length are input from GDF file
!     Origin is at midpoint of midbody
!
!     Patch 1: hemisphere
!     Patch 2: conical tail
!     Patch 3: cylinder (optional)
!
!     Options:
!        Use NPATCH=2 for no midbody (input ZMIDBODY=0.0)
!        Use NPATCH=3 for cylindrical midbody
!
!
!-----------------------------------------------------------------------
      IMPLICIT NONE
      REAL X(3),XU(3),XV(3),U,V
      INTEGER IPI,IBI,IGDFSCRB,IER
      REAL, SAVE :: RADIUS,HALFDCYL,HALFRAD,HALFTAIL
      REAL THETA,CTHETA,STHETA,PSI,CPSI,SPSI,DCYL,DTAIL,DRADDV,RSPSI,R
!-----------------------------------------------------------------------
!   PIO4=pi/4
!-----------------------------------------------------------------------
      REAL, PARAMETER :: PIO4=0.785398163397E0
      IF(IPI.EQ.0) THEN
        READ (IGDFSCRB,*,IOSTAT=IER) RADIUS, DCYL, DTAIL
        IF (IER.NE.0) GOTO 99
        HALFDCYL=0.5E0*DCYL
        HALFRAD =0.5E0*RADIUS
        HALFTAIL=0.5E0*DTAIL
        GOTO 99
      ENDIF
      THETA=PIO4*(U+1.E0)
      CTHETA=COS(THETA)
      STHETA=SIN(THETA)
!-----------------------------------------------------------------------
!  Patch 1: surface of spherical nose
!             Z=HALFDCYL:         PSI=0,  V=-1     (end)
!             Z=HALFDCYL+RADIUS:  PSI=PI/2, V=+1 (junction with midbody)
!-----------------------------------------------------------------------
      IF(IPI.EQ.1) THEN
         PSI=PIO4*(1.0+V)
         CPSI=COS(PSI)
         SPSI=SIN(PSI)
         RSPSI=RADIUS*SPSI
         X(1)=RSPSI*CTHETA
         X(2)=RSPSI*STHETA
         X(3)=HALFDCYL+RADIUS*CPSI
         XU(1)=-PIO4*X(2)
         XU(2)= PIO4*X(1)
         DRADDV=RADIUS*CPSI*PIO4
         XV(1)=DRADDV*CTHETA
         XV(2)=DRADDV*STHETA
         XV(3)=-RSPSI*PIO4
!-----------------------------------------------------------------------
!  Patch 2: conical tail  R=local radius
!-----------------------------------------------------------------------
      ELSEIF(IPI.EQ.2) THEN
         R=HALFRAD*(1.0-V)
         X(1)=R*CTHETA
         X(2)=R*STHETA
         X(3)=-HALFDCYL-HALFTAIL*(1.0+V)
         XU(1)=-PIO4*X(2)
         XU(2)= PIO4*X(1)
         XV(1)=-HALFRAD*CTHETA
         XV(2)=-HALFRAD*STHETA
         XV(3)=-HALFTAIL
!-----------------------------------------------------------------------
!  Patch 3: cylinder midbody
!-----------------------------------------------------------------------
      ELSE
         X(1)=RADIUS*CTHETA
         X(2)=RADIUS*STHETA
         X(3)=-HALFDCYL*V
         XU(1)=-PIO4*X(2)
         XU(2)= PIO4*X(1)
         XV(3)=-HALFDCYL
      ENDIF
 99   RETURN
      END

      SUBROUTINE SPAR2(U,V,IBI,IPI,X,XU,XV,IGDFSCRB,IER)
!-----------------------------------------------------------------------
!
!     COPYRIGHT (C) 2002      WAMIT INCORPORATED
!
!-----------------------------------------------------------------------
!
!     Version : 6.1
!
!-----------------------------------------------------------------------
!
!     Source-code file :  GEOMXACT.F
!
!-----------------------------------------------------------------------
!
!     DESCRIPTION :  (IGDEF=-14)
!
!     Subroutine defines the first quadrant of a spar, consisting of
!     a circular cylinder with a skirt.  Input parameters read from
!     GDF:
!
!     RAD1  radius of upper cylinder
!     RAD2  radius of skirt
!     DRAFT  draft (total to bottom)
!     SKIRT_HEIGHT height of skirt above bottom
!
!     Patch 1: upper side
!     Patch 2: top of skirt
!     Patch 3: skirt side
!     Patch 4: bottom
!     Patch 5: interior free surface (optional)
!
!     Options:
!        Use NPATCH=3 for bottom-mounted cylinder (draft=depth)
!        Use NPATCH=4 for conventional case
!        Use NPATCH=5 for irregular-frequency removal (IRR=1)
!
!-----------------------------------------------------------------------
      IMPLICIT NONE
      REAL X(3),XU(3),XV(3),U,V
      INTEGER IPI,IBI,IGDFSCRB,IER
!-----------------------------------------------------------------------
!  user-assigned definitions -- inputs from GDF must have SAVE attribute
!  HALFDRAFT, HALFRAD are evaluated/saved at first call
!-----------------------------------------------------------------------
      REAL, SAVE :: RAD1, RAD2, DRAFT, SKIRT_HEIGHT, HALFD1, HALFRAD1,  &
     &              HALFRAD2, D1, HALFR12, HALFD2
      REAL THETA,CTHETA,STHETA,SR
!-----------------------------------------------------------------------
!   PIO4=pi/4
!-----------------------------------------------------------------------
      REAL, PARAMETER :: PIO4=0.785398163397E0
!-----------------------------------------------------------------------
!   On initial call read data from GDF and evaluate parameters to SAVE
!   On all subsequent calls evaluate azimuthal angle and cos/sine
!-----------------------------------------------------------------------
      IF(IPI.EQ.0) THEN
        READ (IGDFSCRB,*,IOSTAT=IER) RAD1,RAD2
        READ (IGDFSCRB,*,IOSTAT=IER) DRAFT,SKIRT_HEIGHT
        IF (IER.NE.0) GOTO 99
        D1=DRAFT-SKIRT_HEIGHT
        HALFD1=0.5E0*D1
        HALFD2=0.5E0*SKIRT_HEIGHT
        HALFRAD1=0.5E0*RAD1
        HALFRAD2=0.5E0*RAD2
        HALFR12=0.5E0*(RAD2-RAD1)
      ELSE
         THETA=PIO4*(U+1.E0)
         CTHETA=COS(THETA)
         STHETA=SIN(THETA)
      ENDIF
!-----------------------------------------------------------------------
!  Patch 1: side (main cylinder above skirt)
!-----------------------------------------------------------------------
      IF(IPI.EQ.1) THEN
         X(1)=RAD1*CTHETA
         X(2)=RAD1*STHETA
         X(3)=-HALFD1*(V+1.E0)
         XU(1)=-PIO4*X(2)
         XU(2)= PIO4*X(1)
         XV(3)=-HALFD1
!-----------------------------------------------------------------------
!  Patch 2: top of skirt
!-----------------------------------------------------------------------
      ELSEIF(IPI.EQ.2) THEN
         SR=RAD1+HALFR12*(1.E0+V)
         X(1)=SR*CTHETA
         X(2)=SR*STHETA
         X(3)=-D1
         XU(1)=-PIO4*X(2)
         XU(2)= PIO4*X(1)
         XV(1)= HALFR12*CTHETA
         XV(2)= HALFR12*STHETA
!-----------------------------------------------------------------------
!  Patch 3: side of skirt
!-----------------------------------------------------------------------
      ELSEIF(IPI.EQ.3) THEN
         X(1)=RAD2*CTHETA
         X(2)=RAD2*STHETA
         X(3)=-D1-HALFD2*(V+1.E0)
         XU(1)=-PIO4*X(2)
         XU(2)= PIO4*X(1)
         XV(3)=-HALFD2
!-----------------------------------------------------------------------
!  Patch 4: bottom
!-----------------------------------------------------------------------
      ELSEIF(IPI.EQ.4) THEN
         SR=HALFRAD2*(1.E0-V)
         X(1)=SR*CTHETA
         X(2)=SR*STHETA
         X(3)=-DRAFT
         XU(1)=-PIO4*X(2)
         XU(2)= PIO4*X(1)
         XV(1)=-HALFRAD2*CTHETA
         XV(2)=-HALFRAD2*STHETA
!-----------------------------------------------------------------------
!  Patch 5: interior free surface
!-----------------------------------------------------------------------
      ELSEIF(IPI.EQ.5) THEN
         SR=HALFRAD1*(V+1.E0)
         X(1)=SR*CTHETA
         X(2)=SR*STHETA
         XU(1)=-PIO4*X(2)
         XU(2)= PIO4*X(1)
         XV(1)=HALFRAD1*CTHETA
         XV(2)=HALFRAD1*STHETA
      ENDIF
 99   RETURN
      END
      
      SUBROUTINE SPHERXYZ(U,V,IBI,IPI,IS,X,XU,XV,IGDFSCRB,IER)
!-----------------------------------------------------------------------
!
!     COPYRIGHT (C) 2002-8      WAMIT INCORPORATED
!
!-----------------------------------------------------------------------
!
!     Version : 6.4
!
!-----------------------------------------------------------------------
!
!     Source-code file :  GEOMXACT.F
!
!-----------------------------------------------------------------------
!
!     DESCRIPTION :  (IGDEF=-15)
!
!     Subroutine defines a submerged sphere
!     Radius input from GDF file
!     Center is located at the input point X0(1:3)
!     Error return if X0(3) >= +Radius
!     No error checks for X0(3) = - Radius
!
!     Patch 1: body surface
!
!     Symmetry planes are assumed under the following conditions:
!     (Inputs ISX,ISY must be consistent with X0(1), X0(2) in the GDF
!       file.  There are no error checks for these consistencies.)
!
!     ISX     ISY       |X0(1)|    |X0(2)|
!
!      0       0        >  0.0    >  0.0
!      1       0        =  0.0    >  0.0
!      0       1        >  0.0    =  0.0
!      1       1        =  0.0    =  0.0
!
!   MOD 6/06  sphere can intersect free surface with X0(3) < +Radius
!  11/06 HALFRAD removed
!   9/08 IS(2) added to dummy arguments, code removed to assign ISX,ISY
!        INONUMAP removed (uniform appears beter for drift forces)
!-----------------------------------------------------------------------
      IMPLICIT NONE
      REAL X(3),XU(3),XV(3),U,V
      INTEGER IPI,IBI,IGDFSCRB,IER,IS(2)
      REAL, SAVE :: RADIUS,X0(3),XUFAC,PSIM,PSIV,PSIVV,PSI0
      REAL THETA,CTHETA,STHETA,PSI,RADCPSI,RADSPSI,PRADCPSI,ONEPV
!-----------------------------------------------------------------------
!   PIO4=pi/4, PIO2=pi/2
!-----------------------------------------------------------------------
      REAL, PARAMETER :: PIO4=0.785398163397E0
      REAL, PARAMETER :: PIO2=1.5707963E0
      IF(IPI.EQ.0) THEN
         READ (IGDFSCRB,*,IOSTAT=IER) RADIUS,X0(1:3)
	   IF (IER.NE.0) GOTO 99
         IF (X0(3).GE.RADIUS) GOTO 99
         XUFAC=PIO4*(2-IS(1))*(2-IS(2))
         IF (X0(3).LE.-RADIUS) THEN
            PSI0=0.0
            PSIM=PIO2 
            PSIV=PIO2 
         ELSE
            PSI0=ACOS(-X0(3)/RADIUS)
            PSIM=PIO2+.5*PSI0
            PSIV=PIO2-.5*PSI0
         ENDIF
      ELSE
!-----------------------------------------------------------------------
!  Patch 1: body surface
!    THETA = polar angle about z-axis
!-----------------------------------------------------------------------
         IF (IS(2)==1) THEN
            THETA=XUFAC*(U+1.E0)
         ELSE
            THETA=XUFAC*U
         ENDIF
         CTHETA=COS(THETA)
         STHETA=SIN(THETA)
!-----------------------------------------------------------------------
!  PSI = azimuthal angle below +z-axis
!-----------------------------------------------------------------------
         PSI=PSIM+V*PSIV        
         RADCPSI=RADIUS*COS(PSI)
         RADSPSI=RADIUS*SIN(PSI)
         X(1)=RADSPSI*CTHETA
         X(2)=RADSPSI*STHETA
         X(3)=X0(3)+RADCPSI
         XU(1)=-XUFAC*X(2)
         XU(2)= XUFAC*X(1)
	   X(1:2)=X(1:2)+X0(1:2)
         PRADCPSI=PSIV*RADCPSI
         XV(1)=PRADCPSI*CTHETA
         XV(2)=PRADCPSI*STHETA
         XV(3)=-PSIV*RADSPSI
      ENDIF
 99   RETURN
      END

      SUBROUTINE FPSO2(U,V,IBI,IPI,X,XU,XV,IGDFSCRB,IER)
!-----------------------------------------------------------------------
!
!     COPYRIGHT (C) 2003      WAMIT INCORPORATED
!
!-----------------------------------------------------------------------
!
!     Version : 6.2
!
!-----------------------------------------------------------------------
!
!     Source-code file :  GEOMXACT.F
!
!-----------------------------------------------------------------------
!
!     DESCRIPTION :  (IGDEF=-16)
!
!     Subroutine describes one side of an FPSO hull with elliptic bow,
!     rectangular mid-body, prismatic stern.  Modified from FPSO subroutine
!     FPSO2 uses one extra patch on the bottom of the bow and separates the
!	midbody bottom into a rectangular domain.  Nonuniform mapping is an
!	option on the middle body using the same code as in the extension of
!	BARGE.
!
!   The geometric input parameters:
!     XBOW  length of bow (semi-axis of ellipse)
!     XMID  length of midbody
!     XAFT  length of prismatic afterbody
!     HBEAM half beam
!     HTRANSOM half width of transom
!     DRAFT draft
!     DTRANSOM depth of transom
!
!   General description of patches:
!
!    IPI=1 patch for the bottom of bow
!    IPI=2 patch for side of bow
!    IPI=3 patch for the bottom of midbody
!    IPI=4 patch for side of midbody
!    IPI=5 patch on transom
!    IPI=6 sloping bottom on prismatic stern
!    IPI=7 sloping side on prismatic stern
!    IPI=8,9,10 patches on free surface for IRR=1 (See MOD 1/04 below)
!
!   If NPATCH=5 the prismatic patches are omitted and the following
!   dimensions are assigned: XAFT=0.0, HTRANSOM=HBEAM, DTRANSOM=DRAFT
!
!	mod 1/04 adds patches on interior free surface if NPATCH=10:
!	IPI=8 bow, IPI=9 midbody, IPI=10 stern
!
!-------------------------------------------------------------------------
      IMPLICIT NONE
      REAL X(3),XU(3),XV(3),U,V
      INTEGER IPI,IBI,IGDFSCRB,IER
      REAL, SAVE :: XBOW,XMID,XAFT,HBEAM,HTRANSOM,DRAFT,DTRANSOM,       &
     &              XL,X1,X2,X3,X4,XSIDMID,XUSIDE,HDRAFT,HALFDT,        &
     &              HALFHT,HALFXAFT,XAFTMID,ZBAFTMID,ZBAFTFAC,          &
     &              QHT,QHB,FACUB,YBAFTMID,YSAFTFAC,ZSAFTMID,QDT,QDR,   &
     &              FACUS,THETA,YSAFTMID,HXBOW,HHBEAM,THREEQDR,QXBOW
      REAL CTHETA,STHETA,ONEMV,ONEPV,SRA,SRB,VFAC
      INTEGER, SAVE :: INONUMAP
      REAL, PARAMETER :: PI=3.1415926E0, HALF=0.5E0, ONE=1.E0
      REAL, PARAMETER :: PIO8=PI*0.125E0, PIO4=PI*.25E0
      IF (IPI.EQ.0) THEN
         READ(IGDFSCRB,*,IOSTAT=IER) XBOW,XMID,XAFT
         READ(IGDFSCRB,*,IOSTAT=IER) HBEAM,HTRANSOM
         READ(IGDFSCRB,*,IOSTAT=IER) DRAFT,DTRANSOM
         READ(IGDFSCRB,*,IOSTAT=IER) INONUMAP
         IF (IER.NE.0) GOTO 99
         XL=XBOW+XMID+XAFT
         X4=HALF*XL
         X1=-X4
         X2=X1+XAFT
         X3=X2+XMID
         XSIDMID=HALF*(X2+X3)
         XUSIDE=HALF*(X3-X2)
         HDRAFT=HALF*DRAFT
         HXBOW=HALF*XBOW
         QXBOW=0.25E0*XBOW
         HHBEAM=HALF*HBEAM
         HALFDT=HALF*DTRANSOM
         HALFHT=HALF*HTRANSOM
         HALFXAFT=HALF*XAFT
         XAFTMID=HALF*(X1+X2)
         ZBAFTMID=-HALF*(DRAFT+DTRANSOM)
         ZBAFTFAC=-HALF*(DRAFT-DTRANSOM)
         QHT=0.25*HTRANSOM
         QHB=0.25*HBEAM
         FACUB=QHB-QHT
         YBAFTMID=QHB+QHT
         YSAFTMID=HALF*(HBEAM+HTRANSOM)
         YSAFTFAC=HALF*(HBEAM-HTRANSOM)
         ZSAFTMID=0.25E0*(DTRANSOM+DRAFT)
         QDT=0.25*DTRANSOM
         QDR=0.25*DRAFT
         THREEQDR=0.75*DRAFT
         FACUS=QDR-QDT
      ELSEIF (IPI.EQ.1) THEN
!-----------------------------------------------------------------------
!     IPI=1: patch on bottom of bow
!    If INONUMAP=0 X(2) is a linear function of V
!    If INONUMAP=1 X(2) is quadratic in V with derivative=0 at V=+1
!-----------------------------------------------------------------------
         THETA=PIO4*(U+ONE)
         CTHETA=COS(THETA)
         STHETA=SIN(THETA)
         IF (INONUMAP==0) THEN
            ONEMV=1.E0-V
            SRA=ONEMV*HXBOW
            SRB=ONEMV*HHBEAM
            X(1)=X3+SRA*CTHETA
            X(2)=SRB*STHETA
            XV(1)=-HXBOW*CTHETA
            XV(2)=-HHBEAM*STHETA
         ELSE
           ONEPV=1.E0+V
           VFAC=4.E0-ONEPV*ONEPV
           SRA=QXBOW*VFAC
           SRB=QHB*VFAC
           X(1)=X3+SRA*CTHETA
           X(2)=SRB*STHETA

           XV(1)=-HXBOW*CTHETA*ONEPV
           XV(2)=-HHBEAM*STHETA*ONEPV
         ENDIF
         X(3)=-DRAFT
         XU(1)=-PIO4*SRA*STHETA
         XU(2)= PIO4*SRB*CTHETA
      ELSEIF (IPI.EQ.2) THEN
!-----------------------------------------------------------------------
!   IPI=2: patch on side of bow
!    If INONUMAP=0 X(3) is a linear function of V
!    If INONUMAP=1 X(3) is quadratic in V with derivative=0 at V=+1
!    If INONUMAP=2 X(3) is a cubic in V with derivatives=0 at V=+1,-1
!-----------------------------------------------------------------------
          THETA=PIO4*(U+ONE)
          CTHETA=COS(THETA)
          STHETA=SIN(THETA)
          X(1)=X3+XBOW*CTHETA
          X(2)=HBEAM*STHETA
          XU(1)=-PIO4*XBOW*STHETA
          XU(2)= PIO4*HBEAM*CTHETA
          IF (INONUMAP==0) THEN
            X(3)=-HDRAFT*(V+1.E0)
            XV(3)=-HDRAFT
          ELSEIF (INONUMAP==1) THEN
            ONEMV=1.E0-V
            X(3)=-QDR*(4.E0-ONEMV*ONEMV)
            XV(3)=-HDRAFT*ONEMV
          ELSE
            ONEPV=V+1.E0
            X(3)= QDR*(V-2.E0)*ONEPV*ONEPV
            XV(3)=THREEQDR*ONEPV*(V-1.E0)
          ENDIF
      ELSEIF (IPI.EQ.3) THEN
!-----------------------------------------------------------------------
!   IPI=3: patch on bottom of midbody
!    If INONUMAP=0 X(2) is a linear function of V
!    If INONUMAP=1 X(2) is quadratic in V with derivative=0 at V=+1
!-----------------------------------------------------------------------
         X(1)=XSIDMID+XUSIDE*U
         XU(1)=XUSIDE
         X(3)=-DRAFT
         IF (INONUMAP==0) THEN
            X(2)=(V+1.E0)*HHBEAM
            XV(2)=HHBEAM
         ELSE
            ONEMV=1.E0-V
            X(2)=QHB*(4.E0-ONEMV*ONEMV)
            XV(2)=HHBEAM*ONEMV
         ENDIF
      ELSEIF (IPI.EQ.4) THEN
!-----------------------------------------------------------------------
!   IPI=4: patch on side of midbody
!    If INONUMAP=0 X(3) is a linear function of V
!    If INONUMAP=1 X(3) is quadratic in V with derivative=0 at V=+1
!    If INONUMAP=2 X(3) is a cubic in V with derivatives=0 at V=+1,-1
!-----------------------------------------------------------------------
         X(1)=XSIDMID-XUSIDE*U
         X(2)=HBEAM
         XU(1)=-XUSIDE
         IF (INONUMAP==0) THEN
            X(3)=-HDRAFT*(V+1.E0)
            XV(3)=-HDRAFT
         ELSEIF (INONUMAP==1) THEN
            ONEMV=1.E0-V
            X(3)=-QDR*(4.E0-ONEMV*ONEMV)
            XV(3)=-HDRAFT*ONEMV
         ELSE
            ONEPV=V+1.E0
            X(3)= QDR*(V-2.E0)*ONEPV*ONEPV
            XV(3)=THREEQDR*ONEPV*(V-1.E0)
         ENDIF
      ELSEIF (IPI==5) THEN
!-----------------------------------------------------------------------
!    Transom
!-----------------------------------------------------------------------
         XV(3)=-HALFDT
         X(3)=(V+1.E0)*XV(3)
         XU(2)=-HALFHT
         X(1)=X1
         X(2)=(U-1.E0)*XU(2)
      ELSEIF (IPI==6) THEN
!-----------------------------------------------------------------------
!    Sloping bottom aft
!-----------------------------------------------------------------------
         X(1)=HALFXAFT*U+XAFTMID
         X(2)=(V+ONE)*(FACUB*U+YBAFTMID)
         X(3)=ZBAFTMID+U*ZBAFTFAC
         XU(1)=HALFXAFT
         XU(2)=FACUB*(V+ONE)
         XV(2)=FACUB*U+YBAFTMID
         XU(3)=ZBAFTFAC
      ELSEIF (IPI==7) THEN
!-----------------------------------------------------------------------
!    Sloping side aft
!-----------------------------------------------------------------------
         X(1)=HALFXAFT*U+XAFTMID
         X(2)=YSAFTMID+U*YSAFTFAC
         X(3)=(V-ONE)*(FACUS*U+ZSAFTMID)
         XU(1)=HALFXAFT
         XU(2)=YSAFTFAC
         XU(3)=(V-ONE)*FACUS
         XV(3)=FACUS*U+ZSAFTMID
      ELSEIF (IPI==8) THEN
!-----------------------------------------------------------------------
!    Interior free surface -- bow
!-----------------------------------------------------------------------
         THETA=PIO4*(U+ONE)
         CTHETA=COS(THETA)
         STHETA=SIN(THETA)
         IF (INONUMAP==0) THEN
            ONEPV=1.E0+V
            SRA=ONEPV*HXBOW
            SRB=ONEPV*HHBEAM
            X(1)=X3+SRA*CTHETA
            X(2)=SRB*STHETA
            XV(1)= HXBOW*CTHETA
            XV(2)= HHBEAM*STHETA
         ELSE
           ONEMV=1.E0-V
           VFAC=4.E0-ONEMV*ONEMV
           SRA=QXBOW*VFAC
           SRB=QHB*VFAC
           X(1)=X3+SRA*CTHETA
           X(2)=SRB*STHETA
           XV(1)= HXBOW*CTHETA*ONEMV
           XV(2)= HHBEAM*STHETA*ONEMV
         ENDIF
         XU(1)=-PIO4*SRA*STHETA
         XU(2)= PIO4*SRB*CTHETA
      ELSEIF (IPI==9) THEN
!-----------------------------------------------------------------------
!    Interior free surface -- midbody
!-----------------------------------------------------------------------
         X(1)=XSIDMID+XUSIDE*U
         XU(1)=XUSIDE
         IF (INONUMAP==0) THEN
            X(2)=(-V+1.E0)*HHBEAM
            XV(2)=-HHBEAM
         ELSE
            ONEPV=1.E0+V
            X(2)=QHB*(4.E0-ONEPV*ONEPV)
            XV(2)=-HHBEAM*ONEPV
         ENDIF
      ELSE
!-----------------------------------------------------------------------
!    Interior free surface -- prismatic stern
!-----------------------------------------------------------------------
         X(1)=HALFXAFT*U+XAFTMID
         X(2)=(-V+ONE)*(FACUB*U+YBAFTMID)
         XU(1)=HALFXAFT
         XU(2)=FACUB*(-V+ONE)
         XV(2)=-(FACUB*U+YBAFTMID)
      ENDIF
 99   RETURN
      END
      
      SUBROUTINE FPSO12(U,V,IBI,IPI,X,XU,XV,IGDFSCRB,IER)
!-----------------------------------------------------------------------
!
!     COPYRIGHT (C) 2003      WAMIT INCORPORATED
!
!-----------------------------------------------------------------------
!
!     Version : 6.2
!
!-----------------------------------------------------------------------
!
!     Source-code file :  GEOMXACT.F
!
!-----------------------------------------------------------------------
!
!     DESCRIPTION :  (IGDEF=-17)
!
!	derived from FPSO2 for catamaran study with two different FPSO's
!     NBODY=2 is hardwired, dimensions are input as vectors
!     Option to include interior free surface patches for IRR=1
!
!     Subroutine describes one side of an FPSO hull with elliptic bow,
!     rectangular mid-body, prismatic stern.  Modified from FPSO subroutine
!     FPSO2 uses one extra patch on the bottom of the bow and separates the
!	midbody bottom into a rectangular domain.  Nonuniform mapping is an
!	option on the middle body using the same code as in the extension of
!	BARGE.
!
!   The geometric input parameters:
!     XBOW  length of bow (semi-axis of ellipse)
!     XMID  length of midbody
!     XAFT  length of prismatic afterbody
!     HBEAM half beam
!     HTRANSOM half width of transom
!     DRAFT draft
!     DTRANSOM depth of transom
!
!   General description of patches:
!
!    IPI=1 patch for the bottom of bow
!    IPI=2 patch for side of bow
!    IPI=3 patch for the bottom of midbody
!    IPI=4 patch for side of midbody
!    IPI=5 patch on transom
!    IPI=6 sloping bottom on prismatic stern
!    IPI=7 sloping side on prismatic stern
!    IPI=8 patch for interior free surface of bow
!    IPI=9 patch for interior free surface of midbody
!    IPI=10 patch for interior free surface of stern
!
!   If NPATCH=5 the prismatic patches are omitted and the following
!   dimensions are assigned: XAFT=0.0, HTRANSOM=HBEAM, DTRANSOM=DRAFT
!   If NPATCH=7 interior free surface is omitted
!   NBODY may be modified by user in code
!
!-------------------------------------------------------------------------
      IMPLICIT NONE
      REAL X(3),XU(3),XV(3),U,V
      INTEGER IPI,IBI,IGDFSCRB,IER
!-----------------------------------------------------------------------
!    Real scalar data read or derived from each GDF file, not saved
!-----------------------------------------------------------------------
      REAL XMID,XAFT, HTRANSOM,DTRANSOM,XL,X2,X4
      REAL THETA,CTHETA,STHETA,ONEMV,ONEPV,SRA,SRB,VFAC
      INTEGER, SAVE :: INONUMAP
      INTEGER, PARAMETER :: NBODY=2
!-----------------------------------------------------------------------
!    Read data read from GDF file or derived on initial call, must be
!    saved for each body
!-----------------------------------------------------------------------
      REAL, DIMENSION(NBODY), SAVE :: XBOW, HBEAM, DRAFT,               &
     &               X1, X3, XSIDMID,XUSIDE,HDRAFT,HALFDT,              &
     &              HALFHT,HALFXAFT,XAFTMID,ZBAFTMID,ZBAFTFAC,          &
     &              QHT,QHB,FACUB,YBAFTMID,YSAFTFAC,ZSAFTMID,QDT,QDR,   &
     &              FACUS,YSAFTMID,HXBOW,HHBEAM,THREEQDR,QXBOW
      REAL, PARAMETER :: PI=3.1415926E0, HALF=0.5E0, ONE=1.E0
      REAL, PARAMETER :: PIO8=PI*0.125E0, PIO4=PI*.25E0
      IF (IPI.EQ.0) THEN
         READ(IGDFSCRB,*,IOSTAT=IER) XBOW(IBI),XMID,XAFT
         READ(IGDFSCRB,*,IOSTAT=IER) HBEAM(IBI),HTRANSOM
         READ(IGDFSCRB,*,IOSTAT=IER) DRAFT(IBI),DTRANSOM
         READ(IGDFSCRB,*,IOSTAT=IER) INONUMAP
         IF (IER.NE.0) GOTO 99
         XL=XBOW(IBI)+XMID+XAFT
         X4=HALF*XL
         X1(IBI)=-X4
         X2=X1(IBI)+XAFT
         X3(IBI)=X2+XMID
         XSIDMID(IBI)=HALF*(X2+X3(IBI))
         XUSIDE(IBI)=HALF*(X3(IBI)-X2)
         HDRAFT(IBI)=HALF*DRAFT(IBI)
         HXBOW(IBI)=HALF*XBOW(IBI)
         QXBOW(IBI)=0.25E0*XBOW(IBI)
         HHBEAM(IBI)=HALF*HBEAM(IBI)
         HALFDT(IBI)=HALF*DTRANSOM
         HALFHT(IBI)=HALF*HTRANSOM
         HALFXAFT(IBI)=HALF*XAFT
         XAFTMID(IBI)=HALF*(X1(IBI)+X2)
         ZBAFTMID(IBI)=-HALF*(DRAFT(IBI)+DTRANSOM)
         ZBAFTFAC(IBI)=-HALF*(DRAFT(IBI)-DTRANSOM)
         QHT(IBI)=0.25*HTRANSOM
         QHB(IBI)=0.25*HBEAM(IBI)
         FACUB(IBI)=QHB(IBI)-QHT(IBI)
         YBAFTMID(IBI)=QHB(IBI)+QHT(IBI)
         YSAFTMID(IBI)=HALF*(HBEAM(IBI)+HTRANSOM)
         YSAFTFAC(IBI)=HALF*(HBEAM(IBI)-HTRANSOM)
         ZSAFTMID(IBI)=0.25E0*(DTRANSOM+DRAFT(IBI))
         QDT(IBI)=0.25*DTRANSOM
         QDR(IBI)=0.25*DRAFT(IBI)
         THREEQDR(IBI)=0.75*DRAFT(IBI)
         FACUS(IBI)=QDR(IBI)-QDT(IBI)
      ELSEIF (IPI.EQ.1) THEN
!-----------------------------------------------------------------------
!     IPI=1: patch on bottom of bow
!    If INONUMAP=0 X(2) is a linear function of V
!    If INONUMAP=1 X(2) is quadratic in V with derivative=0 at V=+1
!-----------------------------------------------------------------------
         THETA=PIO4*(U+ONE)
         CTHETA=COS(THETA)
         STHETA=SIN(THETA)
         IF (INONUMAP==0) THEN
            ONEMV=1.E0-V
            SRA=ONEMV*HXBOW(IBI)
            SRB=ONEMV*HHBEAM(IBI)
            X(1)=X3(IBI)+SRA*CTHETA
            X(2)=SRB*STHETA
            XV(1)=-HXBOW(IBI)*CTHETA
            XV(2)=-HHBEAM(IBI)*STHETA
         ELSE
           ONEPV=1.E0+V
           VFAC=4.E0-ONEPV*ONEPV
           SRA=QXBOW(IBI)*VFAC
           SRB=QHB(IBI)*VFAC
           X(1)=X3(IBI)+SRA*CTHETA
           X(2)=SRB*STHETA

           XV(1)=-HXBOW(IBI)*CTHETA*ONEPV
           XV(2)=-HHBEAM(IBI)*STHETA*ONEPV
         ENDIF
         X(3)=-DRAFT(IBI)
         XU(1)=-PIO4*SRA*STHETA
         XU(2)= PIO4*SRB*CTHETA
      ELSEIF (IPI.EQ.2) THEN
!-----------------------------------------------------------------------
!   IPI=2: patch on side of bow
!    If INONUMAP=0 X(3) is a linear function of V
!    If INONUMAP=1 X(3) is quadratic in V with derivative=0 at V=+1
!    If INONUMAP=2 X(3) is a cubic in V with derivatives=0 at V=+1,-1
!-----------------------------------------------------------------------
          THETA=PIO4*(U+ONE)
          CTHETA=COS(THETA)
          STHETA=SIN(THETA)
          X(1)=X3(IBI)+XBOW(IBI)*CTHETA
          X(2)=HBEAM(IBI)*STHETA
          XU(1)=-PIO4*XBOW(IBI)*STHETA
          XU(2)= PIO4*HBEAM(IBI)*CTHETA
          IF (INONUMAP==0) THEN
            X(3)=-HDRAFT(IBI)*(V+1.E0)
            XV(3)=-HDRAFT(IBI)
          ELSEIF (INONUMAP==1) THEN
            ONEMV=1.E0-V
            X(3)=-QDR(IBI)*(4.E0-ONEMV*ONEMV)
            XV(3)=-HDRAFT(IBI)*ONEMV
          ELSE
            ONEPV=V+1.E0
            X(3)= QDR(IBI)*(V-2.E0)*ONEPV*ONEPV
            XV(3)=THREEQDR(IBI)*ONEPV*(V-1.E0)
          ENDIF
      ELSEIF (IPI.EQ.3) THEN
!-----------------------------------------------------------------------
!   IPI=3: patch on bottom of midbody
!    If INONUMAP=0 X(2) is a linear function of V
!    If INONUMAP=1 X(2) is quadratic in V with derivative=0 at V=+1
!-----------------------------------------------------------------------
         X(1)=XSIDMID(IBI)+XUSIDE(IBI)*U
         XU(1)=XUSIDE(IBI)
         X(3)=-DRAFT(IBI)
         IF (INONUMAP==0) THEN
            X(2)=(V+1.E0)*HHBEAM(IBI)
            XV(2)=HHBEAM(IBI)
         ELSE
            ONEMV=1.E0-V
            X(2)=QHB(IBI)*(4.E0-ONEMV*ONEMV)
            XV(2)=HHBEAM(IBI)*ONEMV
         ENDIF
      ELSEIF (IPI.EQ.4) THEN
!-----------------------------------------------------------------------
!   IPI=4: patch on side of midbody
!    If INONUMAP=0 X(3) is a linear function of V
!    If INONUMAP=1 X(3) is quadratic in V with derivative=0 at V=+1
!    If INONUMAP=2 X(3) is a cubic in V with derivatives=0 at V=+1,-1
!-----------------------------------------------------------------------
         X(1)=XSIDMID(IBI)-XUSIDE(IBI)*U
         X(2)=HBEAM(IBI)
         XU(1)=-XUSIDE(IBI)
         IF (INONUMAP==0) THEN
            X(3)=-HDRAFT(IBI)*(V+1.E0)
            XV(3)=-HDRAFT(IBI)
         ELSEIF (INONUMAP==1) THEN
            ONEMV=1.E0-V
            X(3)=-QDR(IBI)*(4.E0-ONEMV*ONEMV)
            XV(3)=-HDRAFT(IBI)*ONEMV
         ELSE
            ONEPV=V+1.E0
            X(3)= QDR(IBI)*(V-2.E0)*ONEPV*ONEPV
            XV(3)=THREEQDR(IBI)*ONEPV*(V-1.E0)
         ENDIF
      ELSEIF (IPI==5) THEN
!-----------------------------------------------------------------------
!    Transom
!-----------------------------------------------------------------------
         XV(3)=-HALFDT(IBI)
         X(3)=(V+1.E0)*XV(3)
         XU(2)=-HALFHT(IBI)
         X(1)=X1(IBI)
         X(2)=(U-1.E0)*XU(2)
      ELSEIF (IPI==6) THEN
!-----------------------------------------------------------------------
!    Sloping bottom aft
!-----------------------------------------------------------------------
         X(1)=HALFXAFT(IBI)*U+XAFTMID(IBI)
         X(2)=(V+ONE)*(FACUB(IBI)*U+YBAFTMID(IBI))
         X(3)=ZBAFTMID(IBI)+U*ZBAFTFAC(IBI)
         XU(1)=HALFXAFT(IBI)
         XU(2)=FACUB(IBI)*(V+ONE)
         XV(2)=FACUB(IBI)*U+YBAFTMID(IBI)
         XU(3)=ZBAFTFAC(IBI)
      ELSEIF (IPI.EQ.7) THEN
!-----------------------------------------------------------------------
!    Sloping side aft
!-----------------------------------------------------------------------
         X(1)=HALFXAFT(IBI)*U+XAFTMID(IBI)
         X(2)=YSAFTMID(IBI)+U*YSAFTFAC(IBI)
         X(3)=(V-ONE)*(FACUS(IBI)*U+ZSAFTMID(IBI))
         XU(1)=HALFXAFT(IBI)
         XU(2)=YSAFTFAC(IBI)
         XU(3)=(V-ONE)*FACUS(IBI)
         XV(3)=FACUS(IBI)*U+ZSAFTMID(IBI)
      ELSEIF (IPI.EQ.8) THEN
!-----------------------------------------------------------------------
!     IPI=8: patch on interior free surface of bow
!    If INONUMAP=0 X(2) is a linear function of V
!    If INONUMAP=1 X(2) is quadratic in V with derivative=0 at V=+1
!-----------------------------------------------------------------------
         THETA=PIO4*(-U+ONE)
         CTHETA=COS(THETA)
         STHETA=SIN(THETA)
         IF (INONUMAP==0) THEN
            ONEMV=1.E0-V
            SRA=ONEMV*HXBOW(IBI)
            SRB=ONEMV*HHBEAM(IBI)
            X(1)=X3(IBI)+SRA*CTHETA
            X(2)=SRB*STHETA
            XV(1)=-HXBOW(IBI)*CTHETA
            XV(2)=-HHBEAM(IBI)*STHETA
         ELSE
           ONEPV=1.E0+V
           VFAC=4.E0-ONEPV*ONEPV
           SRA=QXBOW(IBI)*VFAC
           SRB=QHB(IBI)*VFAC
           X(1)=X3(IBI)+SRA*CTHETA
           X(2)=SRB*STHETA

           XV(1)=-HXBOW(IBI)*CTHETA*ONEPV
           XV(2)=-HHBEAM(IBI)*STHETA*ONEPV
         ENDIF
         XU(1)= PIO4*SRA*STHETA
         XU(2)=-PIO4*SRB*CTHETA
      ELSEIF (IPI.EQ.9) THEN
!-----------------------------------------------------------------------
!   IPI=9: patch on interior free surface of midbody
!    If INONUMAP=0 X(2) is a linear function of V
!    If INONUMAP=1 X(2) is quadratic in V with derivative=0 at V=+1
!-----------------------------------------------------------------------
         X(1)=XSIDMID(IBI)-XUSIDE(IBI)*U
         XU(1)=-XUSIDE(IBI)
         IF (INONUMAP==0) THEN
            X(2)=(V+1.E0)*HHBEAM(IBI)
            XV(2)=HHBEAM(IBI)
         ELSE
            ONEMV=1.E0-V
            X(2)=QHB(IBI)*(4.E0-ONEMV*ONEMV)
            XV(2)=HHBEAM(IBI)*ONEMV
         ENDIF
      ELSEIF (IPI.EQ.10) THEN
!-----------------------------------------------------------------------
!   IPI=10: patch on interior free surface of stern
!-----------------------------------------------------------------------
         X(1)=-HALFXAFT(IBI)*U+XAFTMID(IBI)
         X(2)=-(V+ONE)*(FACUB(IBI)*U+YBAFTMID(IBI))
         XU(1)=-HALFXAFT(IBI)
         XU(2)=-FACUB(IBI)*(V+ONE)
         XV(2)=-FACUB(IBI)*U+YBAFTMID(IBI)
      ENDIF
 99   RETURN
      END
      
      SUBROUTINE TORUS_ELLIP(U,V,IBI,IPI,X,XU,XV,IGDFSCRB,IER)
!-----------------------------------------------------------------------
!
!     COPYRIGHT (C) 2000      WAMIT INCORPORATED
!
!-----------------------------------------------------------------------
!
!     Version : 6.0
!
!-----------------------------------------------------------------------
!
!     Source-code file :  GEOMXACT.F
!
!-----------------------------------------------------------------------
!
!     DESCRIPTION :  (IGDEF=-18)
!
!     Subroutine defines the torus (IPI=1) by one patch
!     Option to define the moonpool free surface by IPI=2 if NPATCH=2
!     Axis is in the free surface, DRAFT is specified separately
!        If DRAFT<>RCIRC the generating sections are ellipses
!
!     Input from GDF file:
!         RCIRC     radius of generating circles (sections of torus)
!         RAXIS     radius of axis (center of generating circles)
!         DRAFT     vertical coordinate of axis (negative if submerged)
!
!     Restrictions:
!         RAXIS > RCIRC   (open hole in center of torus)
!         ZAXIS=0 (center of ellipses is in free surface)
!-----------------------------------------------------------------------
      IMPLICIT NONE
      REAL X(3),XU(3),XV(3),U,V
      INTEGER IPI,IBI,IGDFSCRB,IER
      REAL, SAVE :: RAXIS,RCIRC,ZAXIS,PSI1,HALFRMP,DRAFT,ZFAC
      REAL THETA,CTHETA,STHETA,PSI,RCOSPSI,RSINPSI,SR,AZ,HBWL
!-----------------------------------------------------------------------
!   PIO4=pi/4
!-----------------------------------------------------------------------
      REAL, PARAMETER :: PIO4=0.78539816E0, PI=3.14159263E0
      ZAXIS=0.0E0
      IF(IPI.EQ.0) THEN
        READ (IGDFSCRB,*,IOSTAT=IER) RCIRC,RAXIS,DRAFT
        IF (IER.NE.0) GOTO 99
!-----------------------------------------------------------------------
!  First if block is for floating case, second submerged
!  HBWL = half beam of section at waterline
!  PSI1 = angle to waterline from below axis (if axis in Z=0, PSI1=pi/2,
!                                             if submerged, PSI1=pi)
!  HALFRMP = half * radius of moon pool
!-----------------------------------------------------------------------
        AZ=ABS(ZAXIS)
        IF (AZ < RCIRC) THEN
          HBWL=SQRT(RCIRC*RCIRC-AZ*AZ)
          PSI1=ATAN2(HBWL,ZAXIS)
          HALFRMP=0.5E0*(RAXIS-HBWL)
          ZFAC=DRAFT/RCIRC
        ELSE
          PSI1=PI
        ENDIF
      ELSEIF(IPI.EQ.1) THEN
!-----------------------------------------------------------------------
!  THETA = polar angle about z-axis, 0<THETA<PI/4 for the quadrant:
!     U=-1, THETA=0, Y=0 plane;  U=+1, THETA=pi/2, X=0 plane
!  PSI = angle around semi-circle below free surface, -PSI1<PSI< PSI1
!     V=-1, PSI=-PSI1  is the outer waterline (if floating)
!     V= 0, PSI=0      is the point of maximum draft
!     V=+1, PSI=+PSI1  is the inner waterline (if floating)
!-----------------------------------------------------------------------
         THETA=PIO4*(U+1.E0)
         CTHETA=COS(THETA)
         STHETA=SIN(THETA)
         PSI=PSI1*V
         RCOSPSI=RCIRC*COS(PSI)
         RSINPSI=RCIRC*SIN(PSI)
         SR=RAXIS-RSINPSI
         X(1)=SR*CTHETA
         X(2)=SR*STHETA
         X(3)=-RCOSPSI*ZFAC
         XU(1)=-PIO4*X(2)
         XU(2)= PIO4*X(1)
         XV(1)=-PSI1*RCOSPSI*CTHETA
         XV(2)=-PSI1*RCOSPSI*STHETA
         XV(3)= PSI1*RSINPSI*ZFAC
!-----------------------------------------------------------------------
!  Patch 2: moonpool free surface
!-----------------------------------------------------------------------
      ELSEIF(IPI.EQ.2) THEN
         SR=HALFRMP*(1.E0-V)
         THETA=PIO4*(U+1.E0)
         CTHETA=COS(THETA)
         STHETA=SIN(THETA)
         X(1)=SR*CTHETA
         X(2)=SR*STHETA
         XU(1)=-PIO4*X(2)
         XU(2)= PIO4*X(1)
         XV(1)=-HALFRMP*CTHETA
         XV(2)=-HALFRMP*STHETA
      ENDIF
 99   RETURN
      END

      SUBROUTINE TORUS2(U,V,IBI,IPI,X,XU,XV,IGDFSCRB,IER)
!-----------------------------------------------------------------------
!
!     COPYRIGHT (C) 2000      WAMIT INCORPORATED
!
!-----------------------------------------------------------------------
!
!     Version : 6.0
!
!-----------------------------------------------------------------------
!
!     Source-code file :  GEOMXACT.F
!
!-----------------------------------------------------------------------
!
!     DESCRIPTION :  (IGDEF=-19)
!
!     Subroutine defines two concentric toroids (IPI=1,2)
!     With elliptical sections, NPATCH=2
!
!     Input from GDF file for each torus (1,2)
!         RCIRC1,2  radius of generating circles (sections of torus)
!         RAXIS1,2  radius of axis (center of generating circles)
!         DRAFT1,2  vertical coordinate of axis (negative if submerged)
!
!     Restrictions:
!         RAXIS1 > RCIRC1   (open hole in center of torus 1)
!         RAXIS2 > RAXIS1+RCIRC1+RCIRC2 (open space between toroids)
!-----------------------------------------------------------------------
      IMPLICIT NONE
      REAL X(3),XU(3),XV(3),U,V
      INTEGER IPI,IBI,IGDFSCRB,IER
      REAL, SAVE :: RAXIS1,RCIRC1,DRAFT1,PSI1,DRAT1
      REAL, SAVE :: RAXIS2,RCIRC2,DRAFT2,PSI2,DRAT2
      REAL THETA,CTHETA,STHETA,PSI,RCOSPSI,RSINPSI,SR
!-----------------------------------------------------------------------
!   PIO4=pi/4
!-----------------------------------------------------------------------
      REAL, PARAMETER :: PIO4=0.78539816E0, PI=3.14159263E0
      IF(IPI.EQ.0) THEN
        READ (IGDFSCRB,*,IOSTAT=IER) RCIRC1,RAXIS1,DRAFT1
        READ (IGDFSCRB,*,IOSTAT=IER) RCIRC2,RAXIS2,DRAFT2
        IF (IER.NE.0) GOTO 99
        DRAT1=DRAFT1/RCIRC1
        DRAT2=DRAFT2/RCIRC2
        PSI1=PI*0.5D0
        PSI2=PI*0.5D0
      ELSEIF(IPI.EQ.1) THEN
!-----------------------------------------------------------------------
!  Patch 1: inner toroid
!  THETA = polar angle about z-axis, 0<THETA<PI/2 for the quadrant:
!     U=-1, THETA=0, Y=0 plane;  U=+1, THETA=pi/2, X=0 plane
!  PSI = angle around semi-circle below free surface, -PSI1<PSI< PSI1
!     V=-1, PSI=-PSI1  is the outer waterline (if floating)
!     V= 0, PSI=0      is the point of maximum draft
!     V=+1, PSI=+PSI1  is the inner waterline (if floating)
!-----------------------------------------------------------------------
         THETA=PIO4*(U+1.E0)
         CTHETA=COS(THETA)
         STHETA=SIN(THETA)
         PSI=PSI1*V
         RCOSPSI=RCIRC1*COS(PSI)
         RSINPSI=RCIRC1*SIN(PSI)
         SR=RAXIS1-RSINPSI
         X(1)=SR*CTHETA
         X(2)=SR*STHETA
         X(3)=-RCOSPSI*DRAT1
         XU(1)=-PIO4*X(2)
         XU(2)= PIO4*X(1)
         XV(1)=-PSI1*RCOSPSI*CTHETA
         XV(2)=-PSI1*RCOSPSI*STHETA
         XV(3)= PSI1*RSINPSI*DRAT1
!-----------------------------------------------------------------------
!  Patch 2: outer toroid
!-----------------------------------------------------------------------
      ELSEIF(IPI.EQ.2) THEN
         THETA=PIO4*(U+1.E0)
         CTHETA=COS(THETA)
         STHETA=SIN(THETA)
         PSI=PSI2*V
         RCOSPSI=RCIRC2*COS(PSI)
         RSINPSI=RCIRC2*SIN(PSI)
         SR=RAXIS2-RSINPSI
         X(1)=SR*CTHETA
         X(2)=SR*STHETA
         X(3)=-RCOSPSI*DRAT2
         XU(1)=-PIO4*X(2)
         XU(2)= PIO4*X(1)
         XV(1)=-PSI2*RCOSPSI*CTHETA
         XV(2)=-PSI2*RCOSPSI*STHETA
         XV(3)= PSI2*RSINPSI*DRAT2
      ENDIF
 99   RETURN
      END

      SUBROUTINE CIRCCYLH(U,V,IBI,IPI,X,XU,XV,IGDFSCRB,IER)
!-----------------------------------------------------------------------
!
!     COPYRIGHT (C) 2004      WAMIT INCORPORATED
!
!-----------------------------------------------------------------------
!
!     Version : 6.2
!
!-----------------------------------------------------------------------
!
!     Source-code file :  GEOMXACT.F
!
!-----------------------------------------------------------------------
!
!     DESCRIPTION :  (IGDEF=-20)
!
!     Subroutine defines the first quadrant of a horizontal circular cylinder
!     RADIUS and HALFLEN are read from GDF input file initially
!     INONUMAP             read from GDF input file initially
!
!     Patch 1: side
!     Patch 2: bottom
!     Patch 3: interior free surface (optional)
!
!     Options:
!
!       Use NPATCH=2 for conventional case
!       Use NPATCH=3 for irregular-frequency removal (IRR=1)
!
!
!-----------------------------------------------------------------------
      IMPLICIT NONE
      REAL X(3),XU(3),XV(3),U,V
      INTEGER IPI,IBI,IGDFSCRB,IER
!-----------------------------------------------------------------------
!  user-assigned definitions -- inputs from GDF must have SAVE attribute
!  HALFDRAFT, HALFRAD are evaluated/saved at first call
!-----------------------------------------------------------------------
      REAL, SAVE :: RADIUS,HALFLEN,HALFRAD,QRAD,QLEN
      REAL THETA,CTHETA,STHETA,SR,ONEPMV,XV3FAC
!-----------------------------------------------------------------------
!   PIO4=pi/4
!-----------------------------------------------------------------------
      REAL, PARAMETER :: PIO4=0.785398163397E0
!-----------------------------------------------------------------------
!   Tolerance for zero draft input
!-----------------------------------------------------------------------
      REAL, PARAMETER :: TOLD=1.0E-8
!-----------------------------------------------------------------------
!   On initial call read data from GDF and evaluate parameters to SAVE
!   On all subsequent calls evaluate azimuthal angle and cos/sine
!-----------------------------------------------------------------------
      IF(IPI.EQ.0) THEN
        READ (IGDFSCRB,*,IOSTAT=IER) RADIUS,HALFLEN
        IF (IER.NE.0) GOTO 99
        HALFRAD=0.5E0*RADIUS
        QRAD=0.25E0*RADIUS
        QLEN=0.5E0*HALFLEN
!-----------------------------------------------------------------------
!   For all patches, first evaluate the angle THETA in quadrant one:
!      U=-1 <--> THETA=0,   U=+1 <--> THETA=PI/2
!-----------------------------------------------------------------------
      ELSE
         THETA=PIO4*(U+1.E0)
         CTHETA=COS(THETA)
         STHETA=SIN(THETA)
      ENDIF
!-----------------------------------------------------------------------
!  Patch 1: side/bottom     V=-1 at midship section, V=+1 at end
!     X(1) is a linear function of V
!-----------------------------------------------------------------------
      IF (IPI.EQ.1) THEN
         X(2)=RADIUS*CTHETA
         X(3)=-RADIUS*STHETA
         XU(2)= PIO4*X(3)
         XU(3)=-PIO4*X(2)
         X(1)=QLEN*(V+1.E0)
         XV(1)=QLEN
!-----------------------------------------------------------------------
!  Patch 2: end      V=+1 on axis, -1 at corner
!    If INONUMAP=0 SR (radius) is a linear function of V
!    If INONUMAP=1 SR is quadratic in V with derivative=0 at V=-1
!-----------------------------------------------------------------------
      ELSEIF (IPI.EQ.2) THEN
         SR=HALFRAD*(1.E0-V)
         X(2)=SR*CTHETA
         X(3)=-SR*STHETA
         X(1)=HALFLEN
         XU(2)= PIO4*X(3)
         XU(3)=-PIO4*X(2)
         XV(2)=-HALFRAD*CTHETA
         XV(3)= HALFRAD*STHETA
!-----------------------------------------------------------------------
!  Patch 3: interior free surface   V=-1 on axis,  V=+1 on side
!-----------------------------------------------------------------------
      ELSEIF(IPI.EQ.3) THEN
         X(1)=(U+1.E0)*QLEN
         X(2)=(1.E0-V)*HALFRAD
         XU(1)=QLEN
         XV(2)=-HALFRAD
      ENDIF
 99   RETURN
      END

      SUBROUTINE FPSOINT(U,V,IBI,IPI,IRR,X,XU,XV,IGDFSCRB,IER)
!-----------------------------------------------------------------------
!
!     COPYRIGHT (C) 2004-8      WAMIT INCORPORATED
!
!-----------------------------------------------------------------------
!
!     Version : 6.4
!
!-----------------------------------------------------------------------
!
!     Source-code file :  GEOMXACT.F
!
!-----------------------------------------------------------------------
!
!     DESCRIPTION :  (IGDEF=-21)  ADAPTED FROM FPSO2 to include interior tanks
!
!     Subroutine describes one side of an FPSO hull with elliptic bow,
!     rectangular mid-body, prismatic stern.  Modified from FPSO subroutine
!     FPSO2 uses one extra patch on the bottom of the bow and separates the
!	midbody bottom into a rectangular domain.  Nonuniform mapping is an
!	option on the middle body using the same code as in the extension of
!	BARGE.
!	FPSOINT includes internal tank(s) defined by low-order panels
!        each tank has 4(5) patches: front, side, back, bottom, (free surface)
!        for one tank without free surface patch NPATCH=7+4=11
!     Code for tank patches is copied from bodygeom.f code used for igdef=0
!   Array XVER(4,3,NPAN) is dimensioned with NPAN=8 for two tanks
!		  Change PARAMETER statement if necessary
!
!   The geometric input parameters:
!     XBOW  length of bow (semi-axis of ellipse)
!     XMID  length of midbody
!     XAFT  length of prismatic afterbody
!     HBEAM half beam
!     HTRANSOM half width of transom
!     DRAFT draft
!     DTRANSOM depth of transom
!
!   General description of patches:
!
!    IPI=1 patch for the bottom of bow
!    IPI=2 patch for side of bow
!    IPI=3 patch for the bottom of midbody
!    IPI=4 patch for side of midbody
!    IPI=5 patch on transom
!    IPI=6 sloping bottom on prismatic stern
!    IPI=7 sloping side on prismatic stern
!
!   MOD 2/04 Inserted code for free surface patches (IRR=1), IPI=8,9,10
!   example: if IRR=1 and NTANKS=2:
!    IPI=8,9,10 interior free surface on bow, mid, aft, (same as FPSO2)
!    IPI=11,12,13,14 front,side,back,bottom of tank 1
!    IPI=15,16,17,18 front,side,back,bottom of tank 2
!   NTANKS and IRR are assigned in the GDF input file with XVER data
!      (IRR must be either 0 or 1)
!      no other changes are necessary, if each tank has four patches
!   MODS:
!     7/08 Add IRR to calling arguments, remove from READ 
!     9/08 fixed IRR for change from IRR=1 to IRR=2 or 3 
!-------------------------------------------------------------------------
      IMPLICIT NONE
      REAL X(3),XU(3),XV(3),U,V
      INTEGER IPI,IBI,IGDFSCRB,IER,IRR
      REAL, SAVE :: XBOW,XMID,XAFT,HBEAM,HTRANSOM,DRAFT,DTRANSOM,       &
     &              XL,X1,X2,X3,X4,XSIDMID,XUSIDE,HDRAFT,HALFDT,        &
     &              HALFHT,HALFXAFT,XAFTMID,ZBAFTMID,ZBAFTFAC,          &
     &              QHT,QHB,FACUB,YBAFTMID,YSAFTFAC,ZSAFTMID,QDT,QDR,   &
     &              FACUS,THETA,YSAFTMID,HXBOW,HHBEAM,THREEQDR,QXBOW
      REAL, DIMENSION(:,:,:), ALLOCATABLE, SAVE :: XVER
      REAL UP1,UM1,VP1,VM1
      REAL CTHETA,STHETA,ONEMV,ONEPV,SRA,SRB,VFAC
      INTEGER, SAVE :: INONUMAP,NPAN
      INTEGER N,I,K,NTANKS
      REAL, PARAMETER :: PI=3.1415926E0, HALF=0.5E0, ONE=1.E0
      REAL, PARAMETER :: PIO8=PI*0.125E0, PIO4=PI*.25E0
!-----------------------------------------------------------------------
!   Read data from GDF input file, including tank vertices XVER
!    NTANKS = number of internal tanks, IRR=0 or 2 must agree with
!    IRR in POT or CFG file.  Array XVER is local to this subroutine.
!    It is not the GEOMN XVER array, thus it is not saved in P2F file and
!    not deallocated at the end of the POTEN run.  This subroutine is
!    called again to initialized input data at the start of FORCE run.
!-----------------------------------------------------------------------
      IF (IPI.EQ.0) THEN
         READ(IGDFSCRB,*,IOSTAT=IER) XBOW,XMID,XAFT
         READ(IGDFSCRB,*,IOSTAT=IER) HBEAM,HTRANSOM
         READ(IGDFSCRB,*,IOSTAT=IER) DRAFT,DTRANSOM
         READ(IGDFSCRB,*,IOSTAT=IER) INONUMAP,NTANKS
         NPAN=4*NTANKS
         IF (.NOT.ALLOCATED(XVER)) ALLOCATE (XVER(4,3,NPAN))
         IF (IER.NE.0) GOTO 99
         DO N=1,NPAN
          READ (IGDFSCRB,*)                                             &
     &        ( XVER(I,1,N),XVER(I,2,N),XVER(I,3,N), I=1,4)		
         ENDDO
         XL=XBOW+XMID+XAFT
         X4=HALF*XL
         X1=-X4
         X2=X1+XAFT
         X3=X2+XMID
         XSIDMID=HALF*(X2+X3)
         XUSIDE=HALF*(X3-X2)
         HDRAFT=HALF*DRAFT
         HXBOW=HALF*XBOW
         QXBOW=0.25E0*XBOW
         HHBEAM=HALF*HBEAM
         HALFDT=HALF*DTRANSOM
         HALFHT=HALF*HTRANSOM
         HALFXAFT=HALF*XAFT
         XAFTMID=HALF*(X1+X2)
         ZBAFTMID=-HALF*(DRAFT+DTRANSOM)
         ZBAFTFAC=-HALF*(DRAFT-DTRANSOM)
         QHT=0.25*HTRANSOM
         QHB=0.25*HBEAM
         FACUB=QHB-QHT
         YBAFTMID=QHB+QHT
         YSAFTMID=HALF*(HBEAM+HTRANSOM)
         YSAFTFAC=HALF*(HBEAM-HTRANSOM)
         ZSAFTMID=0.25E0*(DTRANSOM+DRAFT)
         QDT=0.25*DTRANSOM
         QDR=0.25*DRAFT
         THREEQDR=0.75*DRAFT
         FACUS=QDR-QDT
      ELSEIF (IPI.EQ.1) THEN
!-----------------------------------------------------------------------
!     IPI=1: patch on bottom of bow
!    If INONUMAP=0 X(2) is a linear function of V
!    If INONUMAP=1 X(2) is quadratic in V with derivative=0 at V=+1
!-----------------------------------------------------------------------
         THETA=PIO4*(U+ONE)
         CTHETA=COS(THETA)
         STHETA=SIN(THETA)
         IF (INONUMAP==0) THEN
            ONEMV=1.E0-V
            SRA=ONEMV*HXBOW
            SRB=ONEMV*HHBEAM
            X(1)=X3+SRA*CTHETA
            X(2)=SRB*STHETA
            XV(1)=-HXBOW*CTHETA
            XV(2)=-HHBEAM*STHETA
         ELSE
           ONEPV=1.E0+V
           VFAC=4.E0-ONEPV*ONEPV
           SRA=QXBOW*VFAC
           SRB=QHB*VFAC
           X(1)=X3+SRA*CTHETA
           X(2)=SRB*STHETA

           XV(1)=-HXBOW*CTHETA*ONEPV
           XV(2)=-HHBEAM*STHETA*ONEPV
         ENDIF
         X(3)=-DRAFT
         XU(1)=-PIO4*SRA*STHETA
         XU(2)= PIO4*SRB*CTHETA
      ELSEIF (IPI.EQ.2) THEN
!-----------------------------------------------------------------------
!   IPI=2: patch on side of bow
!    If INONUMAP=0 X(3) is a linear function of V
!    If INONUMAP=1 X(3) is quadratic in V with derivative=0 at V=+1
!    If INONUMAP=2 X(3) is a cubic in V with derivatives=0 at V=+1,-1
!-----------------------------------------------------------------------
          THETA=PIO4*(U+ONE)
          CTHETA=COS(THETA)
          STHETA=SIN(THETA)
          X(1)=X3+XBOW*CTHETA
          X(2)=HBEAM*STHETA
          XU(1)=-PIO4*XBOW*STHETA
          XU(2)= PIO4*HBEAM*CTHETA
          IF (INONUMAP==0) THEN
            X(3)=-HDRAFT*(V+1.E0)
            XV(3)=-HDRAFT
          ELSEIF (INONUMAP==1) THEN
            ONEMV=1.E0-V
            X(3)=-QDR*(4.E0-ONEMV*ONEMV)
            XV(3)=-HDRAFT*ONEMV
          ELSE
            ONEPV=V+1.E0
            X(3)= QDR*(V-2.E0)*ONEPV*ONEPV
            XV(3)=THREEQDR*ONEPV*(V-1.E0)
          ENDIF
      ELSEIF (IPI.EQ.3) THEN
!-----------------------------------------------------------------------
!   IPI=3: patch on bottom of midbody
!    If INONUMAP=0 X(2) is a linear function of V
!    If INONUMAP=1 X(2) is quadratic in V with derivative=0 at V=+1
!-----------------------------------------------------------------------
         X(1)=XSIDMID+XUSIDE*U
         XU(1)=XUSIDE
         X(3)=-DRAFT
         IF (INONUMAP==0) THEN
            X(2)=(V+1.E0)*HHBEAM
            XV(2)=HHBEAM
         ELSE
            ONEMV=1.E0-V
            X(2)=QHB*(4.E0-ONEMV*ONEMV)
            XV(2)=HHBEAM*ONEMV
         ENDIF
      ELSEIF (IPI.EQ.4) THEN
!-----------------------------------------------------------------------
!   IPI=4: patch on side of midbody
!    If INONUMAP=0 X(3) is a linear function of V
!    If INONUMAP=1 X(3) is quadratic in V with derivative=0 at V=+1
!    If INONUMAP=2 X(3) is a cubic in V with derivatives=0 at V=+1,-1
!-----------------------------------------------------------------------
         X(1)=XSIDMID-XUSIDE*U
         X(2)=HBEAM
         XU(1)=-XUSIDE
         IF (INONUMAP==0) THEN
            X(3)=-HDRAFT*(V+1.E0)
            XV(3)=-HDRAFT
         ELSEIF (INONUMAP==1) THEN
            ONEMV=1.E0-V
            X(3)=-QDR*(4.E0-ONEMV*ONEMV)
            XV(3)=-HDRAFT*ONEMV
         ELSE
            ONEPV=V+1.E0
            X(3)= QDR*(V-2.E0)*ONEPV*ONEPV
            XV(3)=THREEQDR*ONEPV*(V-1.E0)
         ENDIF
      ELSEIF (IPI==5) THEN
!-----------------------------------------------------------------------
!    Transom
!-----------------------------------------------------------------------
         XV(3)=-HALFDT
         X(3)=(V+1.E0)*XV(3)
         XU(2)=-HALFHT
         X(1)=X1
         X(2)=(U-1.E0)*XU(2)
      ELSEIF (IPI==6) THEN
!-----------------------------------------------------------------------
!    Sloping bottom aft
!-----------------------------------------------------------------------
         X(1)=HALFXAFT*U+XAFTMID
         X(2)=(V+ONE)*(FACUB*U+YBAFTMID)
         X(3)=ZBAFTMID+U*ZBAFTFAC
         XU(1)=HALFXAFT
         XU(2)=FACUB*(V+ONE)
         XV(2)=FACUB*U+YBAFTMID
         XU(3)=ZBAFTFAC
      ELSEIF (IPI==7) THEN
!-----------------------------------------------------------------------
!    Sloping side aft
!-----------------------------------------------------------------------
         X(1)=HALFXAFT*U+XAFTMID
         X(2)=YSAFTMID+U*YSAFTFAC
         X(3)=(V-ONE)*(FACUS*U+ZSAFTMID)
         XU(1)=HALFXAFT
         XU(2)=YSAFTFAC
         XU(3)=(V-ONE)*FACUS
         XV(3)=FACUS*U+ZSAFTMID
      ELSEIF ((IRR==1).AND.(IPI==8)) THEN
!-----------------------------------------------------------------------
!    Interior free surface -- bow
!-----------------------------------------------------------------------
         THETA=PIO4*(U+ONE)
         CTHETA=COS(THETA)
         STHETA=SIN(THETA)
         IF (INONUMAP==0) THEN
            ONEPV=1.E0+V
            SRA=ONEPV*HXBOW
            SRB=ONEPV*HHBEAM
            X(1)=X3+SRA*CTHETA
            X(2)=SRB*STHETA
            XV(1)= HXBOW*CTHETA
            XV(2)= HHBEAM*STHETA
         ELSE
           ONEMV=1.E0-V
           VFAC=4.E0-ONEMV*ONEMV
           SRA=QXBOW*VFAC
           SRB=QHB*VFAC
           X(1)=X3+SRA*CTHETA
           X(2)=SRB*STHETA
           XV(1)= HXBOW*CTHETA*ONEMV
           XV(2)= HHBEAM*STHETA*ONEMV
         ENDIF
         XU(1)=-PIO4*SRA*STHETA
         XU(2)= PIO4*SRB*CTHETA
      ELSEIF ((IRR==1).AND.(IPI==9)) THEN
!-----------------------------------------------------------------------
!    Interior free surface -- midbody
!-----------------------------------------------------------------------
         X(1)=XSIDMID+XUSIDE*U
         XU(1)=XUSIDE
         IF (INONUMAP==0) THEN
            X(2)=(-V+1.E0)*HHBEAM
            XV(2)=-HHBEAM
         ELSE
            ONEPV=1.E0+V
            X(2)=QHB*(4.E0-ONEPV*ONEPV)
            XV(2)=-HHBEAM*ONEPV
         ENDIF
      ELSEIF ((IRR==1).AND.(IPI==10)) THEN
!-----------------------------------------------------------------------
!    Interior free surface -- prismatic stern
!-----------------------------------------------------------------------
         X(1)=HALFXAFT*U+XAFTMID
         X(2)=(-V+ONE)*(FACUB*U+YBAFTMID)
         XU(1)=HALFXAFT
         XU(2)=FACUB*(-V+ONE)
         XV(2)=-(FACUB*U+YBAFTMID)
      ELSE
!-----------------------------------------------------------------------
!    Interior tanks -- code copied from FLATPANL subroutine
!-----------------------------------------------------------------------
         UP1=U+1.0
         UM1=1.0-U
         VP1=V+1.0
         VM1=1.0-V
         IF(IRR==1) THEN
            K=IPI-10
         ELSE
            K=IPI-7
         ENDIF
         X(1:3)=.25*(XVER(1,1:3,K)*UP1*VM1 + XVER(2,1:3,K)*UM1*VM1      &
     &           +XVER(3,1:3,K)*UM1*VP1 + XVER(4,1:3,K)*UP1*VP1)
         XU(1:3)=.25*( XVER(1,1:3,K)*VM1 - XVER(2,1:3,K)*VM1            &
     &             -XVER(3,1:3,K)*VP1 + XVER(4,1:3,K)*VP1)
         XV(1:3)=.25*(-XVER(1,1:3,K)*UP1 - XVER(2,1:3,K)*UM1            &
     &             +XVER(3,1:3,K)*UM1 + XVER(4,1:3,K)*UP1)
      ENDIF
 99   RETURN
      END

      SUBROUTINE CIRCCYL_ARRAY(U,V,IBI,IPI,X,XU,XV,IGDFSCRB,IER)
!-----------------------------------------------------------------------
!
!     COPYRIGHT (C) 2004      WAMIT INCORPORATED
!
!-----------------------------------------------------------------------
!
!     Version : 6.2
!
!-----------------------------------------------------------------------
!
!     Source-code file :  GEOMXACT.F
!
!-----------------------------------------------------------------------
!
!     DESCRIPTION :  (IGDEF=-22)
!
!     Subroutine defines an array of circular cylinders with 2 planes of
!        symmetry (ISX=ISY=1)
!     RADIUS and DRAFT are read from GDF input file initially
!     INONUMAP             read from GDF input file initially
! if DRAFT>0: IPI=odd, side of cylinder j, IPI=2J-1
!             IPI=even, bottom of cylinder j, IPI=2j
! 
!     Patch 1: side of cylinder 1
!     Patch 2: side of cylinder 2
!     Patch N: side of cylinder N
!
!	cylinder axes are spaced at ASPACE apart, with NX*NY cylinders
!     NX must be even, NY may be odd, if odd one row of cylinders are
!         in symmetry plane Y=0 with only half of the cylinder mapped
!
!    Format of GDF file:  (sample with 2 cylinders)
!
!       cyl2W.gdf  2 cylinders at X=+/- 2.0
!       1. 9.80665  ULEN GRAV
!       1  1      ISX  ISY
!       2  -22      NPATCH  IGDEF
!       2           NLINES
!       1.0  1.0  4.0     radius, draft, aspace
!       2  1  0          NX,NY, inonumap
!
!     Options:
!       INONUMAP=0 use uniform mapping on patches
!       INONUMAP=1 use nonuniform mapping on patches (see comments below)
!
!      NB: NPATCH=number of patches in one quadrant  
!-----------------------------------------------------------------------
      IMPLICIT NONE
      REAL X(3),XU(3),XV(3),U,V
      INTEGER IPI,IBI,IGDFSCRB,IER,NX,NY,NXQ1,NYQ1,NCYL,NC,NYMOD2,MY,MX,&
     &        IPIMOD2
!-----------------------------------------------------------------------
!  user-assigned definitions -- inputs from GDF must have SAVE attribute
!  HALFDRAFT, HALFRAD are evaluated/saved at first call
!-----------------------------------------------------------------------
      REAL, SAVE :: RADIUS,DRAFT,HALFDRAFT,HALFRAD,QRAD,QDRAFT,THREEQDR,&
     &              ASPACE
      INTEGER, SAVE :: INONUMAP,NLY0
      REAL THETA,CTHETA,STHETA,SR,ONEPMV,XV3FAC,YC,DTDU
      REAL, PARAMETER :: PI  =3.14159265E0
      REAL, PARAMETER :: PIO2=PI*0.5
!-----------------------------------------------------------------------
!   Tolerance for zero draft input
!-----------------------------------------------------------------------
      REAL, PARAMETER :: TOLD=1.0E-8
      REAL, DIMENSION(:), ALLOCATABLE, SAVE :: XAXIS,YAXIS
!-----------------------------------------------------------------------
!   On initial call read data from GDF and evaluate parameters to SAVE
!   On all subsequent calls evaluate azimuthal angle and cos/sine
!-----------------------------------------------------------------------
      IF(IPI.EQ.0) THEN
        READ (IGDFSCRB,*,IOSTAT=IER) RADIUS,DRAFT,ASPACE
        READ (IGDFSCRB,*,IOSTAT=IER) NX,NY,INONUMAP
!-----------------------------------------------------------------------
! Draft <=0 is not supported
!-----------------------------------------------------------------------
        IF (DRAFT<TOLD) IER=1 
        IF (MOD(NX,2)>0) IER=1
        IF (IER.NE.0) GOTO 99
!-----------------------------------------------------------------------
!  Number of cylinders in quadrant 1
!  NYMOD2=0 if NY=even, no cylinders on x-axis, NLY0=0
!  NYMOD2=1 if NY=odd, half of cylinders on x-axis, NLY0=NXQ1
!  NLY0 is the index of the last cylinder on Y=0
!-----------------------------------------------------------------------
        NXQ1=NX/2
        NYQ1=NY/2
        NYMOD2=MOD(NY,2)        
        IF (NYMOD2==1) THEN
	    NYQ1=NYQ1+1 
          YC=0.0
          NLY0=NXQ1
        ELSE
          YC=0.5*ASPACE
          NLY0=0
        ENDIF		 
        NCYL=NXQ1*NYQ1
!-----------------------------------------------------------------------
!  Allocate arrays for coordinates of cylinder axis and assign arrays
!-----------------------------------------------------------------------
        IF (.NOT.ALLOCATED(XAXIS)) ALLOCATE (XAXIS(NCYL),YAXIS(NCYL))
        NC=0    	   
        DO MY=1,NYQ1        
          DO MX=1,NXQ1
            NC=NC+1
            XAXIS(NC)=(MX-0.5)*ASPACE              		  
            YAXIS(NC)=YC		  
          ENDDO		  
          YC=YC+ASPACE
        ENDDO				     
        HALFRAD=0.5E0*RADIUS
        QRAD=0.25E0*RADIUS
        HALFDRAFT=0.5E0*DRAFT
        QDRAFT=0.25E0*DRAFT
        THREEQDR=0.75E0*DRAFT
        GOTO 99
!-----------------------------------------------------------------------
!   For all patches, first evaluate the angle THETA in quadrant one:
!      U=-1 <--> THETA=0,   U=+1 <--> THETA=PI/2
!   NC=cylinder number
!   IPIMOD2=1: side of cylinder
!   IPIMOD2=0: bottom of cylinder
!-----------------------------------------------------------------------
      ELSE
         NC=(IPI+1)/2
         IPIMOD2=MOD(IPI,2)
         IF (NC>NLY0) THEN
           THETA=PI*U
           DTDU=PI
         ELSE
           THETA=PIO2*(U+1.0) 
           DTDU=PIO2
         ENDIF
         CTHETA=COS(THETA)
         STHETA=SIN(THETA)
      ENDIF
!-----------------------------------------------------------------------
!  side        V=-1 at free surface, V=+1 at corner
!    If INONUMAP=0 X(3) is a linear function of V
!    If INONUMAP=1 X(3) is a cubic in V with derivatives=0 at V=+1,-1
!-----------------------------------------------------------------------
      IF(IPIMOD2.EQ.1) THEN
         X(1)=RADIUS*CTHETA 
         X(2)=RADIUS*STHETA 
         XU(1)=-DTDU*X(2)
         XU(2)= DTDU*X(1)
         X(1)=X(1)+XAXIS(NC)
         X(2)=X(2)+YAXIS(NC)
         IF (INONUMAP==0) THEN
            X(3)=-HALFDRAFT*(V+1.E0)
            XV(3)=-HALFDRAFT
         ELSE
            ONEPMV=V+1.E0
            X(3)= QDRAFT*(V-2.E0)*ONEPMV*ONEPMV
            XV(3)=THREEQDR*ONEPMV*(V-1.E0)
         ENDIF
!-----------------------------------------------------------------------
!  bottom      V=+1 on axis, -1 at corner
!    If INONUMAP=0 SR (radius) is a linear function of V
!    If INONUMAP=1 SR is quadratic in V with derivative=0 at V=-1
!-----------------------------------------------------------------------
      ELSE 
         IF (INONUMAP==0) THEN
            SR=HALFRAD*(1.E0-V)
            XV3FAC=-HALFRAD
         ELSE
            ONEPMV=1.E0+V
            SR=QRAD*(4.E0-ONEPMV*ONEPMV)
            XV3FAC=-HALFRAD*ONEPMV
         ENDIF
         X(1)=SR*CTHETA
         X(2)=SR*STHETA
         X(3)=-DRAFT
         XU(1)=-DTDU*X(2)
         XU(2)= DTDU*X(1)
         X(1)=X(1)+XAXIS(NC)
         X(2)=X(2)+YAXIS(NC)
         XV(1)=XV3FAC*CTHETA
         XV(2)=XV3FAC*STHETA
      ENDIF
 99   RETURN
      END

      SUBROUTINE ELLIPINT(U,V,IBI,IPI,IRR,X,XU,XV,IGDFSCRB,IER)
!-----------------------------------------------------------------------
!
!     COPYRIGHT (C) 2004-8      WAMIT INCORPORATED
!
!-----------------------------------------------------------------------
!
!     Version : 6.4
!
!-----------------------------------------------------------------------
!
!     Source-code file :  GEOMXACT.F
!
!-----------------------------------------------------------------------
!
!     DESCRIPTION :  (IGDEF=-23)
!
!
!     Subroutine defines one side of the ellipsoid (IBI=1) by one patch
!     (ISX=0, ISY=1)
!     Center in free surface
!     Semi-axes A,B,C input from GDF file
!     With internal tanks defined by flat panels
!
!     Patch 1: body surface
!     Patch 2: interior free surface (optional)
!
!     Options:
!        Use NPATCH=1 for conventional case
!        Use NPATCH=2 for irregular-frequency removal (IRR=1)
!     If NPATCH>2 4*NTANKS patches are defined on interior tanks
!        using the code from FPSOINT subroutine
!
!   MODS:
!     7/08 Add IRR to calling arguments, remove from READ
!-----------------------------------------------------------------------
      IMPLICIT NONE
      REAL X(3),XU(3),XV(3),U,V
      INTEGER IPI,IBI,IGDFSCRB,IER,IRR
!-----------------------------------------------------------------------
!  user-assigned definitions -- inputs from GDF must have SAVE attribute
!  PI/4*(A,B,C) are evaluated once to reduce CPU time
!-----------------------------------------------------------------------
      REAL, SAVE :: A,B,C,PIO4A,PIO4B,PIO4C,HALFA,HALFB
      REAL, DIMENSION(:,:,:), ALLOCATABLE, SAVE :: XVER
      REAL THETA,CTHETA,STHETA,PSI,CPSI,SPSI,SR,SRA,SRB
      REAL UP1,UM1,VP1,VM1
      REAL, PARAMETER :: PIO4=0.785398163397E0
      REAL, PARAMETER :: PIO2=1.5707963268E0
      INTEGER, SAVE :: I,N,NPAN
      INTEGER NTANKS,K
      IF (IPI == 0) THEN
         READ (IGDFSCRB,*,IOSTAT=IER) A,B,C
         READ(IGDFSCRB,*,IOSTAT=IER) NTANKS
         NPAN=4*NTANKS
         IF (.NOT.ALLOCATED(XVER)) ALLOCATE (XVER(4,3,NPAN))
         IF (IER.NE.0) GOTO 99
		   DO N=1,NPAN
          READ (IGDFSCRB,*)                                             &
     &        ( XVER(I,1,N),XVER(I,2,N),XVER(I,3,N), I=1,4)		
	   ENDDO
        IF (IER.NE.0) GOTO 99
        PIO4A=PIO4*A
        PIO4B=PIO4*B
        PIO4C=PIO4*C
        HALFA=0.5E0*A
        HALFB=0.5E0*B
!-----------------------------------------------------------------------
!  Patch 1: body surface
!-----------------------------------------------------------------------
      ELSEIF (IPI == 1) THEN
         THETA=PIO2*(U+1.E0)
         CTHETA=COS(THETA)
         STHETA=SIN(THETA)
         PSI=PIO4*(V+3.E0)
         CPSI=COS(PSI)
         SPSI=SIN(PSI)
         X(1)=A*SPSI*CTHETA
         X(2)=B*SPSI*STHETA
         X(3)=C*CPSI
         XU(1)=-PIO4A*SPSI*STHETA*2.0
         XU(2)= PIO4B*SPSI*CTHETA*2.0
         XV(1)=PIO4A*CPSI*CTHETA
         XV(2)=PIO4B*CPSI*STHETA
         XV(3)=-PIO4C*SPSI
!-----------------------------------------------------------------------
!  Patch 2: interior free surface
!-----------------------------------------------------------------------
      ELSEIF ((IRR==1).AND.(IPI==2)) THEN
         THETA=PIO2*(U+1.E0)
         CTHETA=COS(THETA)
         STHETA=SIN(THETA)
         SR=V+1.E0
         SRA=SR*HALFA
         SRB=SR*HALFB
         X(1)=SRA*CTHETA
         X(2)=SRB*STHETA
         XU(1)=-PIO4*SRA*STHETA*2.0
         XU(2)= PIO4*SRB*CTHETA*2.0
         XV(1)=HALFA*CTHETA
         XV(2)=HALFB*STHETA
      ELSE
!-----------------------------------------------------------------------
!    Interior tanks -- code copied from FLATPANL subroutine
!-----------------------------------------------------------------------
         UP1=U+1.0
         UM1=1.0-U
         VP1=V+1.0
         VM1=1.0-V
         IF (IRR==1) THEN
	      K=IPI-2
	   ELSE
	      K=IPI-1
	   ENDIF   
         X(1:3)=.25*(XVER(1,1:3,K)*UP1*VM1 + XVER(2,1:3,K)*UM1*VM1      &
     &           +XVER(3,1:3,K)*UM1*VP1 + XVER(4,1:3,K)*UP1*VP1)
         XU(1:3)=.25*( XVER(1,1:3,K)*VM1 - XVER(2,1:3,K)*VM1            &
     &             -XVER(3,1:3,K)*VP1 + XVER(4,1:3,K)*VP1)
         XV(1:3)=.25*(-XVER(1,1:3,K)*UP1 - XVER(2,1:3,K)*UM1            &
     &             +XVER(3,1:3,K)*UM1 + XVER(4,1:3,K)*UP1)
      ENDIF
 99   RETURN
      END

      SUBROUTINE GAPLID(U,V,IBI,IPI,X,XU,XV,IGDFSCRB,IER)
!-----------------------------------------------------------------------
!
!     COPYRIGHT (C) 2004-2005      WAMIT INCORPORATED
!
!-----------------------------------------------------------------------
!
!     Version : 6.1-6.3
!
!-----------------------------------------------------------------------
!
!     Source-code file :  GEOMXACT.F
!
!-----------------------------------------------------------------------
!
!     DESCRIPTION :  (IGDEF=-24)
!
!     Subroutine defines the rectangular lid with one patch
!
!     Inputs from GDF file:
!        X1,X2 = limits of gap
!        GAP = gap width
!
!      Patch 1: rectanlge in free surface 
!
!       INONUMAP=0 use uniform mapping on patches (default)
!       INONUMAP=1,2 use nonuniform mapping on patches (see comments below)
!       4/03  option added for INONUMAP=1,2
!     
!-----------------------------------------------------------------------
      IMPLICIT NONE
      REAL X(3),XU(3),XV(3),U,V,X1,X2
      INTEGER IPI,IBI,IGDFSCRB,IER 
      REAL, SAVE :: HALFBEAM,GAP,XU0,X0
      INTEGER, SAVE :: INONUMAP
      REAL, PARAMETER :: PIO2=1.570796327
      IF(IPI.EQ.0) THEN
         READ (IGDFSCRB,*,IOSTAT=IER) X1,X2,GAP
         IF (IER.NE.0) GOTO 99
         READ (IGDFSCRB,*,END=2) INONUMAP
         GOTO 4
  2      INONUMAP=0
  4      HALFBEAM=0.5E0*GAP
         X0= 0.5*(X1+X2)
         XU0=0.5*(X2-X1)
         IF(INONUMAP.GT.3.AND.INONUMAP.LT.0) THEN
           IER=1
           GOTO 99
         ENDIF
!-----------------------------------------------------------------------
!    If INONUMAP=0 X(1),X(2) are a linear function of U
!    If INONUMAP=1 X(1) is cosine X(2) linear
!    If INONUMAP=2 X(1) linear and X(2) cosine
!    If INONUMAP=3 X(1),X(2) cosine
!-----------------------------------------------------------------------
      ELSE
         IF (INONUMAP==0.OR.INONUMAP==2) THEN
           X(1)=X0+XU0*U
           XU(1)=XU0
         ELSE
           X(1)=XU0*SIN(PIO2*U)+X0
           XU(1)=XU0*PIO2*COS(PIO2*U)
         ENDIF
         IF(INONUMAP.LE.1) THEN
           X(2)=V*HALFBEAM
           XV(2)=HALFBEAM
         ELSE
           X(2)=HALFBEAM*SIN(PIO2*V)
           XV(2)=HALFBEAM*PIO2*COS(PIO2*V)
         ENDIF
      ENDIF
 99   RETURN
      END

      SUBROUTINE CYLFIN(U,V,IBI,IPI,X,XU,XV,IGDFSCRB,IER)
!-----------------------------------------------------------------------
!
!     COPYRIGHT (C) 2004      WAMIT INCORPORATED
!
!-----------------------------------------------------------------------
!
!     Version : 6.3
!
!-----------------------------------------------------------------------
!
!     Source-code file :  GEOMXACT.F
!
!-----------------------------------------------------------------------
!
!     DESCRIPTION :  (IGDEF=-25)
!
!     Subroutine defines the first 2 quadrants of a 
!          bottom-mounted circular cylinder
!          with symmetric fins in the plane x=0
!     RADIUS and DRAFT are read from GDF input file initially
!     WIDTH             read from GDF input file initially
!
!     Patch 1: side in quadrant 1
!     Patch 2: FIN
!     Patch 3: side in quadrant 2 (optional)
!
!     If RADIUS<TOL the sides are omitted, and NPATCH=1 is input.
!      In this case the results should be the same for ISX=0 and ISX=1  
!
!    INONUMAP = 0 --> UNIFORM SPACING ON FIN, OTHERWISE COSINE MAPPING
!
!-----------------------------------------------------------------------
      IMPLICIT NONE
      REAL X(3),XU(3),XV(3),U,V,TH
      INTEGER IPI,IBI,IGDFSCRB,IER
!-----------------------------------------------------------------------
!  user-assigned definitions -- inputs from GDF must have SAVE attribute
!  HALFDRAFT, HALFRAD are evaluated/saved at first call
!-----------------------------------------------------------------------
      REAL, SAVE :: RADIUS,DRAFT,HALFDRAFT
      REAL, SAVE :: WIDTH,FINRADM,FINRADD
      INTEGER, SAVE :: INONUMAP
      REAL THETA,CTHETA,STHETA
!-----------------------------------------------------------------------
!   PIO4=pi/4
!-----------------------------------------------------------------------
      REAL, PARAMETER :: PIO4=0.785398163397E0, TOL=1.E-6
!-----------------------------------------------------------------------
!   On initial call read data from GDF and evaluate parameters to SAVE
!   On all subsequent calls evaluate azimuthal angle and cos/sine
!-----------------------------------------------------------------------
      IF(IPI.EQ.0) THEN
        READ (IGDFSCRB,*,IOSTAT=IER) RADIUS,DRAFT
        READ (IGDFSCRB,*,IOSTAT=IER) WIDTH
        INONUMAP=0
        READ (IGDFSCRB,*,END=2) INONUMAP
 2      IF (IER.NE.0) GOTO 99
        HALFDRAFT=0.5E0*DRAFT
        FINRADM=RADIUS+.5E0*WIDTH
        FINRADD=.5E0*WIDTH
!-----------------------------------------------------------------------
!   For patches on side evaluate the angle THETA in quadrant one or two
!      IPI=1: U=-1 <--> THETA=0,   U=+1 <--> THETA=PI/2
!      IPI=3: U=-1 <--> THETA=PI/2,   U=+1 <--> THETA=PI 
!   Skip both if radius <=TOL
!           V=-1 at free surface, V=+1 at bottom
!-----------------------------------------------------------------------
      ELSEIF ((RADIUS>TOL).AND.((IPI.EQ.1).OR.(IPI==3))) THEN
         THETA=PIO4*(U+IPI)
         CTHETA=COS(THETA)
         STHETA=SIN(THETA)
         X(1)=RADIUS*CTHETA
         X(2)=RADIUS*STHETA
         XU(1)=-PIO4*X(2)
         XU(2)= PIO4*X(1)
         X(3)=-HALFDRAFT*(V+1.E0)
         XV(3)=-HALFDRAFT
!-----------------------------------------------------------------------
!  Patch 2 (or 1 if radius=0): Fin  V=-1 at free surface, V=+1 at bottom
!  Cosine mapping in radial direction if INONUMAP/=0
!-----------------------------------------------------------------------
      ELSE
         X(1)=0.0
         X(3)=-HALFDRAFT*(V+1.E0)
         XV(3)=-HALFDRAFT
         IF (INONUMAP==0) THEN
            X(2)=FINRADM+U*FINRADD
            XU(2)=FINRADD
         ELSE
            TH=PIO4*(1.+U)
            X(2)=RADIUS+WIDTH*SIN(TH)
            XU(2)=PIO4*WIDTH*COS(TH)
         ENDIF	    
      ENDIF
 99   RETURN
      END
      
      SUBROUTINE CYLFIN4(U,V,IBI,IPI,IS,NPATCH,X,XU,XV,IGDFSCRB,IER)
!-----------------------------------------------------------------------
!
!     COPYRIGHT (C) 2008      WAMIT INCORPORATED
!
!-----------------------------------------------------------------------
!
!     Version : 6.4
!
!-----------------------------------------------------------------------
!
!     Source-code file :  GEOMXACT.F
!
!-----------------------------------------------------------------------
!
!     DESCRIPTION :  (IGDEF=-26)
!
!     Subroutine defines a circular cylinder with 4 symmetric fins
!     in the planes x=0 and y=0 with or without planes of symmetry
!     (used to test new code with dipole patches in planes of symmetry)
!
!     4 options:
!     ISX,ISY     NPATCH (not bottom mounted)
!      1,1          4
!      0,1          6
!      1,0          6
!      0,0          9
!
!     RADIUS and DRAFT are read from GDF input file initially
!     WIDTH             read from GDF input file initially
!
!     Patches 1:NS: sides 
!     Patches NS+1:NPATCH-1:  fins
!     Patch NPATCH:  bottom  (optional, can be bottom-mounted)
!     NS=(NPATCH-1)/2  =  1,2,4
!
!    INONUMAP = 0 --> UNIFORM SPACING ON FINS, OTHERWISE COSINE MAPPING
!
!-----------------------------------------------------------------------
      IMPLICIT NONE
      REAL X(3),XU(3),XV(3),U,V,TH
      INTEGER IPI,IBI,IGDFSCRB,IER,NPATCH,IS(2)
!-----------------------------------------------------------------------
!  user-assigned definitions -- inputs from GDF must have SAVE attribute
!  HALFDRAFT, HALFRAD are evaluated/saved at first call
!-----------------------------------------------------------------------
      REAL, SAVE :: RADIUS,DRAFT,HALFDRAFT
      REAL, SAVE :: WIDTH,FINRADM,FINRADD,THETA0,BFAC,HALFRAD
      INTEGER, SAVE :: INONUMAP,NS,IPIFIN1
      REAL THETA,CTHETA,STHETA,R,RU,SR
!-----------------------------------------------------------------------
!   PIO4=pi/4
!-----------------------------------------------------------------------
      REAL, PARAMETER :: PIO4=0.78539816E0, PIO2=1.57079633E0,          &
     &                   PI=3.1415927E0
      REAL, DIMENSION(0:2), PARAMETER :: BFC=(/ PI,PIO2,PIO4 /)
!-----------------------------------------------------------------------
!   On initial call read data from GDF and evaluate parameters to SAVE
!   On all subsequent calls evaluate azimuthal angle and cos/sine
!-----------------------------------------------------------------------
      IF(IPI.EQ.0) THEN
        READ (IGDFSCRB,*,IOSTAT=IER) RADIUS,DRAFT
        READ (IGDFSCRB,*,IOSTAT=IER) WIDTH
        INONUMAP=0
        READ (IGDFSCRB,*,END=2) INONUMAP
 2      IF (IER.NE.0) GOTO 99
        HALFRAD=0.5E0*RADIUS
        HALFDRAFT=0.5E0*DRAFT
        FINRADM=RADIUS+.5E0*WIDTH
        FINRADD=.5E0*WIDTH
!-----------------------------------------------------------------------
!   NS = number of sides to represent by patches (1,2,4)
!   THETA0 = polar angle of start of first side
!   IPIFIN1 = patch index of fin on X>0 axis
!
!   In all cases except IS=1,0 start at THETA=0 with first fin on X>0 axis
!   For IS=1,0 start at -PI/2 with first fin on Y<0 axis
!-----------------------------------------------------------------------
        NS=(NPATCH-1)/2
        IF ((IS(1)==1).AND.(IS(2)==0)) THEN
           THETA0=-PIO2
           IPIFIN1=NS+2
        ELSE
           THETA0=0.0
           IPIFIN1=NS+1
        ENDIF
        BFAC=BFC(IS(1)+IS(2))
!-----------------------------------------------------------------------
!   For patches on side evaluate the angle THETA  
!      IPI=1: U=-1 <--> THETA=THETA0,   U=+1 <--> THETA=THETA0+PI/2
!      IPI=2,3,4:   THETA= " " " + (IPI-1)*PI/2 
!   Skip both if radius <=TOL
!           V=-1 at free surface, V=+1 at bottom
!-----------------------------------------------------------------------
      ELSEIF (IPI.LE.NS) THEN
         THETA=PIO4*(U+1.0)+THETA0+(IPI-1)*PIO2
         CTHETA=COS(THETA)
         STHETA=SIN(THETA)
         X(1)=RADIUS*CTHETA
         X(2)=RADIUS*STHETA
         XU(1)=-PIO4*X(2)
         XU(2)= PIO4*X(1)
         X(3)=-HALFDRAFT*(V+1.E0)
         XV(3)=-HALFDRAFT
!-----------------------------------------------------------------------
!  Fins.   V=-1 at free surface, V=+1 at bottom
!  Cosine mapping in radial direction if INONUMAP/=0
!  First evaluate for fin on +X-axis  (IPI=NS+1)
!  Transform coordinates for other fins
!-----------------------------------------------------------------------
      ELSEIF (IPI<NPATCH) THEN
         X(3)=-HALFDRAFT*(V+1.E0)
         XV(3)=-HALFDRAFT
         IF (INONUMAP==0) THEN
            R=FINRADM+U*FINRADD
            RU=FINRADD
         ELSE
            TH=PIO4*(1.+U)
            R=RADIUS+WIDTH*SIN(TH)
            RU=PIO4*WIDTH*COS(TH)
         ENDIF
         IF (IPI==IPIFIN1) THEN
            X(1)=R
            XU(1)=RU
         ELSEIF (IPI==IPIFIN1+1) THEN
            X(2)=R
            XU(2)=RU
         ELSEIF (IPI==IPIFIN1+2) THEN
            X(1)=-R
            XU(1)=-RU
         ELSE
            X(2)=-R
            XU(2)=-RU
         ENDIF
!-----------------------------------------------------------------------
!  Patch NPATCH : Bottom  V=+1 on axis, -1 at corner
!-----------------------------------------------------------------------
      ELSEIF (IPI==NPATCH) THEN
         THETA=BFAC*(U+1.E0)+THETA0
         CTHETA=COS(THETA)
         STHETA=SIN(THETA)
         SR=HALFRAD*(1.E0-V)
         X(1)=SR*CTHETA
         X(2)=SR*STHETA
         X(3)=-DRAFT
         XU(1)=-BFAC*X(2)
         XU(2)= BFAC*X(1)
         XV(1)=-HALFRAD*CTHETA
         XV(2)=-HALFRAD*STHETA	   	    
      ENDIF
 99   RETURN
      END

      SUBROUTINE SKEW_SPHERE(U,V,IBI,IPI,X,XU,XV,IGDFSCRB,IER)
!-----------------------------------------------------------------------
!
!     COPYRIGHT (C) 2006      WAMIT INCORPORATED
!
!-----------------------------------------------------------------------
!
!     Version : 6.3
!
!-----------------------------------------------------------------------
!
!     Source-code file :  GEOMXACT.F
!
!-----------------------------------------------------------------------
!
!     DESCRIPTION :  (IGDEF=-27)
!
!     Subroutine defines TWO quadrants of a floating skewed hemisphere
!     Center in free surface
!     Radius and skew factor input from GDF file
!
!     Skew factor S is used to introduce anti-symmetric slope at 
!     waterline.  The X-coordinate of the body is defined by 
!          X(S) = X(0)+S*Z
!     Y,Z coordinates are unchanged.
!
!     Patch 1: body surface
!     Patch 2: interior free surface (optional)
!
!     Options:
!        Use NPATCH=1 for conventional case
!        Use NPATCH=2 for irregular-frequency removal (IRR=1)
!
!-----------------------------------------------------------------------
      IMPLICIT NONE
      REAL X(3),XU(3),XV(3),U,V
      INTEGER IPI,IBI,IGDFSCRB,IER
      REAL, SAVE :: RADIUS,HALFRAD,SKEW
      REAL THETA,CTHETA,STHETA,PSI,RADCPSI,RADSPSI,PRADCPSI,SR
!-----------------------------------------------------------------------
!   PIO4=pi/4
!-----------------------------------------------------------------------
      REAL, PARAMETER :: PIO4=0.785398163397E0, PIO2=1.570796327E0
      IF(IPI.EQ.0) THEN
        READ (IGDFSCRB,*,IOSTAT=IER) RADIUS,SKEW
        IF (IER.NE.0) GOTO 99
        HALFRAD=0.5E0*RADIUS
      ELSE
         THETA=PIO2*(U+1.E0)
         CTHETA=COS(THETA)
         STHETA=SIN(THETA)
      ENDIF
!-----------------------------------------------------------------------
!  Patch 1: body surface
!-----------------------------------------------------------------------
      IF(IPI.EQ.1) THEN
         PSI=PIO4*(V+3.E0)
         RADCPSI=RADIUS*COS(PSI)
         RADSPSI=RADIUS*SIN(PSI)
         X(1)=RADSPSI*CTHETA+SKEW*RADCPSI
         X(2)=RADSPSI*STHETA
         X(3)=RADCPSI
         XU(1)=-PIO2*X(2)
         XU(2)= PIO2*RADSPSI*CTHETA
         PRADCPSI=PIO4*RADCPSI
         XV(1)=PRADCPSI*CTHETA-PIO4*RADSPSI*SKEW
         XV(2)=PRADCPSI*STHETA
         XV(3)=-PIO4*RADSPSI
!-----------------------------------------------------------------------
!  Patch 2: interior free surface
!-----------------------------------------------------------------------
      ELSEIF (IPI.EQ.2) THEN
         SR=HALFRAD*(V+1.)
         X(1)=SR*CTHETA
         X(2)=SR*STHETA
         XU(1)=-PIO2*X(2)
         XU(2)= PIO2*X(1)
         XV(1)=HALFRAD*CTHETA
         XV(2)=HALFRAD*STHETA
      ENDIF
 99   RETURN
      END

      SUBROUTINE CIRCCYL_NOSYM(U,V,IBI,IPI,X,XU,XV,IGDFSCRB,IER)
!-----------------------------------------------------------------------
!
!     COPYRIGHT (C) 2005      WAMIT INCORPORATED
!
!-----------------------------------------------------------------------
!
!     Version : 6.3
!
!-----------------------------------------------------------------------
!
!     Source-code file :  GEOMXACT.F
!
!-----------------------------------------------------------------------
!
!     DESCRIPTION :  (IGDEF=-28)
!
!     Subroutine defines the entire surface of a circular cylinder
!     RADIUS and DRAFT are read from GDF input file initially
!     INONUMAP             read from GDF input file initially
!
!     Patch 1: side
!     Patch 2: bottom
!     Patch 3: interior free surface (optional)
!
!     Options:
!       INONUMAP=0 use uniform mapping on patches
!       INONUMAP=1 use nonuniform mapping on patches (see comments below)
!
!       Use NPATCH=1 for bottom-mounted cylinder (draft=depth)
!                    or for disc on free surface (draft=0.0)
!       Use NPATCH=2 for conventional case
!       Use NPATCH=3 for irregular-frequency removal (IRR=1)
!
!       6/02  Sign corrected for NPATCH=3 and INONUMAP=1
!-----------------------------------------------------------------------
      IMPLICIT NONE
      REAL X(3),XU(3),XV(3),U,V
      INTEGER IPI,IBI,IGDFSCRB,IER
!-----------------------------------------------------------------------
!  user-assigned definitions -- inputs from GDF must have SAVE attribute
!  HALFDRAFT, HALFRAD are evaluated/saved at first call
!-----------------------------------------------------------------------
      REAL, SAVE :: RADIUS,DRAFT,HALFDRAFT,HALFRAD,QRAD,QDRAFT,THREEQDR,&
     &              XSHIFT(3)
      INTEGER, SAVE :: INONUMAP
      REAL THETA,CTHETA,STHETA,SR,ONEPMV,XV3FAC
!-----------------------------------------------------------------------
!   PIO4=pi/4
!-----------------------------------------------------------------------
      REAL, PARAMETER :: PIO4=0.785398163397E0
!-----------------------------------------------------------------------
!   Tolerance for zero draft input
!-----------------------------------------------------------------------
      REAL, PARAMETER :: TOLD=1.0E-8
!-----------------------------------------------------------------------
!   On initial call read data from GDF and evaluate parameters to SAVE
!   On all subsequent calls evaluate azimuthal angle and cos/sine
!-----------------------------------------------------------------------
      IF(IPI.EQ.0) THEN
        READ (IGDFSCRB,*,IOSTAT=IER) RADIUS,DRAFT
        READ (IGDFSCRB,*,IOSTAT=IER) INONUMAP
        READ (IGDFSCRB,*,IOSTAT=IER) XSHIFT(1:3)
        IF (IER.NE.0) GOTO 99
        HALFRAD=0.5E0*RADIUS
        QRAD=0.25E0*RADIUS
        HALFDRAFT=0.5E0*DRAFT
        QDRAFT=0.25E0*DRAFT
        THREEQDR=0.75E0*DRAFT
!-----------------------------------------------------------------------
!   For all patches, first evaluate the angle THETA in quadrant one:
!      U=-1 <--> THETA=0,   U=+1 <--> THETA=PI/2
!-----------------------------------------------------------------------
      ELSE
         THETA=PIO4*(U+1.E0)*4
         CTHETA=COS(THETA)
         STHETA=SIN(THETA)
      ENDIF
!-----------------------------------------------------------------------
!  Patch 1: side        V=-1 at free surface, V=+1 at corner
!    If INONUMAP=0 X(3) is a linear function of V
!    If INONUMAP=1 X(3) is a cubic in V with derivatives=0 at V=+1,-1
!-----------------------------------------------------------------------
      IF((IPI.EQ.1).AND.(DRAFT.GE.TOLD)) THEN
         X(1)=RADIUS*CTHETA
         X(2)=RADIUS*STHETA
         XU(1)=-PIO4*X(2)*4
         XU(2)= PIO4*X(1)*4
         IF (INONUMAP==0) THEN
            X(3)=-HALFDRAFT*(V+1.E0)
            XV(3)=-HALFDRAFT
         ELSE
            ONEPMV=V+1.E0
            X(3)= QDRAFT*(V-2.E0)*ONEPMV*ONEPMV
            XV(3)=THREEQDR*ONEPMV*(V-1.E0)
         ENDIF
!-----------------------------------------------------------------------
!  Patch 2: bottom      V=+1 on axis, -1 at corner
!    If INONUMAP=0 SR (radius) is a linear function of V
!    If INONUMAP=1 SR is quadratic in V with derivative=0 at V=-1
!-----------------------------------------------------------------------
      ELSEIF((IPI.EQ.2).OR.(DRAFT<TOLD)) THEN
         IF (INONUMAP==0) THEN
            SR=HALFRAD*(1.E0-V)
            XV3FAC=-HALFRAD
         ELSE
            ONEPMV=1.E0+V
            SR=QRAD*(4.E0-ONEPMV*ONEPMV)
            XV3FAC=-HALFRAD*ONEPMV
         ENDIF
         X(1)=SR*CTHETA
         X(2)=SR*STHETA
         X(3)=-DRAFT
         XU(1)=-PIO4*X(2)*4
         XU(2)= PIO4*X(1)*4
         XV(1)=XV3FAC*CTHETA
         XV(2)=XV3FAC*STHETA
!-----------------------------------------------------------------------
!  Patch 3: interior free surface   V=-1 on axis,  V=+1 on side
!-----------------------------------------------------------------------
      ELSEIF(IPI.EQ.3) THEN
         IF (INONUMAP==0) THEN
            SR=HALFRAD*(1.E0+V)
            XV3FAC=HALFRAD
         ELSE
            ONEPMV=1.E0-V
            SR=QRAD*(4.E0-ONEPMV*ONEPMV)
            XV3FAC=HALFRAD*ONEPMV
         ENDIF
         X(1)=SR*CTHETA
         X(2)=SR*STHETA
         XU(1)=-PIO4*X(2)*4
         XU(2)= PIO4*X(1)*4
         XV(1)=XV3FAC*CTHETA
         XV(2)=XV3FAC*STHETA
      ENDIF
      X(1:3)=XSHIFT(1:3)+X(1:3)
 99   RETURN
      END

      SUBROUTINE ELLIPSOID_NOSYM_TANK(U,V,IBI,IPI,IRR,                  &
     &                                X,XU,XV,IGDFSCRB,IER)
!-----------------------------------------------------------------------
!
!     COPYRIGHT (C) 2006      WAMIT INCORPORATED
!
!-----------------------------------------------------------------------
!
!     Version : 6.3
!
!-----------------------------------------------------------------------
!
!     Source-code file :  GEOMXACT.F
!
!-----------------------------------------------------------------------
!
!     DESCRIPTION :  (IGDEF=-29)
!
!
!     Subroutine defines the ellipsoid  with tank 
!     Entire ellipsoid is described for the purpose of checking
!     the codes. The center of the ellipsoid can be moved with respect
!     to body coordinate system. The tank can be moved additionally
!     so that the relative position of the tank to the ellipsoid can
!     be changed. 
!     Semi-axes A,B,C input from GDF file
!     XS, YS and ZS are translation of the ellipsoid so that
!     it is not symmetric in body fixed coordinates system
!     XL,XB,XD,SL,SB,SD are tank dimension in entire length beam
!     and draft, SL SB SD are translation of the tank about
!     its center. Whole tank after S* translation will be translated
!     with ellipsoid. 
!
!     Patch 1: body surface
!     Patch 2: interior free surface (optional)
!     Patches 2-6 or 3-7: internal tank sides and bottom
!
!-----------------------------------------------------------------------
      IMPLICIT NONE
      REAL X(3),XU(3),XV(3),U,V
      INTEGER IPI,IBI,IRR,IGDFSCRB,IER
!-----------------------------------------------------------------------
!  user-assigned definitions -- inputs from GDF must have SAVE attribute
!  PI/4*(A,B,C) are evaluated once to reduce CPU time
!-----------------------------------------------------------------------
      REAL, SAVE :: A,B,C,PIO4A,PIO4B,PIO4C,PIA,PIB,PIC,HALFA,HALFB,    &
     &              XS,YS,ZS,XL,XB,XD,SL,SB,SD,HXL,HXB,HXD
      REAL THETA,CTHETA,STHETA,PSI,CPSI,SPSI,SR,SRA,SRB
      REAL, PARAMETER :: PIO4=0.785398163397E0,PI=3.141592654E0
      IF (IPI == 0) THEN
        READ (IGDFSCRB,*,IOSTAT=IER) A,B,C
        READ (IGDFSCRB,*,IOSTAT=IER) XS,YS,ZS
        READ (IGDFSCRB,*,IOSTAT=IER) XL,XB,XD,SL,SB,SD
        IF (IER.NE.0) GOTO 99
        PIA=PI*A
        PIB=PI*B
        PIC=PI*C
        PIO4A=PIO4*A
        PIO4B=PIO4*B
        PIO4C=PIO4*C
        HALFA=0.5E0*A
        HALFB=0.5E0*B
        HXL=0.5*XL
        HXB=0.5*XB
        HXD=0.5*XD
      ELSEIF ((IPI==1).OR.((IPI==2).AND.(IRR==1))) THEN
        THETA=PI*(U+1.E0)
        CTHETA=COS(THETA)
        STHETA=SIN(THETA)
      ENDIF
!-----------------------------------------------------------------------
!  Patch 1: body surface
!-----------------------------------------------------------------------
      IF (IPI == 1) THEN
         PSI=PIO4*(V+3.E0)
         CPSI=COS(PSI)
         SPSI=SIN(PSI)
         X(1)=A*SPSI*CTHETA
         X(2)=B*SPSI*STHETA
         X(3)=C*CPSI
         XU(1)=-PIA*SPSI*STHETA
         XU(2)= PIB*SPSI*CTHETA
         XV(1)=PIO4A*CPSI*CTHETA
         XV(2)=PIO4B*CPSI*STHETA
         XV(3)=-PIO4C*SPSI
!-----------------------------------------------------------------------
!  Patch 2: interior free surface
!-----------------------------------------------------------------------
      ELSEIF ((IPI==2).AND.(IRR==1)) THEN
         SR=V+1.E0
         SRA=SR*HALFA
         SRB=SR*HALFB
         X(1)=SRA*CTHETA
         X(2)=SRB*STHETA
         XU(1)=-PI*SRA*STHETA
         XU(2)= PI*SRB*CTHETA
         XV(1)=HALFA*CTHETA
         XV(2)=HALFB*STHETA
!-----------------------------------------------------------------------
!  Tank, port side, starboard side, aft side, forward side, bottom
!-----------------------------------------------------------------------
      ELSEIF ((IPI==2).OR.((IPI==3).AND.(IRR==1))) THEN
        X(1)=-HXL*U
        X(2)=-HXB
        X(3)=-HXD*(V+1.)
        XU(1)=-HXL
        XV(3)=-HXD
      ELSEIF ((IPI==3).OR.((IPI==4).AND.(IRR==1))) THEN
        X(1)= HXL*U
        X(2)= HXB
        X(3)=-HXD*(V+1.)
        XU(1)= HXL
        XV(3)=-HXD
      ELSEIF ((IPI==4).OR.((IPI==5).AND.(IRR==1))) THEN
        X(1)=-HXL
        X(2)= HXB*U
        X(3)=-HXD*(V+1.)
        XU(2)= HXB
        XV(3)=-HXD
      ELSEIF ((IPI==5).OR.((IPI==6).AND.(IRR==1))) THEN
        X(1)= HXL
        X(2)=-HXB*U
        X(3)=-HXD*(V+1.)
        XU(2)=-HXB
        XV(3)=-HXD
      ELSEIF ((IPI==6).OR.((IPI==7).AND.(IRR==1))) THEN
        X(1)=-HXL*U
        X(2)= HXB*V
        X(3)=-XD
        XU(1)=-HXL
        XV(2)= HXB
	ENDIF
	IF ((IPI>=3).OR.((IPI==2).AND.(IRR.NE.1))) THEN
        X(1)=X(1)+SL
        X(2)=X(2)+SB
        X(3)=X(3)+SD
      ENDIF
      X(1)=X(1)+XS
      X(2)=X(2)+YS
      X(3)=X(3)+ZS
 99   RETURN
      END

      SUBROUTINE BARGE_INT(U,V,IBI,IPI,X,XU,XV,IGDFSCRB,IER)
!-----------------------------------------------------------------------
!
!     COPYRIGHT (C) 2000-03      WAMIT INCORPORATED
!
!-----------------------------------------------------------------------
!
!     Version : 6.2
!
!-----------------------------------------------------------------------
!
!     Source-code file :  GEOMXACT.F
!
!-----------------------------------------------------------------------
!
!     DESCRIPTION :  (IGDEF=-30)
!
!     Subroutine defines the first quadrant of a rectangular tank,
!      equivalent to the interior surface of a rectangular barge
!     This subroutine is adapted from the subroutine BARGE with the
!       following modifications:
!       1) parametric coordinate U is reversed in sign, ibid for XU
!          (U and XU are unchanged on the interior free surface)
!
!     Inputs from GDF file:
!        HALFLEN = .5*Length
!        HALFBEAM= .5*Beam
!        DRAFT   = Draft
!
!      Patch 1: end
!      Patch 2: side
!      Patch 3: bottom
!
!-----------------------------------------------------------------------
      IMPLICIT NONE
      REAL X(3),XU(3),XV(3),U,V
      INTEGER IPI,IBI,IGDFSCRB,IER
      REAL, SAVE :: HALFLEN,HALFBEAM,DRAFT,                             &
     &              HALFDRAFT,QUARTLEN,QUARTBEAM 
      REAL ONEPMV,XV3FAC
      REAL, PARAMETER :: TOL=1.E-6
      IF(IPI.EQ.0) THEN
         READ (IGDFSCRB,*,IOSTAT=IER) HALFLEN,HALFBEAM,DRAFT
         IF (IER.NE.0) GOTO 99
         HALFDRAFT=0.5E0*DRAFT
         QUARTLEN=0.5E0*HALFLEN
         QUARTBEAM=0.5E0*HALFBEAM
!-----------------------------------------------------------------------
!   IPI=1: END
!-----------------------------------------------------------------------
      ELSEIF(IPI.EQ.1) THEN
         X(1)=HALFLEN
         X(2)=(-U+1.E0)*QUARTBEAM
         X(3)=-HALFDRAFT*(V+1.E0)
         XU(2)=-QUARTBEAM
         XV(3)=-HALFDRAFT
!-----------------------------------------------------------------------
!   IPI=2: SIDE
!-----------------------------------------------------------------------
      ELSEIF(IPI.EQ.2) THEN
         X(1)=(1.E0+U)*QUARTLEN
         X(2)=HALFBEAM
         XU(1)= QUARTLEN
         X(3)=-HALFDRAFT*(V+1.E0)
         XV(3)=-HALFDRAFT
!-----------------------------------------------------------------------
!   IPI=3: BOTTOM
!-----------------------------------------------------------------------
      ELSEIF (IPI.EQ.3) THEN
         X(1)=(-U+1.E0)*QUARTLEN
         X(3)=-DRAFT
         XU(1)=-QUARTLEN
         X(2)=(V+1.E0)*QUARTBEAM
         XV(2)=QUARTBEAM         
      ENDIF
 99   RETURN
      END
      
      SUBROUTINE BARGENUC(U,V,IBI,IPI,X,XU,XV,IGDFSCRB,IER)
!-----------------------------------------------------------------------
!
!     COPYRIGHT (C) 2005      WAMIT INCORPORATED
!
!-----------------------------------------------------------------------
!
!     Version : 6.2
!
!-----------------------------------------------------------------------
!
!     Source-code file :  GEOMXACT.F
!
!-----------------------------------------------------------------------
!
!     DESCRIPTION :  (IGDEF=-31)
!
!     Subroutine defines the first quadrant of a rectangular barge
!        with extra patches near the corners for nonuniform mapping
!
!     Inputs from GDF file:
!        HALFLEN = .5*Length
!        HALFBEAM= .5*Beam
!        DRAFT   = Draft
!        STRIP   = distance from the corner Over which nonunifrom
!                   mapping is used along the length 
!        STRIP < Draft, HALFBEAM and HALFLEN
!
!    INONUMAP=1:
!      Patch 1-4: end (x=HALFLEN)
!      Patch 5-8: side(y=HALFBEAM)
!      Patch 9-12: bottom (z=-DRAFT)
!      Patch 13: internal free surface (optional, use if IRR=1)
!
!     Options:
!        Use NPATCH=12 for conventional case
!        Use NPATCH=13 for irregular-frequency removal (IRR=1)
!
!-----------------------------------------------------------------------
      IMPLICIT NONE
      REAL X(3),XU(3),XV(3),U,V
      INTEGER IPI,IBI,IGDFSCRB,IER
      REAL, SAVE :: HALFLEN,HALFBEAM,DRAFT,BS,XLS,SD,BSH,XLSH,SDH,
     &              TWR,STWR,STWR3H,XLH,BH,STRIP
      INTEGER IPJ
      REAL SQRTV,SQRTU
      IF(IPI.EQ.0) THEN
         READ (IGDFSCRB,*,IOSTAT=IER) HALFLEN,HALFBEAM,DRAFT,STRIP
         IF (IER.NE.0) GOTO 99
         GOTO 4
  4      BS=HALFBEAM-STRIP
         XLS=HALFLEN-STRIP
         SD=-DRAFT+STRIP
         BSH=0.5E0*BS
         XLSH=0.5E0*XLS
         SDH=0.5*SD
         XLH=0.5E0*HALFLEN
         BH=0.5E0*HALFBEAM
         TWR=2.*SQRT(2.)
         STWR=STRIP/TWR
         STWR3H=STWR*1.5E0
!-----------------------------------------------------------------------
!   IPI=1-4 End
!-----------------------------------------------------------------------
      ELSEIF(IPI.LE.4)THEN
        X(1)=HALFLEN
        IF(IPI==1.OR.IPI==2) THEN
          X(2)=BSH*(U+1.)
          XU(2)=BSH
        ELSE
          SQRTU=SQRT(1.-U)
          X(2)=-STWR*(1.-U)*SQRTU+HALFBEAM
          XU(2)=STWR3H*SQRTU
        ENDIF
        IF(IPI==1.OR.IPI==4) THEN
          X(3)=SDH*(V+1.)
          XV(3)=SDH
        ELSE
          SQRTV=SQRT(1.-V)
          X(3)=STWR*(1.-V)*SQRTV-DRAFT
          XV(3)=-STWR3H*SQRTV
        ENDIF
!-----------------------------------------------------------------------
!   IPI=5-8 Side
!-----------------------------------------------------------------------
      ELSEIF(IPI.GE.4.AND.IPI.LE.8)THEN
        IPJ=IPI-4
        X(2)=HALFBEAM
        IF(IPJ==1.OR.IPJ==2) THEN
          X(1)=XLSH*(U+1.)
          XU(1)=XLSH
        ELSE
          SQRTU=SQRT(1.-U)
          X(1)=-STWR*(1.-U)*SQRTU+HALFLEN
          XU(1)=STWR3H*SQRTU
        ENDIF
        IF(IPJ==1.OR.IPJ==4) THEN
          X(3)=-SDH*(V-1.)
          XV(3)=-SDH
        ELSE
          SQRTV=SQRT(1.+V)
          X(3)=STWR*(1.+V)*SQRTV-DRAFT
          XV(3)=STWR3H*SQRTV
        ENDIF
!-----------------------------------------------------------------------
!   IPI=9-12 Bottom
!-----------------------------------------------------------------------
      ELSEIF(IPI.GE.9.AND.IPI.LE.12)THEN
        X(3)=-DRAFT
        IPJ=IPI-8
        IF(IPJ==1.OR.IPJ==2) THEN
          X(1)=XLSH*(U+1.)
          XU(1)=XLSH
        ELSE
          SQRTU=SQRT(1.-U)
          X(1)=-STWR*(1.-U)*SQRTU+HALFLEN
          XU(1)=STWR3H*SQRTU
        ENDIF
        IF(IPJ==1.OR.IPJ==4) THEN
          X(2)=BSH*(V+1.)
          XV(2)=BSH
        ELSE
          SQRTV=SQRT(1.-V)
          X(2)=-STWR*(1.-V)*SQRTV+HALFBEAM
          XV(2)=STWR3H*SQRTV
        ENDIF
!-----------------------------------------------------------------------
!   IPI=13 Internal free surface
!-----------------------------------------------------------------------
      ELSE
        X(1)=XLH*(1.+U)
        X(2)=BH*(1.-V)
        XU(1)=XLH
        XV(2)=-BH
      ENDIF
 99   RETURN
      END

      SUBROUTINE CCYLHSP(U,V,IBI,IPI,X,XU,XV,IGDFSCRB,IER)
!-----------------------------------------------------------------------
!
!     COPYRIGHT (C) 2006      WAMIT INCORPORATED
!
!-----------------------------------------------------------------------
!
!     Version : 6.3
!
!-----------------------------------------------------------------------
!
!     Source-code file :  GEOMXACT.F
!
!-----------------------------------------------------------------------
!
!     DESCRIPTION :  (IGDEF=-32)
!
!     Subroutine defines the first quadrant of a vessel with semi-circular 
!        sections and horizontal axis.  The vessel can be sub-divided into
!        separate segments, to permit the analysis of a hinged structure
!        using generalized modes.  The ends of the vessel are spheroidal.
!     The vessel has two planes of symmetry (isx=isy=1)
!
!     Input data from GDF file (IGDEF=-32):
!     NSEGGDF  number of separate segments (cylinders + spheroids)
!     RADIUS (radius of cylindrical sections)
!     XSEG(1:(NSEGGDF+1)/2)  array = x-coordinates of segment ends
!       
!     NSEGGDF is the total number of segments on the body
!     (NSEGGDF+1)/2 is the number of segments in x>0
!     
!     If NPATCH=(NSEGGDF+1)/2 only the submerged portion of the 
!         body is represented, with one patch for each segment
!     If NPATCH=(NSEGGDF+1)/2 + 2  the interior free surface is
!         also represented by the last two patches (the first is 
!         a rectangle covering the interior of all the cylinders, and
!         the second is a semi-ellipse covering the interior of the 
!         spheroidal end)
!
!     Options:
!
!       Use NPATCH=(NSEGGDF+1)/2 for IRR=0
!       Use NPATCH=(NSEGGDF+1)/2 + 2 for IRR=1
!
!-----------------------------------------------------------------------
      IMPLICIT NONE
      REAL X(3),XU(3),XV(3),U,V
      INTEGER IPI,IBI,IGDFSCRB,IER
!-----------------------------------------------------------------------
!  user-assigned definitions -- inputs from GDF must have SAVE attribute
!  HALFDRAFT, HALFRAD are evaluated/saved at first call
!-----------------------------------------------------------------------
      REAL, SAVE :: RADIUS,HALFRAD,QRAD,A,HALFA,PIO4A,PIO4RAD 
      REAL, DIMENSION(:), ALLOCATABLE, SAVE :: XSEG
      INTEGER, SAVE :: NSEGGDF,NXSEG,NXEGEO
      INTEGER N
      REAL THETA,CTHETA,STHETA,HLEN,PSI,CPSI,SPSI,SR,SRA,SRB
!-----------------------------------------------------------------------
!   PIO4=pi/4
!-----------------------------------------------------------------------
      REAL, PARAMETER :: PIO4=0.785398163397E0
!-----------------------------------------------------------------------
!   Tolerance for zero draft input
!-----------------------------------------------------------------------
      REAL, PARAMETER :: TOLD=1.0E-8
!-----------------------------------------------------------------------
!   On initial call read data from GDF and evaluate parameters to SAVE
!     A=semi-axis of spheroidal end
!   XSEG is an array and 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
!-----------------------------------------------------------------------
      IF(IPI.EQ.0) THEN
        READ (IGDFSCRB,*,IOSTAT=IER) NSEGGDF 
        READ (IGDFSCRB,*,IOSTAT=IER) RADIUS
        NXSEG=(NSEGGDF+1)/2 
        IF (.NOT.ALLOCATED(XSEG)) ALLOCATE (XSEG(0:NXSEG))
        XSEG(0)=0.0E0
        READ (IGDFSCRB,*,IOSTAT=IER) (XSEG(N),N=1,NXSEG)
        IF (IER.NE.0) GOTO 99
        HALFRAD=0.5E0*RADIUS
        QRAD=0.25E0*RADIUS
        A=XSEG(NXSEG)-XSEG(NXSEG-1)
        PIO4A=PIO4*A
        PIO4RAD=PIO4*RADIUS
        HALFA=0.5E0*A
        GOTO 99
      ELSEIF (IPI.NE.NXSEG+1) THEN
         THETA=PIO4*(U+1.E0)
         CTHETA=COS(THETA)
         STHETA=SIN(THETA)
      ENDIF
!-----------------------------------------------------------------------
!  Patches on cylindrical segments: 
!      U=-1 <--> THETA=0 (waterline),   U=+1 <--> THETA=PI/2 (bottom)
!      V=-1 at left end, V=+1 at right end
!-----------------------------------------------------------------------
      IF (IPI<NXSEG) THEN
         X(2)=RADIUS*CTHETA
         X(3)=-RADIUS*STHETA
         XU(2)= PIO4*X(3)
         XU(3)=-PIO4*X(2)
         HLEN=0.5E0*(XSEG(IPI)-XSEG(IPI-1))
         X(1)=HLEN*(V+1.E0)+XSEG(IPI-1)
         XV(1)=HLEN
!-----------------------------------------------------------------------
!  Patch on spheroidal end      
!      U=-1 <--> THETA=0 (right end), U=+1 <--> THETA=PI/2 (left end)
!      V=-1 <--> PSI=PI/2 (waterline),   V=+1 <--> PSI=PI (bottom)
!-----------------------------------------------------------------------
      ELSEIF (IPI==NXSEG) THEN
         PSI=PIO4*(V+3.E0)
         CPSI=COS(PSI)
         SPSI=SIN(PSI)
         X(1)=A*SPSI*CTHETA+XSEG(IPI-1)
         X(2)=RADIUS*SPSI*STHETA
         X(3)=RADIUS*CPSI
         XU(1)=-PIO4A*SPSI*STHETA
         XU(2)= PIO4RAD*SPSI*CTHETA
         XV(1)=PIO4A*CPSI*CTHETA
         XV(2)=PIO4RAD*CPSI*STHETA
         XV(3)=-PIO4RAD*SPSI
!-----------------------------------------------------------------------
!  Patch on interior free surface (midbody) V=-1 on side, V=+1 on y=0
!-----------------------------------------------------------------------
      ELSEIF(IPI==NXSEG+1) THEN
         HLEN=0.5E0*XSEG(IPI-2) 
         X(1)=(U+1.E0)*HLEN
         X(2)=(1.E0-V)*HALFRAD
         XU(1)=HLEN
         XV(2)=-HALFRAD
!-----------------------------------------------------------------------
!  Patch on interior free surface (spheroidal end) 
!      U=-1 <--> THETA=0 (right end), U=+1 <--> THETA=PI/2 (left end)
!      V=-1 <--> PSI=PI/2 (bottom),   V=+1 <--> PSI=PI (waterline)
!       V=-1 on axis, V=+1 on waterline
!-----------------------------------------------------------------------
      ELSEIF(IPI==NXSEG+2) THEN
         SR=V+1.E0
         SRA=SR*HALFA
         SRB=SR*HALFRAD
         X(1)=SRA*CTHETA+XSEG(NXSEG-1)
         X(2)=SRB*STHETA
         XU(1)=-PIO4*SRA*STHETA
         XU(2)= PIO4*SRB*CTHETA
         XV(1)=HALFA*CTHETA
         XV(2)=HALFRAD*STHETA
      ENDIF
 99   RETURN
      END      

      SUBROUTINE CCYL_CS(U,V,IBI,IPI,X,XU,XV,IGDFSCRB,IER)
!-----------------------------------------------------------------------
!
!     COPYRIGHT (C) 2005-08      WAMIT INCORPORATED
!
!-----------------------------------------------------------------------
!
!     Version : 6.4
!
!-----------------------------------------------------------------------
!
!     Source-code file :  GEOMXACT.F
!
!-----------------------------------------------------------------------
!
!     DESCRIPTION :  (IGDEF=-1001)
!
!     Subroutine defines the first quadrant of a control surface
!     of circular cylinder
!     RADIUS, Radius of inclosed cylinder and DRAFT are read 
!             GDF input file initially
!     INONUMAP             read from GDF input file initially
!
!     Patch 1: free surface
!     Patch 2: side
!     Patch 3: bottom 
!
!     Options:
!       INONUMAP=0 use uniform mapping on patches
!       INONUMAP=1 use nonuniform mapping on patches (see comments below)
! 
!       2/08  order of patches changed as above
!       
!-----------------------------------------------------------------------
      IMPLICIT NONE
      REAL X(3),XU(3),XV(3),U,V
      INTEGER IPI,IBI,IGDFSCRB,IER
!-----------------------------------------------------------------------
!  user-assigned definitions -- inputs from GDF must have SAVE attribute
!  HALFDRAFT, HALFRAD are evaluated/saved at first call
!-----------------------------------------------------------------------
      REAL, SAVE :: RADIUS,DRAFT,HALFDRAFT,HALFRAD,QRAD,QDRAFT,THREEQDR
      REAL, SAVE :: RADIUSI,HALFTHICK
      INTEGER, SAVE :: INONUMAP
      REAL THETA,CTHETA,STHETA,SR,ONEPMV,XV3FAC
!-----------------------------------------------------------------------
!   PIO4=pi/4
!-----------------------------------------------------------------------
      REAL, PARAMETER :: PIO4=0.785398163397E0
!-----------------------------------------------------------------------
!   On initial call read data from GDF and evaluate parameters to SAVE
!   On all subsequent calls evaluate azimuthal angle and cos/sine
!-----------------------------------------------------------------------
      IF(IPI.EQ.0) THEN
        READ (IGDFSCRB,*,IOSTAT=IER) RADIUS,DRAFT,RADIUSI
        READ (IGDFSCRB,*,IOSTAT=IER) INONUMAP
        IF (IER.NE.0) GOTO 99
        HALFRAD=0.5E0*RADIUS
        QRAD=0.25E0*RADIUS
        HALFDRAFT=0.5E0*DRAFT
        QDRAFT=0.25E0*DRAFT
        THREEQDR=0.75E0*DRAFT
        HALFTHICK=0.5*(RADIUS-RADIUSI)
        IF(HALFTHICK.LE.0.) THEN
          IER=1
          GOTO 99
        ENDIF
!-----------------------------------------------------------------------
!   For all patches, first evaluate the angle THETA in quadrant one:
!      U=-1 <--> THETA=0,   U=+1 <--> THETA=PI/2
!-----------------------------------------------------------------------
      ELSE
         THETA=PIO4*(U+1.E0)
         CTHETA=COS(THETA)
         STHETA=SIN(THETA)
      ENDIF
!-----------------------------------------------------------------------
!  Patch 1: free surface   V=-1 RadiusI and V=+1 on side=RASIUS
!         r=(r2-r1)/2*(v+1)+r1 r=(r2-r1)/4*(v+1)**2+r1
!-----------------------------------------------------------------------
      IF(IPI.EQ.1) THEN
         IF (INONUMAP==0) THEN
            SR=HALFTHICK*(1.E0+V)+RADIUSI
            XV3FAC=HALFTHICK
         ELSE
            ONEPMV=1.E0+V
            SR=0.5*HALFTHICK*ONEPMV*ONEPMV+RADIUSI
            XV3FAC=HALFTHICK*ONEPMV
         ENDIF
         X(1)=SR*CTHETA
         X(2)=SR*STHETA
         XU(1)=-PIO4*X(2)
         XU(2)= PIO4*X(1)
         XV(1)=XV3FAC*CTHETA
         XV(2)=XV3FAC*STHETA
!-----------------------------------------------------------------------
!  Patch 2: side        V=-1 at free surface, V=+1 at corner
!    If INONUMAP=0 X(3) is a linear function of V
!    If INONUMAP=1 X(3) is a cubic in V with derivatives=0 at V=+1,-1
!-----------------------------------------------------------------------
      ELSEIF (IPI.EQ.2) THEN
         X(1)=RADIUS*CTHETA
         X(2)=RADIUS*STHETA
         XU(1)=-PIO4*X(2)
         XU(2)= PIO4*X(1)
         IF (INONUMAP==0) THEN
            X(3)=-HALFDRAFT*(V+1.E0)
            XV(3)=-HALFDRAFT
         ELSE
            ONEPMV=V+1.E0
            X(3)= QDRAFT*(V-2.E0)*ONEPMV*ONEPMV
            XV(3)=THREEQDR*ONEPMV*(V-1.E0)
         ENDIF
!-----------------------------------------------------------------------
!  Patch 3: bottom      V=+1 on axis, -1 at corner
!    If INONUMAP=0 SR (radius) is a linear function of V
!    If INONUMAP=1 SR is quadratic in V with derivative=0 at V=-1
!-----------------------------------------------------------------------
      ELSE
         IF (INONUMAP==0) THEN
            SR=HALFRAD*(1.E0-V)
            XV3FAC=-HALFRAD
         ELSE
            ONEPMV=1.E0+V
            SR=QRAD*(4.E0-ONEPMV*ONEPMV)
            XV3FAC=-HALFRAD*ONEPMV
         ENDIF
         X(1)=SR*CTHETA
         X(2)=SR*STHETA
         X(3)=-DRAFT
         XU(1)=-PIO4*X(2)
         XU(2)= PIO4*X(1)
         XV(1)=XV3FAC*CTHETA
         XV(2)=XV3FAC*STHETA
      ENDIF
 99   RETURN
      END

      SUBROUTINE CCYL_CS_NOSYSM(U,V,IBI,IPI,X,XU,XV,IGDFSCRB,IER)
!-----------------------------------------------------------------------
!
!     COPYRIGHT (C) 2005      WAMIT INCORPORATED
!
!-----------------------------------------------------------------------
!
!     Version : 6.3
!
!-----------------------------------------------------------------------
!
!     Source-code file :  GEOMXACT.F
!
!-----------------------------------------------------------------------
!
!     DESCRIPTION :  (IGDEF=-1002)
!
!     Subroutine defines entire control surface
!     of circular cylinder
!     RADIUS, Radius of inclosed cylinder and DRAFT are read 
!             GDF input file initially
!     INONUMAP             read from GDF input file initially
!
!     Patch 1: side
!     Patch 2: free surface
!     Patch 3: bottom 
!
!     Options:
!       INONUMAP=0 use uniform mapping on patches
!       INONUMAP=1 use nonuniform mapping on patches (see comments below)
!
!
!       6/02  Sign corrected for NPATCH=3 and INONUMAP=1
!-----------------------------------------------------------------------
      IMPLICIT NONE
      REAL X(3),XU(3),XV(3),U,V
      INTEGER IPI,IBI,IGDFSCRB,IER
!-----------------------------------------------------------------------
!  user-assigned definitions -- inputs from GDF must have SAVE attribute
!  HALFDRAFT, HALFRAD are evaluated/saved at first call
!-----------------------------------------------------------------------
      REAL, SAVE :: RADIUS,DRAFT,HALFDRAFT,HALFRAD,QRAD,QDRAFT,THREEQDR
      REAL, SAVE :: RADIUSI,HALFTHICK,XSHIFT(3)
      INTEGER, SAVE :: INONUMAP
      REAL THETA,CTHETA,STHETA,SR,ONEPMV,XV3FAC
!-----------------------------------------------------------------------
!   PIO4=pi/4
!-----------------------------------------------------------------------
      REAL, PARAMETER :: PI=3.14159263E0
!-----------------------------------------------------------------------
!   Tolerance for zero draft input
!-----------------------------------------------------------------------
      REAL, PARAMETER :: TOLD=1.0E-8
!-----------------------------------------------------------------------
!   On initial call read data from GDF and evaluate parameters to SAVE
!   On all subsequent calls evaluate azimuthal angle and cos/sine
!-----------------------------------------------------------------------
      IF(IPI.EQ.0) THEN
        READ (IGDFSCRB,*,IOSTAT=IER) RADIUS,DRAFT,RADIUSI
        READ (IGDFSCRB,*,IOSTAT=IER) INONUMAP
        READ (IGDFSCRB,*,IOSTAT=IER) XSHIFT(1:3)
        IF (IER.NE.0) GOTO 99
        HALFRAD=0.5E0*RADIUS
        QRAD=0.25E0*RADIUS
        HALFDRAFT=0.5E0*DRAFT
        QDRAFT=0.25E0*DRAFT
        THREEQDR=0.75E0*DRAFT
        HALFTHICK=0.5*(RADIUS-RADIUSI)
        IF(HALFTHICK.LE.0.) THEN
          IER=1
          GOTO 99
        ENDIF
!-----------------------------------------------------------------------
!   For all patches, first evaluate the angle THETA 
!      U=-1 <--> THETA=0,   U=+1 <--> THETA=2*PI
!-----------------------------------------------------------------------
      ELSE
         THETA=PI*(U+1.E0)
         CTHETA=COS(THETA)
         STHETA=SIN(THETA)
      ENDIF
!-----------------------------------------------------------------------
!  Patch 1: side        V=-1 at free surface, V=+1 at corner
!    If INONUMAP=0 X(3) is a linear function of V
!    If INONUMAP=1 X(3) is a cubic in V with derivatives=0 at V=+1,-1
!-----------------------------------------------------------------------
      IF(IPI.EQ.1) THEN
         X(1)=RADIUS*CTHETA
         X(2)=RADIUS*STHETA
         XU(1)=-PI*X(2)
         XU(2)= PI*X(1)
         IF (INONUMAP==0) THEN
            X(3)=-HALFDRAFT*(V+1.E0)
            XV(3)=-HALFDRAFT
         ELSE
            ONEPMV=V+1.E0
            X(3)= QDRAFT*(V-2.E0)*ONEPMV*ONEPMV
            XV(3)=THREEQDR*ONEPMV*(V-1.E0)
         ENDIF
!-----------------------------------------------------------------------
!  Patch 3: bottom      V=+1 on axis, -1 at corner
!    If INONUMAP=0 SR (radius) is a linear function of V
!    If INONUMAP=1 SR is quadratic in V with derivative=0 at V=-1
!-----------------------------------------------------------------------
      ELSEIF(IPI.EQ.2) THEN
         IF (INONUMAP==0) THEN
            SR=HALFRAD*(1.E0-V)
            XV3FAC=-HALFRAD
         ELSE
            ONEPMV=1.E0+V
            SR=QRAD*(4.E0-ONEPMV*ONEPMV)
            XV3FAC=-HALFRAD*ONEPMV
         ENDIF
         X(1)=SR*CTHETA
         X(2)=SR*STHETA
         X(3)=-DRAFT
         XU(1)=-PI*X(2)
         XU(2)= PI*X(1)
         XV(1)=XV3FAC*CTHETA
         XV(2)=XV3FAC*STHETA
!-----------------------------------------------------------------------
!  Patch 2: free surface   V=-1 RadiusI and V=+1 on side=RASIUS
!         r=(r2-r1)/2*(v+1)+r1 r=(r2-r1)/4*(v+1)**2+r1
!-----------------------------------------------------------------------
      ELSEIF(IPI.EQ.3) THEN
         IF (INONUMAP==0) THEN
            SR=HALFTHICK*(1.E0+V)+RADIUSI
            XV3FAC=HALFTHICK
         ELSE
            ONEPMV=1.E0+V
            SR=0.5*HALFTHICK*ONEPMV*ONEPMV+RADIUSI
            XV3FAC=HALFTHICK*ONEPMV
         ENDIF
         X(1)=SR*CTHETA
         X(2)=SR*STHETA
         XU(1)=-PI*X(2)
         XU(2)= PI*X(1)
         XV(1)=XV3FAC*CTHETA
         XV(2)=XV3FAC*STHETA
      ENDIF
      X(1:3)=X(1:3)+XSHIFT(1:3)
 99   RETURN
      END

      SUBROUTINE ELLIPSOID_CS(U,V,IBI,IPI,X,XU,XV,IGDFSCRB,IER)
!-----------------------------------------------------------------------
!
!     COPYRIGHT (C) 2006      WAMIT INCORPORATED
!
!-----------------------------------------------------------------------
!
!     Version : 6.3
!
!-----------------------------------------------------------------------
!
!     Source-code file :  GEOMXACT.F
!
!-----------------------------------------------------------------------
!
!     DESCRIPTION :  (IGDEF=-1003)
!
!
!     Subroutine defines the ellipsoid (IBI=403) by one patch
!     Center on free surface
!     Semi-axes A,B,C input from CSF file
!
!     Patch 1: body surface
!     Patch 2: free surface
!
!-----------------------------------------------------------------------
      IMPLICIT NONE
      REAL X(3),XU(3),XV(3),U,V
      INTEGER IPI,IBI,IGDFSCRB,IER
!-----------------------------------------------------------------------
!  user-assigned definitions -- inputs from GDF must have SAVE attribute
!  PI/4*(A,B,C) are evaluated once to reduce CPU time
!-----------------------------------------------------------------------
      REAL, SAVE :: A,B,C,PIO4A,PIO4B,PIO4C,HALFA,HALFB,AI,BI,CI,
     &              HALFAI,HALFBI,DIFA,SUMA,DIFB,SUMB,PIO4AI,PIO4BI
      REAL THETA,CTHETA,STHETA,PSI,CPSI,SPSI,SR,SRA,SRB
      REAL, PARAMETER :: PIO4=0.785398163397E0
      IF (IPI == 0) THEN
        READ (IGDFSCRB,*,IOSTAT=IER) A,B,C
        IF(IER/=0) GOTO 99
        READ (IGDFSCRB,*,IOSTAT=IER) AI,BI
        IF(IER/=0) GOTO 99
        PIO4A=PIO4*A
        PIO4B=PIO4*B
        PIO4C=PIO4*C
        HALFA=0.5E0*A
        HALFB=0.5E0*B
        PIO4AI=PIO4*AI
        PIO4BI=PIO4*BI
        HALFAI=0.5E0*AI
        HALFBI=0.5E0*BI
        DIFA=HALFA-HALFAI
        SUMA=HALFA+HALFAI
        DIFB=HALFB-HALFBI
        SUMB=HALFB+HALFBI
      ELSE
        THETA=PIO4*(U+1.E0)
        CTHETA=COS(THETA)
        STHETA=SIN(THETA)
      ENDIF
!-----------------------------------------------------------------------
!  Patch 1: body surface
!-----------------------------------------------------------------------
      IF (IPI == 1) THEN
         PSI=PIO4*(V+3.E0)
         CPSI=COS(PSI)
         SPSI=SIN(PSI)
         X(1)=A*SPSI*CTHETA
         X(2)=B*SPSI*STHETA
         X(3)=C*CPSI
         XU(1)=-PIO4A*SPSI*STHETA
         XU(2)= PIO4B*SPSI*CTHETA
         XV(1)=PIO4A*CPSI*CTHETA
         XV(2)=PIO4B*CPSI*STHETA
         XV(3)=-PIO4C*SPSI
!-----------------------------------------------------------------------
!  Patch 2: free surface
!-----------------------------------------------------------------------
      ELSEIF (IPI == 2) THEN
         SRA=DIFA*V+SUMA
         SRB=DIFB*V+SUMB
         X(1)=SRA*CTHETA
         X(2)=SRB*STHETA
         XU(1)=-PIO4*SRA*STHETA
         XU(2)= PIO4*SRB*CTHETA
         XV(1)=DIFA*CTHETA
         XV(2)=DIFB*STHETA
      ENDIF
 99   RETURN
      END

      SUBROUTINE ELLIPSOID_CS_NOSYM(U,V,IBI,IPI,X,XU,XV,IGDFSCRB,IER)
!-----------------------------------------------------------------------
!
!     COPYRIGHT (C) 2006      WAMIT INCORPORATED
!
!-----------------------------------------------------------------------
!
!     Version : 6.3
!
!-----------------------------------------------------------------------
!
!     Source-code file :  GEOMXACT.F
!
!-----------------------------------------------------------------------
!
!     DESCRIPTION :  (IGDEF=-1004)
!
!
!     Subroutine defines the ellipsoid (IBI=403) by one patch
!     Center on free surface
!     Semi-axes A,B,C input from CSF file
!
!     Patch 1: body surface
!     Patch 2: free surface
!
!-----------------------------------------------------------------------
      IMPLICIT NONE
      REAL X(3),XU(3),XV(3),U,V
      INTEGER IPI,IBI,IGDFSCRB,IER
!-----------------------------------------------------------------------
!  user-assigned definitions -- inputs from GDF must have SAVE attribute
!  PI/4*(A,B,C) are evaluated once to reduce CPU time
!-----------------------------------------------------------------------
      REAL, SAVE :: A,B,C,PIO4A,PIO4B,PIO4C,HALFA,HALFB,AI,BI,CI,
     &              HALFAI,HALFBI,DIFA,SUMA,DIFB,SUMB,PIO4AI,PIO4BI,
     &              XS,YS,ZS
      REAL THETA,CTHETA,STHETA,PSI,CPSI,SPSI,SR,SRA,SRB
      REAL, PARAMETER :: PIO4=0.785398163397E0,PI=3.141592654E0
      IF (IPI == 0) THEN
        READ (IGDFSCRB,*,IOSTAT=IER) A,B,C
        IF(IER/=0) GOTO 99
        READ (IGDFSCRB,*,IOSTAT=IER) AI,BI
        IF(IER/=0) GOTO 99
        READ (IGDFSCRB,*,IOSTAT=IER) XS,YS,ZS
        IF(IER/=0) GOTO 99
        PIO4A=PIO4*A
        PIO4B=PIO4*B
        PIO4C=PIO4*C
        HALFA=0.5E0*A
        HALFB=0.5E0*B
        PIO4AI=PIO4*AI
        PIO4BI=PIO4*BI
        HALFAI=0.5E0*AI
        HALFBI=0.5E0*BI
        DIFA=HALFA-HALFAI
        SUMA=HALFA+HALFAI
        DIFB=HALFB-HALFBI
        SUMB=HALFB+HALFBI
      ELSE
        THETA=PI*(U+1.E0)
        CTHETA=COS(THETA)
        STHETA=SIN(THETA)
      ENDIF
!-----------------------------------------------------------------------
!  Patch 1: body surface
!-----------------------------------------------------------------------
      IF (IPI == 1) THEN
         PSI=PIO4*(V+3.E0)
         CPSI=COS(PSI)
         SPSI=SIN(PSI)
         X(1)=A*SPSI*CTHETA
         X(2)=B*SPSI*STHETA
         X(3)=C*CPSI
         XU(1)=-PI*A*SPSI*STHETA
         XU(2)= PI*B*SPSI*CTHETA
         XV(1)=PIO4A*CPSI*CTHETA
         XV(2)=PIO4B*CPSI*STHETA
         XV(3)=-PIO4C*SPSI
!-----------------------------------------------------------------------
!  Patch 2: free surface
!-----------------------------------------------------------------------
      ELSEIF (IPI == 2) THEN
         SRA=DIFA*V+SUMA
         SRB=DIFB*V+SUMB
         X(1)=SRA*CTHETA
         X(2)=SRB*STHETA
         XU(1)=-PI*SRA*STHETA
         XU(2)= PI*SRB*CTHETA
         XV(1)=DIFA*CTHETA
         XV(2)=DIFB*STHETA
      ENDIF
      X(1)=X(1)+XS
      X(2)=X(2)+YS
      X(3)=X(3)+ZS
 99   RETURN
      END