[Fortran] 纯文本查看 复制代码
PROGRAM BOUSSINESQ
REAL,ALLOCATABLE,DIMENSION(:)::YITA1,YITA2,XU1,XU2,B1,B2,YITAF,XUF,YITAT,XUT,XUD,YITAD
!YITA为波高,XU为速度,B为方程组的向量,YITAF,XUF为上一步的初值
INTEGER::N1=100,K=1000
!定义变量部分
REAL::XL,DETX,DETT,XK,OMEGA,U0,AA,XD,CK,CC,T,XX,IMAX,TOL,DX2
!定义波长XL,观察点的距离间隔DETX和时间间隔DETT,波数XK,角速度OMEGA,初速度U0,振幅AA,水深XD,CK=DETT/DETX,
!CC为Sommerfled边界里的系数,T为周期,XX为观察点,IMAX为迭代最大次数,TOL为迭代判断精读
REAL::PI,GRA
PARAMETER(PI=3.1415926,GRA=9.81)
!GRA为重力加速度
PARAMETER(T=8.0,XD=1.5,AA=0.1,TOL=1E-5,IMAX=200)
DETT=T/20.0
CC=(GRA*XD)**(1.0/2.0)
OMEGA=T/(2.0*PI)
XK=((OMGRA*OMGRA)/(GRA*XD))**(1.0/2.0)
XL=2.0*PI/XK
DETX=XL/10.0
XX=100
!边界条件的设定
!时间为0时的初始条件
DO I=1,N1
XU1(I)=0.0
YITA1(I)=0.0
ENDDO
!差分计算部分
DO J=1,K
!空间为0的边界条件
XU1(1)=(PI*2.0*AA*COSH(XK*XD))*COS(OMEGA*(J-1)*DETT)/(T*SINH(XK*XD))
YITA1(1)=AA*COS(OMEGA*(J-1)*DETT)
DO I=1,N1
XUF(I)=XU1(I)
YITAF(I)=YITA1(I)
ENDDO
!预测步
!切记！！数字写成小数形式,不要写成整数（如1要写成1.0）
DO I=1,N1
SELECT CASE(I)
CASE(1)
YITA2(I)=(PI*2.0*AA*COSH(XK*XD))*COS(OMEGA*(J-(1.0))*DETT)/(T*SINH(XK*XD))
XU2(I)=AA*COS(OMEGA*(J-(1.0))*DETT)
CASE(2:99)
YITA2(I)=YITA1(I)-CK*XU1(I)*(YITA1(I+1)-YITA1(I))-CK*(XD+YITA1(I))*(XU1(I+1)-XU1(I))
XU2(I)=XU1(I)-CK*XU1(I)*(XU1(I+1)-XU1(I))-CK*GRA*(YITA1(I+1)-YITA1(I))
CASE(100)
YITA2(I)=YITA1(I)-CK*CC*(YITA1(I)-YITA1(I-1))
XU2(I)=XU1(I)-CK*CC*(XU1(I)-XU1(I-1))
END SELECT
ENDDO
DO I=1,N1
XUD(I)=XU1(I)
YITAD(I)=YITA1(I)
ENDDO
DO I=1,N1
SELECT CASE(I)
CASE(1)
B1(I)=YITAF(I)-CK*XUF(I)*(YITAF(I+1)-YITAF(I))-CK*(XD+YITAF(I))*(XUF(I+1)-XUF(I))-CK*(XD+YITA1(I))*(XU1(I+1)-XU1(I))
CASE(2:99)
B1(I)=YITAF(I)-CK*XUF(I)*(YITAF(I+1)-YITAF(I))-CK*(XD+YITAF(I))*(XUF(I+1)-XUF(I))-CK*(XD+YITA1(I))*(XU1(I+1)-XU1(I))
CASE(100)
B1(I)=YITAF(I)-CK*CC*(YITAF(I)-YITAF(I-1))
END SELECT
ENDDO
DO I=1,N1
YITA1(I)=B1(I)
ENDDO
DO I=1,N1
SELECT CASE(I)
CASE(1)
B2(I)=XUF(I)-CK*XUF(I)*(XUF(I+1)-XUF(I))-CK*GRA*(YITAF(I+1)-YITAF(I))-CK*GRA*(YITA1(I+1)-YITA1(I))
CASE(2:99)
B2(I)=XUF(I)-CK*XUF(I)*(XUF(I+1)-XUF(I))-CK*GRA*(YITAF(I+1)-YITAF(I))-CK*GRA*(YITA1(I+1)-YITA1(I))
CASE(100)
B2(I)=XUF(I)-CK*CC*(XUF(I)-XUF(I-1))
END SELECT
ENDDO
DO I=1,N1
XU1(I)=B2(I)
ENDDO
DX2=0.0
DO I=1,N1
IF(ABS(XU1(I)-XUD(I)).GT.DX2)THEN
DX2=ABS(XU1(I)-XUD(I))
ENDIF
IF(ABS(YITA1(I)-YITAD(I)).GT.DX2)THEN
DX2=ABS(YITA1(I)-YITAD(I))
ENDIF
ENDDO
IF(DX2.LT.TOL)EXIT
YITAT(J)=YITA1(XX)
XUT(J)=XU1(XX)
ENDDO
!结果输出部分
DO I=1,N1
WRITE(*,*)YITA1(I)
ENDDO
DO I=1,N1
WRITE(*,*)XU1(I)
ENDDO
END PROGRAM BOUSSINESQ