SUBROUTINE DISP(U,KSTEP,KINC,TIME,NODE,NOEL,JDOF,COORDS)
C
      INCLUDE 'ABA_PARAM.INC'
C
      DIMENSION U(3),TIME(2),COORDS(3)
C
      IF(KSTEP.EQ.1) THEN
        U(1)=45.*TIME(1)
        U(2)=0.0
        U(3)=0.0
      ELSE
        TMP1=1.
        TMP4=4.
        PI=ATAN(TMP1)*TMP4
        TMP22=.22222222     
        TMP45=45.
        U(1)=2.*SIN(TMP22*PI*TIME(1))+TMP45
        U(2)=(2.*TMP22*PI)*COS(TMP22*PI*TIME(1))
        U(3)=-(2.*(TMP22*PI)**2)*SIN(TMP22*PI*TIME(1))
      END IF
      RETURN
      END
C     
C**************************
C  USER DEFINED WAVE THEORY
C**************************
C     
      subroutine uwave(v,a,pdyn,dpdyndz,surf,lpdyn,
     1        lrecompute,luplocal,lupglobal,
     1        lsurf,ndim,xcur,xintmed,
     2        grav,density,elevb,elevs,
     3        seed,nspectrum,wamp,
     2        time,dtime,noel,npt,kstep,kinc)              
C                                                                               
      include 'aba_param.inc'
C                         
      dimension v(ndim),a(ndim),xcur(ndim),xintmed(ndim)
      dimension time(2),wamp(2,nspectrum)
  
C
      parameter(pi=3.14159265358979d0,two=2.0d0,abig=1.d36)
      parameter (const2=2.d10, twopi = 2.d0*pi )
C
      luplocal=0
      lupglobal=0
      if(lrecompute.ne.0) then
C      Only stochastic analysis with UWAVE can have lrecompute=1
C      The user must set other flags accordingly; see User's manual.
      else
C      For regular Aqua analysis with UWAVE, lrecompute=0 always
C    
C
C   Wave definition for a single Airy wave component:
C      Phase angle of waves: in radians
       phase=-54.026d0*pi/180.d0
C
C      Wave travel direction:
       xdir=1.0d0
       ydir=0.0d0
C
C      Period, wavelength, wave number, wave height, frequency:
C
       period=9.d0
       waveln=415.0d0
       wavenum=twopi/waveln
       waveamp=10.0d0
       freq=twopi/period
C
       if (lsurf.eq.1) then
C         Calculate the instantaneous water surface only, no
C         wave kinematics are required:

          wtp=-freq*time(2)+phase
          sn=xdir*xcur(1)
          if (ndim.eq.3) sn=sn+ydir*xcur(2)
          termt=wavenum*sn+wtp
          surf=elevs-waveamp*cos(termt)

       else

C         Calculate wave kinematics:
C
C             lpdyn = 0 - velocity and acceleration for
C                         drag or inertia loads
C             lpdyn = 1 - dynamic pressure and gradient of
C                         dynamic pressure for buoyancy loads


C         For AIRY waves -
C            Extrapolation is used above the mean water line.
C            Therefore, if the current integration point is above the
C            mean water line, move it back to the still water surface
C            in order to calculate the wave kinematics:
C 

          if (ndim.eq.3) then
            if (xcur(3).gt.elevs) xcur(3)=elevs
          else
            if (xcur(2).gt.elevs) xcur(2)=elevs
          endif

C         Calculate some terms needed in the computations:
          tacc=twopi*grav
          rhog = - density*grav
          twopirhog = twopi * rhog
          w10=tanh(twopi*(elevs-elevb)/waveln)
          wtp=-freq*time(2)+phase
          sn=xdir*xcur(1)
          if (ndim.eq.3) sn=sn+ydir*xcur(2)
          termt=wavenum*sn+wtp
C     
C         Watch out for overflow:
C     
          if (abs(xcur(ndim)-elevb).gt.
     1        abs((elevs-elevb)/two)) then
             termz=wavenum*(xcur(ndim)-elevs)
             if (abs(termz) .gt. log(const2)) then
                coshm = const2/two
                sinhm = sign(const2/two,termz)
                coshz = coshm +w10*sinhm
                sinhz = sinhm +w10*coshm
             else
                coshm = cosh(termz)
                sinhm = sinh(termz)
                coshz = coshm +w10*sinhm
                sinhz = sinhm +w10*coshm
             end if
          else
             termz=wavenum*(xcur(ndim)-elevb)
             term2 = wavenum*(elevs-elevb)
             if (abs(term2) .gt. log(abig)) then
                term2 = log(abig)
             end if
             if (abs(termz) .gt. log(const2)) then
                coshm = const2/two
                sinhm = sign(const2/two,termz)
                coshz = coshm/cosh(term2)
                sinhz = sinhm/cosh(term2)
             else
                coshz = cosh(termz)/cosh(term2)
                sinhz = sinh(termz)/cosh(term2)
             end if
          end if
C
          cost=cos(termt)
          sint=sin(termt)
          tcoshcos=coshz*cost
          if (lpdyn.eq.0) then

C            Drag or inertia loads:

C            - g * tauN * aN * amp / lambdaN
             term1=-grav*period*waveamp/waveln
             term11=term1*tcoshcos
C            Velocity:
             v(1)=term11*xdir + v(1)
             if (ndim.eq.3) v(2)=term11*ydir + v(2)
             v(ndim)=term1*sinhz*sint + v(ndim)
             term1=tacc*waveamp/waveln
             term11=term1*coshz*sint
C            Acceleration:
             a(1)=-term11*xdir
             if(ndim.eq.3) a(2)=-term11*ydir
             a(ndim)=term1*sinhz*cost
          else
C            Dynamic pressure for buoyancy loads:
C            -rho*g*aN*cosh*cos*amp:
             pDyn  = rhog*waveamp*tcoshcos

C            Gradient of the dynamic pressure,
C            -2pi*rho*g*aN/LambdaN*sinh*cos*amp 
             DpdynDz = twopirhog*waveamp/waveln*sinhz*cost

C            Airy wave extrapolation:  Above the mean water line
C            the gradient of the dynamic pressure is zero.
             if (ndim.eq.3) then
                if (xcur(3).eq.elevs) dpdyndz=0.0d0
             else
                if (xcur(2).eq.elevs) dpdyndz=0.0d0
             end if
          end if
       endif
      end if
C
      return
      end