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 xindir*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 lxindir = 0 xfname =' ' xindir =' ' call getjobname(xfname,lxfname) !input file name call getoutdir(xindir,lxindir) !input 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 results file C The process involves transfering each record in the current C uwavexx2.fil to the new uwavexx2.fin file, modifying the C NFORC record in Step 2 to store the intermed config data. C noel1=noel kinc1=kinc npt1=npt if(lrdflg.ne.1) then FNAMEX=dmkname(xfname(1:lxfname),xindir(1:lxindir),' ') LRUNIT(1,1)=8 C ! Fortran unit 8 ! LRUNIT(2,1)=2 C ! Binary format ! Input LOUTF=2 C ! Binary format ! Output INRU=1 NRU=1 C ! to read only one file ! CALL INITPF (FNAMEX, NRU, LRUNIT, LOUTF) JUNIT = LRUNIT(1,1) CALL DBRNU (JUNIT) I2001 = 0 KEYPRV = 0 KSTEP1 = 0 C C Covers a maximum of 10 million records in each file. C DO 80 IXX2 = 1, 100 DO 80 IXX = 1, 99999 jrcd=0 CALL DBFILE(0,ARRAY,JRCD) WRITE(6,*) 'KEY/RECORD LENGTH = ', JRRAY(1,2),JRRAY(1,1) IF (JRCD .NE. 0 .AND. KEYPRV .EQ. 2001) THEN WRITE(6,*) 'END OF FILE #', INRU CLOSE (JUNIT) GOTO 100 ELSE IF (JRCD .NE. 0) THEN WRITE(6,*) 'ERROR READING FILE #', FNAMEX CLOSE (JUNIT) GOTO 110 ENDIF C C Initialize the flag to write a record to the file: C LWRITE=0 -- write disabled C LWRITE=1 -- write enabled C LWRITE=1 C C If this is the first input file, or if the write flag has not been C disabled for records in subsequent files, then write the data to C the output file. C IF (INRU .EQ. 1 .OR. LWRITE .EQ. 1) THEN KEY=JRRAY(1,2) if(key.eq.2000) then C Get step and inc number kstep1=jrray(1,8) kinc1=jrray(1,9) npt1=1 end if if(key.eq.1) then C Get elem number noel1=jrray(1,3) end if if(key.eq.15 .and. kstep1.eq.2) then C Get node number and set NFORC record to zero nnum1=jrray(1,3) call initialize(array(4),zero,jrray(1,1)-3) jrray(1,3)=npt1 call copy(xintm(1,npt1,noel1,kinc1),array(4),ndim) npt1=npt1+1 if(npt1.gt.2) npt1=npt1-2 end if CALL DBFILW(1,ARRAY,JRCD) IF (JRCD .NE. 0) THEN call stdb_abqerr(-2,'Error writing file %S with '// $ 'unit number %I ',junit,zero,xfname(1:lxfname)) CLOSE (JUNIT) GOTO 110 ENDIF ENDIF KEYPRV = JRRAY(1,2) 80 CONTINUE go to 201 100 CONTINUE 110 CONTINUE 201 continue lrdflg=1 C Done - That's it for Step 3 C call xit 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 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