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 common/crdflg/xintm(2,2,10,144),lrdflg,lwrtflg C 2-ndim x (2-node x 10-elem) x 144-inc DIMENSION ARRAY(513), JRRAY(NPRECD,513) EQUIVALENCE (ARRAY(1), JRRAY(1,1)) CHARACTER xoutdir*255, xfname*80 CHARACTER dmkname*255, FNAMEX*80 PARAMETER (MXUNIT=2, zero=0.d0) INTEGER LRUNIT(2,MXUNIT),LUNIT(10) 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 lxfname = 0 lxoutdir = 0 xfname =' ' xoutdir =' ' call getjobname(xfname,lxfname) !input file name call getoutdir(xoutdir,lxoutdir) !output directory luplocal=0 lupglobal=0 if(kstep.eq.3) then C Step 3 is a dummy step set up to C Copy intermed config data to the external file uwavexx3.017. C These data are reprocessed in the .fil file post-processing C program uwavexx4.f and stored as NFORC records to facilitate QA. C noel1=noel kinc1=kinc npt1=npt if(lrdflg.ne.1) then fnamex=dmkname(xfname(1:lxfname),xoutdir(1:lxoutdir),'.017') open(unit=17,file=fnamex,status='unknown', 1 form='unformatted') write(17)(xintm(ix,1,1,1),ix=1,5760) close(17) lrdflg=1 end if C else 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. write(7,*) 'time = ',time write(7,*) 'increment number = ',kinc write(7,*) 'LRECOMPUTE, LUPLOCAL, LUPGLOBAL ', 1 LRECOMPUTE, LUPLOCAL, LUPGLOBAL write(7,*) 'element, point = ', noel, npt write(7,*) 'current - ',(xcur(ii), ii=1,ndim) write(7,*) 'intermed - ',(xintmed(ii), ii=1,ndim) write(7,*) ' ' if(kstep.eq.2 ) then C Write user specifed data to jobid.96 file for verification - if(lwrtflg.eq.1 .and. (kinc.eq.1 .or. kinc.eq.144)) then write(96,111)kstep,kinc,noel,npt,seed,wamp 111 format(' From uwave-kstep,kinc,noel,npt,seed,wamp :', 1 4i5,1pg10.3,/,4(1pg10.3,',',2x) ) end if C Exercise global and local update requests - if(kinc.eq.1 .or.kinc.eq.141) then C Global Update of interm config to current config lupglobal=1 else if(kinc.eq.11 .and. noel.eq.1) then C Local Update of interm config to current config: Elem 1 lupglobal=0 luplocal=1 else if(kinc.eq.21 .and. noel.eq.2) then C Local Update of interm config to current config: Elem 2 lupglobal=0 luplocal=1 else if(kinc.eq.31 .and. noel.eq.3) then C Local Update of interm config to current config: Elem 3 lupglobal=0 luplocal=1 else if(kinc.eq.41 .and. noel.eq.4) then C Local Update of interm config to current config: Elem 4 lupglobal=0 luplocal=1 else if(kinc.eq.51 .and. noel.eq.5) then C Local Update of interm config to current config: Elem 5 lupglobal=0 luplocal=1 else if(kinc.eq.61 .and. noel.eq.6) then C Local Update of interm config to current config: Elem 6 lupglobal=0 luplocal=1 else if(kinc.eq.71 .and. noel.eq.7) then C Local Update of interm config to current config: Elem 7 lupglobal=0 luplocal=1 else if(kinc.eq.81 .and. noel.eq.8) then C Local Update of interm config to current config: Elem 8 lupglobal=0 luplocal=1 else if(kinc.eq.91 .and. noel.eq.9) then C Local Update of interm config to current config: Elem 9 lupglobal=0 luplocal=1 else if(kinc.eq.101 .and. noel.eq.10) then C Local Update of interm config to current config: Elem 10 lupglobal=0 luplocal=1 end if C Store intermed config to temporary array- call copy(xintmed(1),xintm(1,npt,noel,kinc),ndim) end if C 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 wavehgt=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-wavehgt*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*wavehgt/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*wavehgt/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*wavehgt*tcoshcos C Gradient of the dynamic pressure, C -2pi*rho*g*aN/LambdaN*sinh*cos*amp DpdynDz = twopirhog*wavehgt/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 (lrecompute.ne.0) end if C (kstep.eq.3) C return end C User subroutine to manipulate user external files subroutine uexternaldb(lop,lrestart,time,dtime,kstep,kinc) C include 'aba_param.inc' C Common block used in uwave.f for storage of intermediate configuration common/crdflg/xintm(2,2,10,144),lrdflg,lwrtflg C character xoutdir*255, xfname*80 character dmkname*255, fnamex*80 C parameter(nspectrum=6, zero=0.d0) C Dummy wave spectrum data - 6 (freq,ampl) data pairs dimension wamp(2,nspectrum) data wamp/0.,0., 0.05,1.0, 0.075,120., 0.12,20., 1 0.2,1.0, 0.5,0.0/ dimension time(2) C lxfname = 0 lxoutdir = 0 xfname =' ' xoutdir =' ' call getjobname(xfname,lxfname) !input file name call getoutdir(xoutdir,lxoutdir) !output directory if(lop.eq.0) then C C During first visit: create, open and write to file jobid.96 C via Fortran Unit 96. C Fortran unit advisory: Fortran unit should be greated than 100 C to avoid conflicts with internal Fortran unit numbers C The data written to the file jobid.96 are identical to the C dummy wave spectrum data specified for this analysis. C fnamex=dmkname(xfname(1:lxfname),xoutdir(1:lxoutdir),'.96') open(unit=96,file=fnamex,status='unknown',form='formatted') write(96,*)'Opening new user external file...' write(96,*)'Writing dummy data to this file...' write(96,*)' ' write(96,*)'Wave spectrum data (freq vs. ampl)' write(96,*)' ' write(96,1) wamp 1 format(4(1pg10.3,',',2x)) C C Initialize read and print flags and common block used in wave.f C lrdflg=0 lwrtflg=1 call initialize(xintm,zero,5760) write(96,2)LOP,LRESTART,KSTEP,KINC,TIME,DTIME 2 format('LOP,LRESTART,KSTEP,KINC,TIME,DTIME: ', 1 4i5,1p3g11.3) C else if(lop.eq.1 .or. lop.eq.2) then C C write(96,*)'LOP = 1 or 2 = beginning or end of inc: ' write(96,2)LOP,LRESTART,KSTEP,KINC,TIME,DTIME if(lrestart.eq.1) then write(96,*)' Restart file will be written for this inc:' else if(lrestart.eq.2) then write(96,*)' Restart file will be OVERLAYed for this inc:' else if(lrestart.eq.3) then write(96,*)' Beginning of a restart analysis: ' end if else if(lop.eq.3) then C C During last visit: close file; closed file will be C Copied to job origin directory C write(96,*)'Closing external file ( unit 96)...' write(96,2)LOP,LRESTART,KSTEP,KINC,TIME,DTIME close(96) else if(lop.eq.4) then C C write(96,*)'LOP = 4 = beginning of restart analysis: ' write(96,2)LOP,LRESTART,KSTEP,KINC,TIME,DTIME end if C return end C subroutine copy(a,b,n) C include 'aba_param.inc' dimension a(*),b(*) C do 1 i=1,n 1 b(i)=a(i) return end C subroutine initialize(a,b,n) C include 'aba_param.inc' dimension a(*) C do 1 i=1,n 1 a(i)=b return end c c Compose a filename directory/jobname.exten character*(*) function dmkname(fname,dname,exten) C character*(*) fname,dname,exten C fname I jobname C dname I directory C exten I extension C dmkname O directory/jobname.exten ltot = len(fname) lf = 0 do k1 = ltot,2,-1 if (lf.eq.0.and.fname(k1:k1).ne.' ') lf = k1 end do ltot = len(dname) ld = 0 do k1 = ltot,2,-1 if (ld.eq.0.and.dname(k1:k1).ne.' ') ld = k1 end do ltot = len(exten) le = 0 do k1 = ltot,2,-1 if (le.eq.0.and.exten(k1:k1).ne.' ') le = k1 end do if ((lf + ld + le) .le. len(dmkname)) then dmkname = dname(1:ld)//'/'//fname(1:lf) ltot = ld + lf + 1 if ( le.gt.0) then dmkname = dmkname(1:ltot)//exten(1:le) end if end if C return end