C*********************************************************************** C* * C* N-linear Stream Function Theory * C* * C* This program computes the symmetric nonlinear wave travelling * C* with constant form on a shear current represented by n straight * C* line segments. * C* * C*********************************************************************** PARAMETER(IX=104,IC=62,NNI=3) DIMENSION X(IX),XF1(IC,IX),XF2(IC,IX),U(0:NNI+1),UIN(0:NNI+1) DIMENSION THETA(IC),ZETA(NNI+1,IC),D(0:NNI+1) COMMON X,XF1,XF2 COMMON /DATA/NN,NTHTS,QBAR,T,D,H,U,UIN,THETA,DTHETA,NI, 1ZETA,HT PI=3.1415927 WRITE(*,3453) 3453 FORMAT('********************************************************** 1****'/ 2'* *'/ 3'*--- N - L I N E A R W A V E T H E O R Y ----*'/ 4'* *'/ 5'* *'/ 6'* Developed by: Robert A. Dalrymple *'/ 7'* 18 Harvest Lane *'/ 8'* Newark, DE 19711 *'/ 2'* References: *'/ 9'* *') WRITE(*,3454) 3454 FORMAT('* 1 *'/ 2'* Dalrymple, R.A., "Water Waves on a Bilinear *'/ 3'* Shear Current," Proc. 14th Coastal Engrg. Conf, *'/ 4'* ASCE, 1974. *'/ 5'* *'/ 6'* Dalrymple, R.A. and J.C. Heidemann, "Nonlinear *'/ 7'* Water Waves on a Vertically-sheared Current," *'/ 8'* Proc. Workshop on Wave and Current Kinematics *'/ 9'* and Loading, Oil Industry Intl. Exploration *') WRITE(*,3455) 3455 FORMAT('* Production Forum, Rpt. 3.12/156, London 1 *'/ 2'* 69-92, 1989. *'/ 3'* *'/ 4'* Program Input *'/ 5'* *'/ 6'* Wave Height, Period, Water Depth, Interface Depths *'/ 6'* Number of interfaces *'/ 7'* N-linear Current: Surface, Interface and Bottom *'/ 8'* Velocities *'/ 3'* Theory Order, Number of Iterations (try 30) *'/ 4'* *'/ 9'**************************************************************') 888 WRITE(*,*) ' INPUT WAVE HEIGHT, PERIOD, WATER DEPTH' READ(*,*) HT,T,DPT WRITE(*,*) ' INPUT NUMBER OF INTERFACES' READ(*,*) NI IF(NI.GT.NNI) THEN WRITE(*,*) '** NUMBER OF INTERFACES EXCEEDS PROGRAM CAPACITY**' WRITE(*,*) ' INCREASE PARAMETER NNI THROUGHOUT PROGRAM' GO TO 888 ENDIF WRITE(*,*) ' INPUT SURFACE AND BOTTOM VELOCITIES' READ(*,*) US,UB C C DEFINE VELOCITIES AT BOTTOM, INTERFACES AND THE SURFACE C U IS INCREMENTAL CHANGE IN VELOCITY WHILE UIN IS TOTAL C VELOCITY. C U(0)=UB UIN(0)=UB UIN(NI+1)=US D(0)=DPT D(NI+1)=0. WRITE(*,*) ' FROM THE LOWEST TO THE UPPERMOST INTERFACE:' WRITE(*,*) ' INPUT DEPTH, VELOCITY' DO 7 IN=1,NI READ(*,*) D(IN),UTEMPN UIN(IN)=UTEMPN U(IN)=UTEMPN-UIN(IN-1) UTEMP=UTEMPN 7 CONTINUE U(NI+1)=US-UIN(NI) WRITE(*,*) 'DEPTHS, VELOCITIES AND INCREMENTAL CHANGES' WRITE(*,107) (D(J),UIN(J),U(J),J=0,NI+1) 107 FORMAT(3F10.2) NTHTS=31 WRITE(*,*) ' INPUT WAVE THEORY ORDER, NUMBER OF ITERATIONS' READ(*,77) NORDER,KMAX 77 FORMAT(2I5,A60) NN=NORDER*2+2 IF((NI*NN+NN/2+3).GT.IX) THEN WRITE(*,*) ' *** WAVE THEORY ORDER EXCEEDS PROGRAM CAPACITY ***' WRITE(*,*) ' INCREASE PARAMETER IX THROUGHOUT PROGRAM ' C NOTE THAT BRACKET TERM IS THE MAX. NO. OF COEFFS + MULTIPLIERS GO TO 888 ENDIF WRITE(*,*) ' INPUT INITIAL DAMPING, INCREMENT XS(Try 1, or 0)' READ(*,*) DAMP,INC AUTOFF=1 100 FORMAT(2I3,7F6.2,2I3,6X,5A4) XN=NTHTS-1 DTHETA=PI/XN IF(KMAX.EQ.0) KMAX=25 WRITE(6,733) 733 FORMAT('1', ' THEORETICAL WAVE ON N-LINEAR SHEAR CURRENT: ') WRITE(6,102) HT,T,DPT 102 FORMAT(///' WAVE HEIGHT = ',F6.2,' FEET, PERIOD = ',F6.2,' SECONDS 1'/' WATER DEPTH = ',F6.2,' FEET') WRITE(6,103) NORDER,KMAX 103 FORMAT(/' WAVE THEORY ORDER = ',I3,', TOTAL NUMBER OF ITERATIONS A 1LLOWED =',I3) THETA(1)=0.0 DO 10 I=2,NTHTS 10 THETA(I)=THETA(I-1)+DTHETA HBT2=HT/(T*T) DBT2=DPT/(T*T) IF(DBT2.LT.0.075.AND.DAMP.GT.0.3) DAMP=.3 C DETERMINE BREAKING INDEX HBT2T=.873*SINH(.8935*DBT2)/COSH(.8935*DBT2) IF((HBT2/HBT2T).GT.0.75.AND.DAMP.GT.0.3) DAMP=.3 IF((HBT2-.01-HBT2T).GT.0.0) GO TO 6 IF(ABS(HBT2-HBT2T)-.01) 3,3,15 3 WRITE(6,153) WRITE(6,201) 153 FORMAT(' ********************************************************* 1*******************************************************') 201 FORMAT(' **** INPUT WAVE PARAMETERS ARE WITHIN THE ACCURACY OF TH 1E BREAKING INDICATOR; IF SOLUTION IS NOT SUCCESSFUL,'/' CHEC 2K INPUT PARAMETERS. ****' ) WRITE(6,153) GO TO15 6 WRITE(6,153) HB=HBT2T*T*T WRITE(6,200) HB 200 FORMAT(' WAVE HEIGHT EXCEEDS BREAKING WAVE HEIGHT OF APPROXIMATE 1LY',F8.2,' FEET.'/' CHECK INPUT PARAMETERS.' ) WRITE(6,153) GO TO 888 15 CONTINUE 35 CALL WVLEN(DPT,T,X(1)) X(1)=1.1*X(1) C=X(1)/T NX=NI*NN+NN/2-NI+1 DO 5 N=2,NX 5 X(N)=0. IF(INC.EQ.1) THEN NNT=6 ELSE NNT=NN ENDIF A1=2.0*PI*DPT/X(1) A=-C*HT/(2.0*SINH(2.0*PI*(DPT+HT/2.0)/X(1))) X(2)=0. X(3)=A*COSH(A1) X(4)=A*SINH(A1) DO 545 IN=1,NI NINTER=(NI-IN+1)*NNT+(IN+1-NI) X(NINTER)=-D(IN)*((C-UIN(IN))-U(IN+1)*(D(IN)/2.)/(D(IN)-D(IN+1))) 545 CONTINUE WRITE(*,101)(X(I),I=1,NX) 101 FORMAT(//' INITIAL STREAM FUNCTION COEFFICIENTS:',6E12.3/ 1(10E12.3)) CALL STREAM(KMAX,DAMP,*35,*888,AUTOFF,INC) CALL VELOC GO TO 888 8000 STOP END C*********************************************************************** SUBROUTINE STREAM(KMAX,DAMP,*,*,AUTOFF,INC) PARAMETER(IX=104,IC=62,NNI=3) DIMENSION X(IX),XF1(IC,IX),XF2(IC,IX),U(0:NNI+1),UIN(0:NNI+1) DIMENSION THETA(IC),ZETA(NNI+1,IC),D(0:NNI+1),XN(IX) COMMON X,XF1,XF2 COMMON /DATA/NN,NTHTS,QBAR,T,D,H,U,UIN,THETA,DTHETA,NI, 1ZETA,HT DO 6 I=1,NTHTS DO 6 N=1,(NN/2-1) XF1(I,N)=COS(THETA(I)*FLOAT(N)) 6 XF2(I,N)=SIN(THETA(I)*FLOAT(N)) INP1=NI+1 C INCREMENT X'S FROM SECOND ORDER UP TO SPECIFIED ORDER IF INC=1 NFINAL=NN IF(INC.EQ.1) THEN NN=6 ENDIF CALL SCALC(TOTERR,XH) DO 20 KI=1,KMAX CALL CFF(DAMP,0,KI,TOTERR,XH) CALL SCALC(TOTERR,XH) IF(KI.LT.((NFINAL-2)/2).OR.INC.NE.1) GO TO 321 C C PROCEDURE TO INCREMENT COEFFICIENTS C INCREMENT NO. OF COEFFIENTS AND ROLL DOWN IF(NN.LT.NFINAL) THEN NX=NI*(NN+2)+NN/2+3 DO 55 N=1,NX 55 XN(N)=0.0 DO 155 N=1,NN 155 XN(N)=X(N) NLAST=(NN-2)/2+1 DO 221 IN=1,NI NCB=NI-IN+1 NCOEF=NCB*NN+(IN+1-NI)+NLAST DO 223 L=1,NLAST NCOEF=NCOEF-1 XN(NCOEF+2*NCB)=X(NCOEF) 223 CONTINUE NLAST=NN-1 221 CONTINUE NN=NN+2 DO 66 N=1,NX 66 X(N)=XN(N) ENDIF 321 IF(AUTOFF.EQ.1) GO TO 21 IF(TOTERR.LT..01) RETURN 21 IF(KI.EQ.1) E1=TOTERR IF(TOTERR.GT.E1) GO TO 25 20 CONTINUE CALL CFF(DAMP,1,KMAX,TOTERR,XH) RETURN 25 DAMP=DAMP-.1 IF(DAMP.LE.0.0) GO TO 30 WRITE(6,110) 110 FORMAT(' SOLUTION DIVERGING---PROGRAM IS RESTARTING WITH MORE DAM 1PING') WRITE(6,115) DAMP 115 FORMAT(' DAMPING = ', F5.2) RETURN 1 30 WRITE(6,120) 120 FORMAT(//' *****************************************'/' ********** 1 CONVERGENCE IMPOSSIBLE *****'/' ********************************* 2*************'/' CHECK INPUT VARIABLES ') RETURN 2 END C*********************************************************************** SUBROUTINE SCALC(ERRORH,XH) PARAMETER(IX=104,IC=62,NNI=3) DIMENSION X(IX),XF1(IC,IX),XF2(IC,IX),U(0:NNI+1),UIN(0:NNI+1) DIMENSION THETA(IC),ZETA(NNI+1,IC),D(0:NNI+1) COMMON X,XF1,XF2 COMMON /DATA/NN,NTHTS,QBAR,T,D,H,U,UIN,THETA,DTHETA,NI, 1ZETA,HT NNP1=NN+1 DO 5 IN=1,NI NINTER=(NI-IN+1)*NN+IN-NI+1 CALL SURFAS(IN,NINTER,1) 5 CONTINUE NIP1=NI+1 C COMPUTE SURFACE DISPLACEMENTS CALL SURFAS(NIP1,2,1) XH=ZETA(NIP1,1)-ZETA(NIP1,NTHTS) ERRORH=ABS(XH-HT) WRITE(6,105)(THETA(I),ZETA(NIP1,I),I=1,NTHTS) 105 FORMAT(' THETAS AND CALCULATED ETAS'/(1P8E10.3)) DO 56 IN=NI,1,-1 WRITE(6,106)(THETA(I),ZETA(IN,I),I=1,NTHTS) 106 FORMAT(' THETAS AND CALCULATED ZETAS'/(1P8E10.3)) 56 CONTINUE RETURN END C*********************************************************************** SUBROUTINE SURFAS(IN,NINTER,ITP) C C COMPUTE THE INTERFACIAL ELEVATIONS USING NEWTON-RAPHSON METHOD AND C THE TOTAL ELEVATION, TOTH C C*********************************************************************** PARAMETER(IX=104,IC=62,NNI=3) DIMENSION X(IX),XF1(IC,IX),XF2(IC,IX),U(0:NNI+1),UIN(0:NNI+1) DIMENSION THETA(IC),ZETA(NNI+1,IC),D(0:NNI+1) COMMON X,XF1,XF2 COMMON /DATA/NN,NTHTS,QBAR,T,D,H,U,UIN,THETA,DTHETA,NI, 1ZETA,HT CON=3.1415927/X(1) C=X(1)/T PSIY=X(NINTER) IF(NINTER.NE.2) NINTER=NINTER-(NN-1) DPT=D(IN) NIPTS=NTHTS IF (ITP.EQ.2) THEN NIPTS=1 ZETA(IN,I)=0. ENDIF DO 21 I=1,NIPTS YI=ZETA(IN,I) TOTH=YI-DPT DO 20 IT=1,20 SUM=0.0 SUM1=0.0 DO 18 N=1,NN-2,2 XKN=FLOAT(N+1)*CON L=(N+1)/2 ZN=XKN*TOTH IF(ABS(ZN).LT. 15.0 ) GO TO 118 WRITE(6,29) I,IN 29 FORMAT(' SINH OR COSH ARG. TOO LARGE FOR THETA ',2I4) GO TO 21 118 NK=NINTER+N SUM=SUM+(X(NK)*SINH(ZN)+X(NK+1)*COSH(ZN))*XF1(I,L) SUM1=SUM1+(X(NK)*COSH(ZN)+X(NK+1)*SINH(ZN))*XF1(I,L)*XKN 18 CONTINUE IF(IN.EQ.(NI+1))THEN UTOT=UIN(IN-1)-C US=U(IN)/(D(IN-1)-D(IN)) F=PSIY+UTOT*TOTH+US*(TOTH*TOTH/2.+D(IN-1)*TOTH)-SUM FP=UTOT+US*(TOTH+D(IN-1))-SUM1 ELSE UTOT=UIN(IN)-C US=U(IN+1)/(D(IN)-D(IN+1)) F=PSIY+UTOT*TOTH+US*(TOTH*TOTH/2.+D(IN)*TOTH)-SUM FP=UTOT+US*(TOTH+D(IN))-SUM1 ENDIF TOTHN=TOTH-F/FP IF(ABS((TOTHN-TOTH)).LT.0.0001) THEN GO TO 25 ELSE TOTH=TOTHN ENDIF 20 CONTINUE WRITE(*,*) 'TWENTY ITERATIONS FAILED FOR INTERFACE',IN 25 ZETA(IN,I)=TOTHN+DPT 21 CONTINUE RETURN END C*********************************************************************** SUBROUTINE CFF(DAMP,IFSET,KI,TOTERR,XH) PARAMETER(IX=104,IC=62,NNI=3) DIMENSION X(IX),XF1(IC,IX),XF2(IC,IX),U(0:NNI+1),UIN(0:NNI+1) DIMENSION THETA(IC),ZETA(NNI+1,IC),D(0:NNI+1),B(IX*IX) COMMON X,XF1,XF2 COMMON /DATA/NN,NTHTS,QBAR,T,D,H,U,UIN,THETA,DTHETA,NI, 1ZETA,HT COMMON /PARAM/DCDX(IX),DETADX(IC,IX),DUDX(IX),DVDX(IX),DQBAR(IX) 1,DQDX(IC,IX),Q(IC),ERROR(IC),UI1(NNI,IC),UI2(NNI,IC) COMMON /PARAM2/ERROR2(NNI,IC),DU1DX(NNI,IC,IX),DU2DX(NNI,IC,IX), 1DZETDX(NNI,IC,IX),A(IX),DV1DX(NNI,IC,IX),DV2DX(NNI,IC,IX), 2VI2(NNI,IC),VI1(NNI,IC),ERROR3(NNI,IC) COMMON /PARAM3/DEPTHI(NNI+1),DU1DZ(NNI),DU2DZ(NNI),DV2DZ(NNI), 1DV1DZ(NNI),DZETADX(NNI,IC,IX),UZETA(NNI) CON=3.1415927/X(1) CO=6.2831853/X(1) DTHETL=DTHETA/3.1415927 G=32.17 INP1=NI+1 XNTHTS=NTHTS NTHTM1=NTHTS-1 NLAST1=(NN-2)/2 NLAST2=NN-2 C MAXIMUM NO. OF EQUATIONS, NMAX+NI+2 NX=NI*NN+NN/2+3 NX2=NX*NX DO 551 L=1,NX2 551 B(L)=0.0 C NO. OF LAGRANGE MULTIPLIERS NL=NI+2 C MAXIMUM NO. OF X COEFFICIENTS, NXML NXML=NX-NL NXMLP1=NXML+1 C=X(1)/T C INITIALIZE ALL VARIABLES TO ZERO DO 1 N=1,NX DCDX(N)=0.0 DQBAR(N)=0.0 DO 1 I=1,NTHTS DQDX(I,N)=0.0 DETADX(I,N)=0.0 DO 1 IN=1,NI DZETADX(IN,I,N)=0.0 DU1DX(IN,I,N)=0. DU2DX(IN,I,N)=0. DV1DX(IN,I,N)=0. DV2DX(IN,I,N)=0. 1 CONTINUE C INITIALIZE ADDITIONAL VARIABLES DCDX(1)=1.0/T USURF=U(NI+1)/(D(NI)-D(NI+1)) C DO 20 I=1,NTHTS ETA=ZETA(NI+1,I) DETADX(I,1)=ETA/T DETADX(I,2)=-1.0 UU=0.0 VV=0.0 DUDX(1)=0.0 DVDX(1)=0.0 DUDX(2)=0.0 DVDX(2)=0.0 DUDETA=USURF DVDETA=0.0 DO 5 IN=1,INP1 UI1(IN,I)=0.0 UI2(IN,I)=0.0 VI1(IN,I)=0.0 VI2(IN,I)=0.0 DU1DX(IN,I,1)=-1./T DU2DX(IN,I,1)=-1./T DU2DZ(IN)=U(IN)/(D(IN-1)-D(IN)) IF(IN.GT.NI) GO TO 155 DU1DZ(IN)=U(IN+1)/(D(IN)-D(IN+1)) 155 CONTINUE DV1DZ(IN)=0.0 DV2DZ(IN)=0.0 DEPTHI(IN)=-D(IN)+ZETA(IN,I) DZETADX(IN,I,1)=DEPTHI(IN)/T NCB=(NI-IN+1)*NN+(IN+1-NI) DZETADX(IN,I,NCB)=-1. 5 CONTINUE C COMPUTE SURFACE VALUES OF VELOCITY DO 14 N=1,NLAST2,2 XN=N+1 ZN=CON*XN L=(N+1)/2 SHS=SINH(ZN*ETA) CHS=COSH(ZN*ETA) COEF1=X(N+2)*CHS+X(N+3)*SHS COEF2=X(N+2)*SHS+X(N+3)*CHS UU=UU-ZN*COEF1*XF1(I,L) VV=VV-ZN*COEF2*XF2(I,L) DETADX(I,1)=DETADX(I,1)-ZN*ETA*COEF1*XF1(I,L)/X(1) DETADX(I,N+2)=SHS*XF1(I,L) DETADX(I,N+3)=CHS*XF1(I,L) DUDX(1)= DUDX(1)+(ZN* COEF1+ZN*ZN*ETA*COEF2)*XF1(I,L)/X(1) DUDX(N+2)=-ZN*CHS*XF1(I,L) DUDX(N+3)=-ZN*SHS*XF1(I,L) DVDX(1)=DVDX(1)+(ZN*COEF2+ZN*ZN*ETA*COEF1)*XF2(I,L)/X(1) DVDX(N+2)=-ZN*SHS*XF2(I,L) DVDX(N+3)=-ZN*CHS*XF2(I,L) DUDETA=DUDETA-ZN*ZN*COEF2*XF1(I,L) DVDETA=DVDETA-ZN*ZN*COEF1*XF2(I,L) C COMPUTE INTERFACIAL VELOCITIES DO 6 IN=2,INP1 SHI=SINH(ZN*DEPTHI(IN)) CHI=COSH(ZN*DEPTHI(IN)) C NUMBER OF THE COEFFICIENTS BELOW THE `IN' INTERFACE NCB=(NI-IN+1)*NN+(IN+1-NI)+N COEFI1=X(NCB)*CHI+X(NCB+1)*SHI COEFI2=X(NCB)*SHI+X(NCB+1)*CHI UI2(IN,I)=UI2(IN,I)-ZN*COEFI1*XF1(I,L) VI2(IN,I)=VI2(IN,I)-ZN*COEFI2*XF2(I,L) DU2DX(IN,I,NCB)=-ZN*CHI*XF1(I,L) DU2DX(IN,I,NCB+1)=-ZN*SHI*XF1(I,L) DU2DX(IN,I,1)=(ZN*COEFI1+ZN*ZN*DEPTHI(IN)*COEFI2)*XF1(I,L)/ 1 X(1)+DU2DX(IN,I,1) DV2DX(IN,I,1)=DV2DX(IN,I,1)+(ZN*ZN*DEPTHI(IN)*COEFI1+COEFI2*ZN) 1 *XF2(I,L)/X(1) DV2DX(IN,I,NCB)=-ZN*SHI*XF2(I,L) DV2DX(IN,I,NCB+1)=-ZN*CHI*XF2(I,L) DU2DZ(IN)=DU2DZ(IN)-ZN*ZN*COEFI2*XF1(I,L) DV2DZ(IN)=DV2DZ(IN)-ZN*ZN*COEFI1*XF2(I,L) C C COMPUTE UI1,VI1 FOR THE LOWER INTERFACE INM1=IN-1 SHI=SINH(ZN*DEPTHI(INM1)) CHI=COSH(ZN*DEPTHI(INM1)) COEFI1=X(NCB)*CHI+X(NCB+1)*SHI COEFI2=X(NCB)*SHI+X(NCB+1)*CHI UI1(INM1,I)=UI1(INM1,I)-ZN*COEFI1*XF1(I,L) VI1(INM1,I)=VI1(INM1,I)-ZN*COEFI2*XF2(I,L) DU1DX(INM1,I,NCB)=-ZN*CHI*XF1(I,L) DU1DX(INM1,I,NCB+1)=-ZN*SHI*XF1(I,L) DU1DX(INM1,I,1)=DU1DX(INM1,I,1)+(ZN*COEFI1+ZN*ZN*DEPTHI(INM1)* 1 COEFI2)*XF1(I,L)/X(1) DV1DX(INM1,I,1)=DV1DX(INM1,I,1)+(ZN*ZN*DEPTHI(INM1)*COEFI1+ 1 COEFI2*ZN)*XF2(I,L)/X(1) DV1DX(INM1,I,NCB)=-ZN*SHI*XF2(I,L) DV1DX(INM1,I,NCB+1)=-ZN*CHI*XF2(I,L) DU1DZ(INM1)=DU1DZ(INM1)-ZN*ZN*COEFI2*XF1(I,L) DV1DZ(INM1)=DV1DZ(INM1)-ZN*ZN*COEFI1*XF2(I,L) C COMPUTE PARAMETERS FOR THE `IN' INTERFACE DZETADX(INM1,I,1)=DZETADX(INM1,I,1)-ZN*DEPTHI(INM1)*COEFI1* 1 XF1(I,L)/X(1) DZETADX(INM1,I,NCB)=SHI*XF1(I,L) DZETADX(INM1,I,NCB+1)=CHI*XF1(I,L) 6 CONTINUE 14 CONTINUE C C FOR IN=1, COMPUTE VARIOUS VALUES DO 13 N=1,NLAST1 XN=N ZN=XN*CO DPT=DEPTHI(1)+D(0) SH1=SINH(ZN*DPT) CH1=COSH(ZN*DPT) NCB=NI*NN-NI+2+N UI2(1,I)=UI2(1,I)-ZN*X(NCB)*CH1*XF1(I,N) VI2(1,I)=VI2(1,I)-ZN*X(NCB)*SH1*XF2(I,N) DU2DX(1,I,1)=DU2DX(1,I,1)+ZN/X(1)*X(NCB)*XF1(I,N)* 1 (CH1+SH1*ZN*DPT) DV2DX(1,I,1)=DV2DX(1,I,1)+ZN/X(1)*X(NCB)*XF2(I,N)* 1 (SH1+CH1*ZN*DPT) DU2DX(1,I,NCB)=-ZN*CH1*XF1(I,N) DV2DX(1,I,NCB)=-ZN*SH1*XF2(I,N) DU2DZ(1)=DU2DZ(1)-ZN*ZN*X(NCB)*SH1*XF1(I,N) DV2DZ(1)=DV2DZ(1)-ZN*ZN*X(NCB)*CH1*XF2(I,N) 13 CONTINUE C COMPUTE SURFACE VELOCITY UETA=-C+UIN(NI)+USURF*(ETA+D(NI)) C FOR EACH INTERFACE COMPUTE THE TOP AND BOTTOM VELOCITY USING THE C STREAM FUNCTIONS FROM OVER (IN+1) AND UNDER THE INTERFACE (IN). DO 17 IN=1,INP1 IF(IN.EQ.INP1) GO TO 117 UZETA(IN)=UIN(IN)-C +U(IN+1)*(DEPTHI(IN)+D(IN))/(D(IN)-D(IN+1)) UI1(IN,I)=UZETA(IN)+UI1(IN,I) 117 CONTINUE UI2(IN,I)=UIN(IN-1)-C+U(IN)*(DEPTHI(IN)+D(IN-1))/(D(IN-1)-D(IN)) 1 +UI2(IN,I) 17 CONTINUE DO 15 N=1,NX DETADX(I,N)=DETADX(I,N)/UETA DO 15 IN=1,NI DZETADX(IN,I,N)=DZETADX(IN,I,N)/UZETA(IN) DU1DX(IN,I,N)=DU1DX(IN,I,N)+DU1DZ(IN)*DZETADX(IN,I,N) DV1DX(IN,I,N)=DV1DX(IN,I,N)+DV1DZ(IN)*DZETADX(IN,I,N) DU2DX(IN,I,N)=DU2DX(IN,I,N)+DU2DZ(IN)*DZETADX(IN,I,N) DV2DX(IN,I,N)=DV2DX(IN,I,N)+DV2DZ(IN)*DZETADX(IN,I,N) 15 CONTINUE UETA=UETA+UU C CHECK VELOCITIES C WRITE(*,*) UETA,VV,(UI1(L,I),VI1(L,I),UI2(L,I),VI2(L,I),L=1,INP1) DQDETA=1.0 DQDU=UETA/G DQDC=-DQDU DQDV=VV/G C BEGIN SETTING UP MATRIX EQUATION DO 18 N=1,NN 18 DQDX(I,N)=DQDETA*DETADX(I,N)+DQDU*(DUDX(N)+DUDETA*DETADX(I,N))+ 1 DQDV*(DVDX(N)+DVDETA*DETADX(I,N))+DQDC*DCDX(N) Q(I)=ETA+1.0/(2.0*G)*(UETA*UETA+VV*VV) 20 CONTINUE C QBAR=0. DO 41 I=2,NTHTM1,2 41 QBAR=QBAR+Q(I-1)+4.0*Q(I)+Q(I+1) QBAR=QBAR*DTHETL/3.0 WRITE(6,19) QBAR 19 FORMAT(/' BERNOULLI CONSTANT = ',F7.2,' FT.') WRITE(6,112) XH,TOTERR 112 FORMAT(' PRESENT WAVE HEIGHT=',F8.3,' FT.; ERROR= ',1PE11.4, 1' FT.') WRITE(6,113) X(1) 113 FORMAT(' PRESENT WAVE LENGTH=',F8.3,' FT.') DO 22 N=1,NN DQBAR(N)=0.0 DO 21 I=2,NTHTM1,2 21 DQBAR(N)=DQBAR(N)+DQDX(I-1,N)+4.*DQDX(I,N)+DQDX(I+1,N) 22 DQBAR(N)=DQBAR(N)*DTHETL/3. C C COMPUTE THE VARIOUS ERRORS ERRORV=0.0 ERRORU=0.0 ERRORQ=0.0 DO 24 I=1,NTHTS ERROR(I)=Q(I)-QBAR DO 24 IN=1,NI ERROR2(IN,I)=(UI2(IN,I)-UI1(IN,I))**2 ERROR3(IN,I)=(VI2(IN,I)-VI1(IN,I))**2 24 CONTINUE DO 123 I=2,NTHTM1,2 ERRORQ=ERRORQ+ ERROR(I-1)*ERROR(I-1)+4.*ERROR(I)*ERROR(I) 1 +ERROR(I+1)*ERROR(I+1) DO 123 IN=1,NI ERRORU=ERRORU+ERROR2(IN,I-1)+4.*ERROR2(IN,I)+ERROR2(IN,I+1) ERRORV=ERRORV+ERROR3(IN,I-1)+4.*ERROR3(IN,I)+ERROR3(IN,I+1) 123 CONTINUE ERRORV=SQRT(ERRORV*DTHETL/3.0) ERRORU=SQRT(ERRORU*DTHETL/3.0) ERRORQ=SQRT(ERRORQ*DTHETL/3.0) WRITE(6,224) ERRORQ,ERRORU,ERRORV 224 FORMAT(' RMS ERROR IN DFSBC = ',1PE9.3,' FT.'/' TOTAL RMS HORIZONT 1AL VELOCITY ERROR AT INTERFACES =',1PE9.3/' TOTAL RMS VERTICAL V 2ELOCITY ERROR AT INTERFACES =',1PE9.3) C WRITE(6,124) (ERROR(I),I=1,NTHTS) C 124 FORMAT(' THE DISTRIBUTED DFSBC ERRORS:'/(1P10E11.3)) NNP3=NN+3 C SET UP LOWER HALF OF MATRIX K=0 DO 30 N=1,NXML K=K+(N-1) DO 30 M=N,NX K=K+1 B(K)=0.0 IF(M.GT.NXML) GO TO 50 DO 26 I=2,NTHTM1,2 IF(M.GT.NN.OR.N.GT.NN) GO TO 131 B(K)=B(K)+(DQDX(I-1,N)-DQBAR(N))*(DQDX(I-1,M)-DQBAR(M)) 1 +4.*(DQDX(I,N)-DQBAR(N))*(DQDX(I,M)-DQBAR(M)) 2 +(DQDX(I+1,N)-DQBAR(N))*(DQDX(I+1,M)-DQBAR(M)) 131 DO 26 IN=1,NI B(K)=B(K)+ 1 (DU2DX(IN,I-1,N)-DU1DX(IN,I-1,N)) 2 *(DU2DX(IN,I-1,M)-DU1DX(IN,I-1,M)) 3 +(DV2DX(IN,I-1,N)-DV1DX(IN,I-1,N)) 4 *(DV2DX(IN,I-1,M)-DV1DX(IN,I-1,M)) 5 +4.*(DU2DX(IN,I,N)-DU1DX(IN,I,N))*(DU2DX(IN,I,M)-DU1DX(IN,I,M)) 6 +4.*(DV2DX(IN,I,N)-DV1DX(IN,I,N))*(DV2DX(IN,I,M)-DV1DX(IN,I,M)) 7 +(DU2DX(IN,I+1,N)-DU1DX(IN,I+1,N))*(DU2DX(IN,I+1,M)- 8 DU1DX(IN,I+1,M))+(DV2DX(IN,I+1,N)-DV1DX(IN,I+1,N)) 9 *(DV2DX(IN,I+1,M)-DV1DX(IN,I+1,M)) 26 CONTINUE B(K)=B(K)*DTHETL/3.0 GO TO 30 50 IF(M.GT.(NXMLP1)) GO TO 59 DO 51 I=2,NTHTM1,2 51 B(K)=B(K)+DETADX(I-1,N)+4.0*DETADX(I,N)+DETADX(I+1,N) B(K)=B(K)*DTHETL/3.0 GO TO 30 59 IF(M.EQ.NX) GO TO 69 IN=M-(NXMLP1) DO 60 I=2,NTHTM1,2 60 B(K)=B(K)+DZETADX(IN,I-1,N)+4.0*DZETADX(IN,I,N)+DZETADX(IN,I+1,N) B(K)=B(K)*DTHETL/3.0 GO TO 30 69 B(K)=DETADX(1,N)-DETADX(NTHTS,N) 30 CONTINUE C USING SYMMETRY, SET UP THE UPPER HALF OF THE MATRIX DO 140 N=2,NX NM1=N-1 K=NX*NM1 DO 140 M=1,NM1 K=K+1 B(K)=B(N+(M-1)*NX) 140 CONTINUE C SET UP RIGHT HAND SIDE, THE A VECTOR DO 58 M=1,NX-1 A(M)=0.0 IF(M.GT.NXML) GO TO 53 DO 57 I=2,NTHTM1,2 IF(M.GT.NN) GO TO 233 A(M)=A(M)+(QBAR-Q(I-1))*(DQDX(I-1,M)-DQBAR(M)) 1 +4.*(QBAR-Q(I))*(DQDX(I,M)-DQBAR(M)) 2 +(QBAR-Q(I+1))*(DQDX(I+1,M)-DQBAR(M)) C A(M)=A(M)+(QBAR-Q(I-1))*DQDX(I-1,M) C 1 +4.*(QBAR-Q(I))*DQDX(I,M) C 2 +(QBAR-Q(I+1))*DQDX(I+1,M) 233 DO 57 IN=1,NI A(M)=A(M)+ 1 (UI1(IN,I-1)-UI2(IN,I-1))*(DU2DX(IN,I-1,M)-DU1DX(IN,I-1,M)) 2 +(VI1(IN,I-1)-VI2(IN,I-1))*(DV2DX(IN,I-1,M)-DV1DX(IN,I-1,M)) 3 +4.*(UI1(IN,I)-UI2(IN,I))*(DU2DX(IN,I,M)-DU1DX(IN,I,M)) 4 +4.*(VI1(IN,I)-VI2(IN,I))*(DV2DX(IN,I,M)-DV1DX(IN,I,M)) 5 +(UI1(IN,I+1)-UI2(IN,I+1))*(DU2DX(IN,I+1,M)-DU1DX(IN,I+1,M)) 6 +(VI1(IN,I+1)-VI2(IN,I+1))*(DV2DX(IN,I+1,M)-DV1DX(IN,I+1,M)) 57 CONTINUE GO TO 58 53 IF(M.GT.(NXMLP1)) GO TO 159 DO 61 I=2,NTHTM1,2 A(M)=A(M)+ZETA(INP1,I-1)+4.0*ZETA(INP1,I)+ZETA(INP1,I+1) 61 CONTINUE GO TO 58 159 IN=M-(NXMLP1) DO 62 I=2,NTHTM1,2 62 A(M)=A(M)+ZETA(IN,I-1)+4.0*ZETA(IN,I)+ZETA(IN,I+1) 58 A(M)=A(M)*DTHETL/3.0 C 66 A(NX)=HT-(ZETA(INP1,1)-ZETA(INP1,NTHTS)) 75 WRITE(6,96) A(NXMLP1) 96 FORMAT(' MSL=',1PE12.4) A(NXMLP1)=-A(NXMLP1) WRITE(6,102) 102 FORMAT(' MEAN INTERFACE DISPLACEMENTS:') ERRORI=0.0 DO 101 IN=1,IN NCB=NXML+IN+1 WRITE(6,103) IN,A(NCB) A(NCB)=-A(NCB) ERRORI=ERRORI+ABS(A(NCB)) 103 FORMAT(I3,1PE12.4) 101 CONTINUE IF(IFSET.EQ.1) RETURN TOTERR=ERRORQ+ERRORV+ERRORU+ERRORI+ABS(A(NXMLP1))+TOTERR C WRITE(*,*) TOTERR NX2=NX*NX CALL SIMQ(B,A,NX,KS) WRITE(6,1191) KI 1191 FORMAT(///' ITERATION',I3) IF(KS.EQ.1) WRITE(6,87) 87 FORMAT(' CFF MATRIX IS SINGULAR') DO 40 N=1,NXML 40 X(N)=X(N)+A(N)*DAMP WRITE(6,91) (X(N),N=1,NXML) 91 FORMAT(' X(N)=',1P6E12.4) RETURN END C*********************************************************************** SUBROUTINE SIMQ(A,B,N,KS) PARAMETER(IX=104,IC=62) DIMENSION A(IX*IX),B(IX) C C FORWARD SOLUTION C TOL=0.0 KS=0 JJ=-N DO 65 J=1,N JY=J+1 JJ=JJ+N+1 BIGA=0 IT=JJ-J DO 30 I=J,N C C SEARCH FOR MAXIMUM COEFFICIENT IN COLUMN C IJ=IT+I IF(ABS(BIGA)-ABS(A(IJ))) 20,30,30 20 BIGA=A(IJ) IMAX=I 30 CONTINUE C C TEST FOR PIVOT LESS THAN TOLERANCE (SINGULAR MATRIX) C IF(ABS(BIGA)-TOL) 35,35,40 35 KS=1 RETURN C C INTERCHANGE ROWS IF NECESSARY C 40 I1=J+N*(J-2) IT=IMAX-J DO 50 K=J,N I1=I1+N I2=I1+IT SAVE=A(I1) A(I1)=A(I2) A(I2)=SAVE C C DIVIDE EQUATION BY LEADING COEFFICIENT C 50 A(I1)=A(I1)/BIGA SAVE=B(IMAX) B(IMAX)=B(J) B(J)=SAVE/BIGA C C ELIMINATE NEXT VARIABLE C IF(J-N) 55,70,55 55 IQS=N*(J-1) DO 65 IIX=JY,N IXJ=IQS+IIX IT=J-IIX DO 60 JX=JY,N IXJX=N*(JX-1)+IIX JJX=IXJX+IT 60 A(IXJX)=A(IXJX)-(A(IXJ)*A(JJX)) 65 B(IIX)=B(IIX)-(B(J)*A(IXJ)) C C BACK SOLUTION C 70 NY=N-1 IT=N*N DO 80 J=1,NY IA=IT-J IB=N-J IIC=N DO 80 K=1,J B(IB)=B(IB)-A(IA)*B(IIC) IA=IA-N 80 IIC=IIC-1 RETURN END C*********************************************************************** SUBROUTINE WVLEN(DPT,PER,XLENTH) C THIS CALCULATES LINEAR WAVE LENGTH BY NEWTON'S METHOD XKHO=(6.2831853/PER)**2.*DPT/32.2 IF(XKHO-6.3) 2,1,1 1 XKH=XKHO GO TO 9 2 XKH=SQRT(XKHO) 3 SH=SINH(XKH) CH=COSH(XKH) EPS=XKHO-XKH*SH/CH SLOPE=-XKH/CH**2.-SH/CH DXKH=-EPS/SLOPE IF(ABS(DXKH/XKH)-0.0001) 9,9,4 4 XKH=XKH+DXKH GO TO 3 9 XLENTH=6.283185*DPT/XKH RETURN END C*********************************************************************** SUBROUTINE VELOC C C COMPUTES THE VELOCITY PROFILE FROM TOP TO BOTTOM FOR A GIVEN PHASE C POSITION C PARAMETER(IX=104,IC=62,NNI=3) DIMENSION X(IX),XF1(IC,IX),XF2(IC,IX),U(0:NNI+1),UIN(0:NNI+1) DIMENSION THETA(IC),ZETA(NNI+1,IC),D(0:NNI+1) COMMON X,XF1,XF2 COMMON /DATA/NN,NTHTS,QBAR,T,D,H,U,UIN,THETA,DTHETA,NI, 1ZETA,HT PI=3.1415927 NX=NN+(NN-2)/2+4 NXM3=NX-3 NNP1=NN+1 NNP2=NN+2 NTHTM1=NTHTS-1 C=X(1)/T CO=2.*PI/X(1) CON=PI/X(1) 450 WRITE(*,*) ' TO DETERMINE THE VELOCITIES AND ACCELERATIONS:' 500 WRITE(*,*) ' INPUT X POSITION AND NUMBER OF VERTICAL SEGMENTS:' 606 WRITE(*,*) ' Use 999 for X to guit.' READ(*,*) XX,NDY IF(XX.EQ.999) GO TO 1166 PX=XX*2.*PI/X(1) WRITE(*,*) 'X POSITION AND PHASE:',XX,PX I=1 DO 6 N=1,(NN/2-1) XF1(I,N)=COS(PX*FLOAT(N)) 6 XF2(I,N)=SIN(PX*FLOAT(N)) C DETERMINE THE SURFACE AND INTERFACIAL DISPLACEMENTS AT XX C DO 5 IN=1,NI NINTER=(NI-IN+1)*NN+IN-NI+1 CALL SURFAS(IN,NINTER,2) 5 CONTINUE CALL SURFAS(NI+1,2,2) DY=(ZETA(NI+1,1)+D(0))/(NDY-1) WRITE(6,74) WRITE(6,90) 74 FORMAT(//' VELOCITIES AND ACCELERATIONS ') 90 FORMAT( /' ELEVATION 1 U V DUDT DVDT PSI') IN=NI+1 Y=ZETA(IN,I)+DY KOUNT=1 1 Y=Y-DY INP=IN DO 555 INN=1,INP IF(Y.LT.(-D(IN-1)+ZETA(IN-1,I))) THEN IN=IN-1 ELSE GO TO 556 ENDIF 555 CONTINUE 556 CONTINUE IF (IN.EQ.1) GO TO 50 C COMPUTE THE UPPER LAYER VELOCITIES UU=0.0 V=0.0 PSI=0.0 DUDX=0.0 DEN=D(IN-1)-D(IN) DUDZ=U(IN)/DEN NCB=(NI-IN+1)*NN+(IN+1-NI) DO 14 N=1,(NN-2),2 L=(N+1)/2 ZN=FLOAT(N+1)*CON SH=SINH(ZN*Y) CH=COSH(ZN*Y) COEF1=X(NCB+N)*CH+X(NCB+N+1)*SH COEF2=X(NCB+N)*SH+X(NCB+N+1)*CH UU=UU-ZN*COEF1*XF1(I,L) V=V-ZN*COEF2*XF2(I,L) PSI=PSI+COEF2*XF1(I,L) DUDX=DUDX+ZN*ZN*COEF1*XF2(I,L) 14 DUDZ=DUDZ-ZN*ZN*COEF2*XF1(I,L) UU=UIN(IN-1)+U(IN)*(Y+D(IN-1))/DEN +UU PSI=PSI-(UIN(IN)-C)*Y-U(IN)*(D(IN-1)*Y+Y*Y/2.0)/DEN DUDT=(UU-C)*DUDX+V*DUDZ DVDT=(UU-C)*(DUDZ-U(IN)/DEN)-V*DUDX WRITE(6,100) Y,UU,V,DUDT,DVDT,PSI C FIND THE ELEVATION IN NEAREST MULTIPLE OF DY C IF (KOUNT.EQ.1) THEN C Y=FLOAT(INT(ZETA(NI+1,I)/DY)+1)*DY C IF (ZETA(NI+1,I).LT.0.0) Y=Y-DY C ELSE C ENDIF KOUNT=KOUNT+1 100 FORMAT(6F12.2) GO TO 1 60 Y=Y-DY IF(Y.LT.-D(0)) GO TO 500 C COMPUTE KINEMATICS IN THE BOTTOM LAYER 50 UU=0.0 V=0.0 PSI=0.0 DUDX=0.0 DPT=D(0) DEN=DPT-D(1) DUDZ=U(1)/DEN NLAST=(NN-2)/2 DO 15 N=1,NLAST XN=N ZN=CO*XN ELEV=Y+DPT NCB=NI*NN-NI+2+N SH=SINH(ZN*ELEV) CH=COSH(ZN*ELEV) COEF1=X(NCB)*CH COEF2=X(NCB)*SH UU=UU-ZN*COEF1*XF1(I,N) V=V-ZN*COEF2*XF2(I,N) PSI=PSI+COEF2*XF1(I,N) DUDX=DUDX+ZN*ZN*COEF1*XF2(I,N) 15 DUDZ=DUDZ-ZN*ZN*COEF2*XF1(I,N) UU=UU+UIN(0)+U(1)*(Y+DPT)/DEN PSI=PSI-(UIN(1)-C)*Y-U(1)*Y*Y/(2.0*DEN) DUDT=(UU-C)*DUDX+V*(DUDZ) DVDT=(UU-C)*(DUDZ-U(1)/DEN)-V*DUDX WRITE(6,100) Y,UU,V,DUDT,DVDT,PSI GO TO 60 1166 RETURN END