戻る

C ********************************************************************* C
C C
C F2D - FDM ( uvp) C
C C
C By Fourth ordered FDM C
C C
C Equal length & Unequal length of mesh sizes C
C C
C ( V. 4.0 ) C
C C
C Copyright : Yasuhiro MATSUDA C
C C
C ********************************************************************* C
C
PROGRAM MAIN
C
IMPLICIT REAL*8 ( A-H,O-Z)
C
C ----- DATA ( 1 ) -----------------------------------------
C
*** PARAMETER ( NX= 41, NY= 41 )
C
PARAMETER ( NX= 21, NY= 21) ! Re= 100 UE: DT = 0.014 EM: DT=0.04
C
C ----------------------------------------------------------
C
PARAMETER( CAV= 1.D0, LX= NX-1, LY= NY-1, NNP= NX*NY, NEL= LX*LY)
C
DIMENSION U ( NX,NY), UA( NX,NY), UB( NX,NY), OU( NX,NY),
1 V ( NX,NY), VA( NX,NY), VB( NX,NY), OV( NX,NY),
2 P ( NX,NY), PA( NX,NY), B ( NX,NY),
C
3 UN( NNP), VN( NNP), PN( NNP), DP( NX,NY),
4 X ( NNP), Y ( NNP), HX( LX), HY( LY),
5 XZ( NX), YZ( NY) , XL( LX), YL( LY),
6 NDE( NEL,4)
C
DIMENSION AX( NX,NY), BX( NX,NY), FU ( NX+3), FV( NX+3),
1 DX( NX), DY( NY), SSP( NX,NY)
C
C ----------------------------------------------------------
C
COMMON /VRD/ AP
COMMON /SIZ/ II, JJ, IIP1, JJP1, IIM1, JJM1
COMMON /OUT/ DT, RE, WEI, NPR, LOPE, MSH
COMMON /TES/ A1, A2, A3, A4, A5, A6, A7, A8
C
OPEN ( 8, FILE='EPS.DAT' )
C
C ----- DATA ( 2 ) -----------------------------------------
C
*** OPEN (11, FILE='EM20_F2D_FDM_R100.DAT' ) ! Re=100 DT=0.04
C
OPEN (11, FILE='UE20_F2D_FDM_R100.DAT' ) ! Re=100 DT=0.014
C
* OPEN (11, FILE='UE40_F2D_FDM_R400_009.DAT' ) ! Re=400 DT=0.009
* OPEN (11, FILE='UE40_F2D_FDM_R400_008.DAT' ) ! Re=400 DT=0.008
C
* OPEN (11, FILE='UE40_F2D_FDM_R400_007.DAT' ) ! Re=400 DT=0.007
C
*** OPEN (11, FILE='UE40_F2D_FDM_R1000.DAT' )
C
C ----------------------------------------------------------
C
CALL GMSH ( XL, YL, HX, HY, LX, LY)
C
CALL STND ( X, Y, XZ, YZ, NX, NY, LX, LY, NNP, NEL,
1 CAV )
C
CALL SBIN ( U, UA, UB, V, VA, VB, B, OU, OV, P,
1 PA, DP, HX, HY, XZ, YZ, NX, NY, LX )
C
LOP = 0
C
C --------------------------------------------------------------------
9999 LOP = LOP + 1
C
CALL FVLB ( U, V, UB, VB, HX, HY, NX, NY)
C
CALL SORP ( UB, VB, PA, HX, HY, NX, NY, IT1, B, DX, DY,
1 SSP, LOP )
C
CALL VLSX ( U, V, UA, VA, HX, NX, NY, AX, BX, FU, FV)
C
CALL VLSY ( U, V, UA, VA, HY, NX, NY, AY, BY, FU, FV)
C
CALL VLSF ( U, V, UA, VA, PA, HX, HY, NX, NY)
C
CALL SORP ( UA, VA, DP, HX, HY, NX, NY, IT2, B, DX, DY,
1 SSP, LOP )
C
CALL PRVL ( U, V, UA, VA, P, PA, DP, OU, OV, HX, HY,
1 NX, NY )
C
CALL CHCK ( U, V, OU, OV, P, UN, VN, PN, NX, NY,
1 NNP, LOP, IT1, IT2, HX, *9999 )
C
CALL PLTO ( X, Y, NDE, P, PN, U, V, UN, VN, XZ, YZ,
1 HX, HY, NX, NY, NNP, NEL, LX, LY, LOP )
C
CLOSE ( 8)
CLOSE (11)
C
STOP
END
C **********************************************************************
C
SUBROUTINE PLTO ( X, Y, NDE, P, PN, U, V, UN, VN,
1 XZ, YZ, HX, HY, IX, IY, NNP, NEL, LX,
2 LY, LOP )
C
IMPLICIT REAL*8 ( A-H,O-Z)
C
COMMON /DEF/ VUV(10000), VRP(10000)
COMMON /OUT/ DT, RE, WEI, NPR, LOPE, MSH
C
DIMENSION X ( NNP), Y ( NNP), NDE( NEL,4), U ( IX,IY),
1 V ( IX,IY), UN( NNP), VN ( NNP), P ( IX,IY),
2 PN( NNP), XZ( IX), YZ ( IY), HX( LX),
3 HY( LY)
C
C ----------------------------------------------------------
C
CALL CNVT ( U, V, P, UN, VN, PN, IX, IY, NNP )
C
C ----------------------------------------------------------
DO 100 J = 1, IY-1
DO 100 I = 1, IX-1
C ----------------------------------------------------------
NUM = ( J-1)* ( IX-1) + I
NDE( NUM,1) = ( J-1)*( IX) + I
NDE( NUM,2) = ( J-1)*( IX) + I + 1
NDE( NUM,3) = ( J-1)*( IX) + I + ( IX + 1)
NDE( NUM,4) = ( J-1)*( IX) + I + ( IX)
C
100 CONTINUE
C ----------------------------------------------------------
C
WRITE(11,2000) RE, DT, LOPE, MSH, LOP
WRITE(8, 2000) RE, DT, LOPE, MSH, LOP
2000 FORMAT(F8.2,F7.5,I5,I2,I5)
C
WRITE(11,2200) ( PN ( I), I = 1, NNP)
WRITE(11,2200) ( VUV( I), I = 1, LOP)
WRITE(11,2200) ( VRP( I), I = 1, LOP)
WRITE(11,2200) ( X( I), I = 1, NNP)
WRITE(11,2200) ( Y( I), I = 1, NNP)
C
WRITE(11,2200) (( U( I,J), J = 1, IY), I = 1, IX )
WRITE(11,2200) (( V( I,J), J = 1, IY), I = 1, IX )
WRITE(11,2200) ( HX( I), I = 1, IX-1 )
WRITE(11,2200) ( HY( I), I = 1, IY-1 )
2200 FORMAT(10F8.3)
C
WRITE(11,2400) (( NDE( I,J), J = 1, 4), I = 1, NEL )
2400 FORMAT(20I4)
C
WRITE(6,*) ' '
WRITE(6,*) ' *** End ***'
WRITE(8,*) ' '
WRITE(8,*) ' *** End ***'
C ----------------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE CHCK ( U, V, OU, OV, P, UN, VN, PN, NX, NY,
1 NNP, LOP, IT1, IT2, HX, * )
C
IMPLICIT REAL*8 ( A-H,O-Z)
C
DIMENSION U ( NX,NY), V ( NX,NY), OU( NX,NY), OV( NX,NY),
1 P ( NX,NY), UN( NNP), VN( NNP), PN( NNP),
2 HX( NX-1)
C
COMMON /DEF/ VUV(10000), VRP(10000)
COMMON /OUT/ DT, RE, WEI, NPR, LOPE, MSH
C
PMN = 1000.D0
PMX = - 1000.D0
UMX = - 1000.D0
C --------------------------------------
DO 100 J = 1, NY
DO 100 I = 1, NX
C
PMN = DMIN1( PMN, P( I,J))
PMX = DMAX1( PMX, P( I,J))
UMX = DMAX1( UMX, U( I,J))
C
100 CONTINUE
C --------------------------------------
EPV = 0.D0
C ----------------------------------------------------------
DO 200 J = 2, NY-1
DO 200 I = 2, NX-1
C
IF ( DABS( U( I,J)).GT. UMX* 0.01D0)
1 EPV = DMAX1( EPV, DABS(( U(I,J) - OU(I,J))/ U(I,J)))
C
IF ( DABS( V( I,J)).GT. UMX*0.01D0)
1 EPV = DMAX1( EPV, DABS(( V(I,J) - OV(I,J))/ V(I,J)))
C
200 CONTINUE
C ----------------------------------------------------------
C
VUV( LOP) = EPV* 100.D0 ! (%)
C ============================================
WRITE(8,*) EPV* 100.D0, LOP
C
C ----------------------------
IF ( MOD( LOP,NPR).EQ.0) THEN ! PLOT/LOP
WRITE(6,2000) LOP, IT1, IT2, VUV( LOP), PMN, PMX
2000 FORMAT(' * LP=',I5,' IT1=',I4,' IT2=',I4,' Vr_UV=',D11.4,
1 ' Pmi=',F6.2,' Pmx=',F6.2)
END IF
C ----------------------------
C
C ----- Stationarity Check etc.---
C
IF ( UMX.GT.1.5D0) GO TO 5200
IF ( LOP.EQ.LOPE) GO TO 5700
IF ( VUV( LOP).GT.0.1D0) RETURN 1 ! 0.1(%)
C
** IF ( EPV.GT.0.001D0) RETURN 1 ! 0.1(%)
C =======================================================
GO TO 5600 ! Converged.
C
5200 WRITE(6,2200) LOP
2200 FORMAT(/,' * Diverged at Loop =', I6,' U-max. > 1.5 * STOP *')
STOP
C
5600 CONTINUE
C
WRITE(6,2300) LOP
2300 FORMAT(/,' * Stationary condition reached at Loop =',I4,/)
C
C ------------------------------------------------
C
BMX = 0.D0
DO 300 I = 1, NX-1
C
BVAR = UMX* DT/ HX( I) ! Max. bx.-value
IF ( BVAR.GE.BMX) BMX = BVAR
C
300 CONTINUE
C ------------------------------------------------
C
WRITE(6,2500) LOP, IT1, IT2, VUV( LOP), PMN, PMX, BMX
2500 FORMAT(' * LP=',I5,' IT1=',I4,' IT2=',I4,' Vr_UV=',D11.4,
1 /,11X,' Pmin=',F6.2,' Pmax=',F6.2,' Max.bx=',F7.3)
C -----------------------------------------------------------------
C
RETURN
C
5700 CONTINUE
WRITE(6,2400) LOP
2400 FORMAT(/,' * Simulation end : Loop =',I6)
WRITE(6,2500) LOP, IT1, IT2, VUV( LOP), PMN, PMX, BMX
C
RETURN
END
C **********************************************************************
C
SUBROUTINE SBIN ( U, UA, UB, V, VA, VB, B, OU, OV, P,
1 PA, DP, HX, HY, XZ, YZ, NX, NY, LX )
C
IMPLICIT REAL*8 ( A-H,O-Z)
C

DIMENSION U ( NX,NY), UA( NX,NY), UB( NX,NY), V ( NX,NY),
1 VA( NX,NY), VB( NX,NY), B ( NX,NY), OU( NX,NY),
2 OV( NX,NY), P ( NX,NY), PA( NX,NY), DP( NX,NY),
3 HX( NX), HY( NY), XZ( NX), YZ( NY)
C
C ----- Data Input ---------------------
C
CALL SINP
C
C --------------------------------------
DO 100 I = 1, LX
C
HX( I) = XZ( I+1) - XZ( I)
HY( I) = YZ( I+1) - YZ( I)
C
100 CONTINUE
C --------------------------------------
HX( NX) = HX( NX-1)
HY( NY) = HY( NY-1)
C
C --------------------------------------
DO 200 J = 1, NY
DO 200 I = 1, NX
C --------------------------------------
C
U ( I,J) = 0.D0
UA( I,J) = 0.D0
UB( I,J) = 0.D0
V ( I,J) = 0.D0
VA( I,J) = 0.D0
VB( I,J) = 0.D0
C
B ( I,J) = 0.D0
OU( I,J) = 0.D0
OV( I,J) = 0.D0
C
C ----- Attention ------
C
P ( I,J) = 0.D0
C ======================
PA( I,J) = 0.D0
DP( I,J) = 0.D0
C
200 CONTINUE
C --------------------------------------
DO 300 I = 1, NX
C
U ( I,NY) = 1.D0
UA( I,NY) = 1.D0
UB( I,NY) = 1.D0
C
300 CONTINUE
C --------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE GMSH ( XL, YL, HX, HY, LX, LY)
C
IMPLICIT REAL*8 ( A-H,O-Z)
C
COMMON /OUT/ DT, RE, WEI, NPR, LOPE, MSH
C
DIMENSION XL( LX), YL ( LY), HX( LX), HY( LY)
C
C ----------------------------------------------------------
C
CALL TTLE
C
C ----------------------------------------------------------
C
C MSH = 1 : Meshes of equal length ----- EM
C
C MSH = 2 : JSME-Paper Version --------- UEM ( UM)
C
C ----------------------------------------------------------
WRITE(6,*) ' '
WRITE(6,*) ' ** MSH = 1 ( Meshes of equal length ) or 2 ( Meshes o
1f unequal length ) ?'
*** READ(5,*) MSH
C
MSH = 2
WRITE(6,2000) MSH
2000 FORMAT(' ** MSH =',I2,/)
C
C --------------------------------------
IF ( MSH.EQ.1) GO TO 150
IF ( MSH.EQ.2) GO TO 170
C --------------------------------------
150 CONTINUE
C
WRITE(6,*) '* Meshes of equal length '
C ------------------------------------------------
DO 100 I = 1, LX
C
XL( I) = DFLOAT( I-1)/ DFLOAT( LX-1) ! Attention
100 CONTINUE
C ------------------------------------------------
DO 200 J = 1, LY
C
YL( J) = DFLOAT( J-1)/ DFLOAT( LY-1) ! Attention
200 CONTINUE
C ---------------------------------------------
GO TO 900
C
170 CONTINUE
C ----------------------------------------------------------
C
WRITE(6,*) '*** Meshes of unequal length ( UEM ) ***'
C
C ----------------------------------------------------------
C
C ----- JSME - Paper Ref. ---
C
AA = 2.D0
DX = 2.D0/ DFLOAT( LX-1)
DY = 2.D0/ DFLOAT( LY-1)
C
C ----- Attention --- Important ---
C
MDX = IDNINT( DFLOAT( LX+1)/ 2.D0)
MX1 = MDX - 1
MDY = IDNINT( DFLOAT( LY+1)/ 2.D0)
MY1 = MDY - 1
C
XL( MDX) = 0.5D0
YL( MDY) = 0.5D0
C
C ----------------------------------------------------
DO 300 I = 1, MX1
C
XL( I) = ( DEXP( AA* DFLOAT( I-1)* DX) - 1.D0)
1 /( 2.D0* ( DEXP( AA) - 1.D0))
C
XL( LX+1-I) = 1.D0 - XL ( I)
C
300 CONTINUE
C ----------------------------------------------------
DO 400 J = 1, MY1
C
YL( J) = ( DEXP( AA* DFLOAT( J-1)* DY) - 1.D0)
1 /( 2.D0* ( DEXP( AA ) - 1.D0))
C
YL( LY+1-J) = 1.D0 - YL( J)
C
400 CONTINUE
C -----------------------------------------------------------
C
900 CONTINUE
C ------------------------------------------------
DO 500 I = 1, LX-1
C
HX( I) = DABS( XL( I+1) - XL( I))
500 CONTINUE
C ------------------------------------------------
DO 600 J = 1, LY-1
C
HY( J) = DABS( YL( J+1) - YL( J))
600 CONTINUE
C ------------------------------------------------
C
RETURN
END
C ***********************************************************************
C
BLOCK DATA
C
IMPLICIT REAL*8 ( A-H,O-Z)
C
COMMON /TES/ A1, A2, A3, A4, A5, A6, A7, A8
C
C ----- Attention : For Meshes of Equal length ---
C
DATA A1/ 0.583333333333333D0/, ! = 7/12
1 A2/ 0.625D0/, ! = 5/8
2 A3/ -8.33333333333333D-2/, ! = -1/12
3 A4/ -0.125D0/, ! = -1/8
C
4 A5/ -8.33333333333333D-2/, ! = -1/12
5 A6/ -4.16666666666667D-2/, ! = -1/24
6 A7/ 8.33333333333333D-2/, ! = 1/12
7 A8/ 4.16666666666667D-2/ ! = 1/24
C
END
C ***********************************************************************
C
SUBROUTINE STND ( X, Y, XZ, YZ, NX, NY, LX, LY, NNP,
1 NEL, CAV )
C
IMPLICIT REAL*8 ( A-H,O-Z)
C
COMMON /OUT/ DT, RE, WEI, NPR, LOPE, MSH
C
DIMENSION X( NNP), Y( NNP), XZ( NX), YZ( NY)
C
IF ( MSH.NE.1) GO TO 120
C
C ----- In case of Equal length of Meshes ---
C
DO 100 I = 1, NX
C
XZ( I) = CAV/ DFLOAT( LX)* DFLOAT( I-1)
YZ( I) = XZ( I)
C
100 CONTINUE
C -------------------------------------------
GO TO 333
C
C ----- In case of Unequal length of Meshes ---
C
120 MNX = ( NX+1)/ 2
MNY = ( NY+1)/ 2
C ---------------------------------
DO 200 I = 1, MNX - 1
C
XZ( I) = ( DEXP( 2.D0* DFLOAT( I-1)* 2.D0/ DFLOAT( LX)) - 1.D0)
1 /( 2.D0* ( DEXP( 2.D0) - 1.D0))
C
200 CONTINUE
XZ( MNX) = 0.5D0
C ---------------------------------
DO 300 I = MNX + 1, LX
C
XZ( I) = 1.D0 - XZ( NX+1-I)
300 CONTINUE
XZ( NX) = 1.D0
C ---------------------------------
DO 400 I = 1, NX
C
YZ( I) = XZ( I)
400 CONTINUE
C --------------------------------
C
333 J = 1
C ---------------------------------
DO 500 IY = 1, NY
DO 500 IX = 1, NX
C
X( J) = XZ( IX)
Y( J) = YZ( IY)
J = J + 1
C
500 CONTINUE
C ---------------------------------
C
RETURN
END
C ***********************************************************************
C
SUBROUTINE CNVT ( U, V, P, UN, VN, PN, NX, NY, NNP )
C
IMPLICIT REAL*8( A-H,O-Z)
C
DIMENSION U ( NX, NY), V ( NX, NY), P ( NX, NY),
1 UN( NNP), VN( NNP), PN( NNP)
C
J = 1
C ---------------------------------
DO 100 IY = 1, NY
DO 100 IX = 1, NX
C
UN( J) = U( IX,IY)
VN( J) = V( IX,IY)
J = J + 1
C
100 CONTINUE
C ---------------------------------
J = 1
DO 200 IY = NY, 1, -1
DO 200 IX = 1, NX
C
PN( J) = P( IX,IY)
J = J + 1
C
200 CONTINUE
C ---------------------------------
C
RETURN
END
C ***********************************************************************
C
SUBROUTINE TTLE
C
WRITE(6,*) '*****************************************************'
WRITE(6,*) '* *'
WRITE(6,*) '* F 2 D - F D M (uvp) *'
WRITE(6,*) '* *'
WRITE(6,*) '* 2D Fluid Flow Analysis *'
WRITE(6,*) '* *'
WRITE(6,*) '* by Fourth-ordered F.D.M. *'
WRITE(6,*) '* *'
WRITE(6,*) '* ( uvp ) *'
WRITE(6,*) '* *'
WRITE(6,*) '* ( V.4.0 ) *'
WRITE(6,*) '* *'
WRITE(6,*) '* Copyright 2011 : Yasuhiro MATSUDA *'
WRITE(6,*) '* *'
WRITE(6,*) '*****************************************************'
C
RETURN
END
C **********************************************************************
C
SUBROUTINE VLSX ( U, V, UA, VA, HX, NX, NY, AX, BX, FU, FV )
C
IMPLICIT REAL*8 ( A-H,O-Z)
C
DIMENSION U ( NX,NY), V ( NX,NY), UA( NX,NY), VA( NX,NY),
1 HX( NX), AX( NX,NY), BX( NX,NY),
2 FU( NX+3), FV( NX+3)
C
COMMON /OUT/ DT, RE, WEI, NPR, LOPE, MSH
COMMON /TES/ A1, A2, A3, A4, A5, A6, A7, A8
C
C --------------------------------------------------------------------
DO 100 J = 2, NY-1
DO 100 I = 2, NX-1
C --------------------------------------------------------------------
C
AX( I,J) = U( I,J)* DT
C ===========================
C
IF ( I.EQ.2.OR.I.EQ.NX-1) GO TO 70
IF ( I.EQ.3) GO TO 70
IF ( I.EQ.4) GO TO 70
GO TO 80
C
C ----- Meshes with Equal length ---
C
** 50 CONTINUE
C
** UA( I,J) = U( I,J)-BX( I,J)/ 2.D0* (U( I+1,J)- U( I-1,J))
** VA( I,J) = V( I,J)-BX( I,J)/ 2.D0* (V( I+1,J)- V( I-1,J))
C
** GO TO 100
C
** 60 CONTINUE
C
** DO 200 M = I, I+1
C
** FU(M) =
** 1 BX( I,J) * ( A1* (U(M-1,J) +U(M,J)) + A5* (U(M-2,J) +U(M+1,J)))
** 2 + BX( I,J)**2* ( A2* (U(M-1,J)- U(M,J)) + A6* (U(M-2,J)- U(M+1,J)))
** 3 + BX( I,J)**3* ( A3* (U(M-1,J) +U(M,J)) + A7* (U(M-2,J) +U(M+1,J)))
** 4 + BX( I,J)**4* ( A4* (U(M-1,J)- U(M,J)) + A8* (U(M-2,J)- U(M+1,J)))
C
** FV(M) =
** 1 BX( I,J) * ( A1* (V(M-1,J) +V(M,J)) + A5* (V(M-2,J) +V(M+1,J)))
** 2 + BX( I,J)**2* ( A2* (V(M-1,J)- V(M,J)) + A6* (V(M-2,J)- V(M+1,J)))
** 3 + BX( I,J)**3* ( A3* (V(M-1,J) +V(M,J)) + A7* (V(M-2,J) +V(M+1,J)))
** 4 + BX( I,J)**4* ( A4* (V(M-1,J)- V(M,J)) + A8* (V(M-2,J)- V(M+1,J)))
C
** 200 CONTINUE
C
** UA( I,J) = U( I,J) + WEI* ( FU( I)-FU( I+1))
** VA( I,J) = V( I,J) + WEI* ( FV( I)-FV( I+1))
** GO TO 100
C
C ----- Meshes of Unequal length ---
C
70 CONTINUE
C
F1 = ( - HX( I)/ HX( I-1)/( HX( I-1) + HX( I)))
F2 = (( HX( I) - HX( I-1))/ HX( I-1)/ HX( I))
F3 = ( HX( I-1)/ HX( I)/( HX( I-1) + HX( I)))
C
UA( I,J) = U( I,J)
1 - AX( I,J)* ( F1* U( I-1,J) + F2* U( I,J) + F3* U( I+1,J))
C
VA( I,J) = V( I,J)
1 - AX( I,J)* ( F1* V( I-1,J) + F2* V( I,J) + F3* V( I+1,J))
C
GO TO 100
C
80 CONTINUE
C
AX1 = HX( I-1)* HX( I)* ( HX( I) + HX( I+1))
1 / HX( I-2)/( HX( I-2) + HX( I-1))
2 /( HX( I-2) + HX( I-1) + HX( I))
3 /( HX( I-2) + HX( I-1) + HX( I) + HX( I+1))
C
AX2 = - 2.D0* ( HX( I-1)* ( 2.D0* HX( I) + HX( I+1))
1 - HX( I)* ( HX( I) + HX( I+1)))
2 / HX( I-2)/( HX( I-2) + HX( I-1))
3 /( HX( I-2) + HX( I-1) + HX( I))
4 /( HX( I-2) + HX( I-1) + HX( I) + HX( I+1))
C
AX3 = 6.D0* HX( I)/ HX( I-2)/( HX( I-2) + HX( I-1))
1 /( HX( I-2) + HX( I-1) + HX( I))
2 /( HX( I-2) + HX( I-1) + HX( I) + HX( I+1))
C
AX4 = 24.D0* ( HX( I-1)- 2.D0* HX( I) - HX( I+1))
1 / HX( I-2)/( HX( I-2) + HX( I-1))
2 /( HX( I-2) + HX( I-1) + HX( I))
3 /( HX( I-2) + HX( I-1) + HX( I) + HX( I+1))
C
AX5 = HX( I-2)* HX( I-1)* ( HX( I-1) + HX( I))
1 / HX( I-3)/( HX( I-3) + HX( I-2))
2 /( HX( I-3) + HX( I-2) + HX( I-1))
3 /( HX( I-3) + HX( I-2) + HX( I-1) + HX( I))
C
AX6 = - 2.D0* ( HX( I-2)* ( 2.D0* HX( I-1) + HX( I))
1 - HX( I-1)* ( HX( I-1) + HX( I)))
2 / HX( I-3)/( HX( I-3) + HX( I-2))
3 /( HX( I-3) + HX( I-2) + HX( I-1))
4 /( HX( I-3) + HX( I-2) + HX( I-1) + HX( I))
C
AX7 = 6.D0* HX( I-1)/ HX( I-3)/( HX( I-3) + HX( I-2))
1 /( HX( I-3) + HX( I-2) + HX( I-1))
2 /( HX( I-3) + HX( I-2) + HX( I-1) + HX( I))
C
AX8 = 24.D0* ( HX( I-2)- 2.D0* HX( I-1) - HX( I))
1 / HX( I-3)/( HX( I-3) + HX( I-2))
2 /( HX( I-3) + HX( I-2) + HX( I-1))
3 /( HX( I-3) + HX( I-2) + HX( I-1) + HX( I))
C
AX9 = HX( I-3)* HX( I-2)* ( HX( I-2) + HX( I-1))
1 / HX( I-4)/( HX( I-4) + HX( I-3))
2 /( HX( I-4) + HX( I-3) + HX( I-2))
3 /( HX( I-4) + HX( I-3) + HX( I-2) + HX( I-1))
C
AX10= - 2.D0* ( HX( I-3)* ( 2.D0* HX( I-2) + HX( I-1))
1 - HX( I-2)* ( HX( I-2) + HX( I-1)))
2 / HX( I-4)/( HX( I-4) + HX( I-3))
3 /( HX( I-4) + HX( I-3) + HX( I-2))
4 /( HX( I-4) + HX( I-3) + HX( I-2) + HX( I-1))
C
AX11 = 6.D0* HX( I-2)/ HX( I-4)/( HX( I-4) + HX( I-3))
1 /( HX( I-4) + HX( I-3) + HX( I-2))
2 /( HX( I-4) + HX( I-3) + HX( I-2) + HX( I-1))
C
AX12 = 24.D0* ( HX( I-3)- 2.D0* HX( I-2) - HX( I-1))
1 / HX( I-4)/( HX( I-4) + HX( I-3))
2 /( HX( I-4) + HX( I-3) + HX( I-2))
3 /( HX( I-4) + HX( I-3) + HX( I-2) + HX( I-1))
C
AX13 = HX( I)* HX( I+1)* ( HX( I+1) + HX( I+2))
1 / HX( I-1)/( HX( I-1) + HX( I))
2 /( HX( I-1) + HX( I) + HX( I+1))
3 /( HX( I-1) + HX( I) + HX( I+1) + HX( I+2))
C
AX14 = - 2.D0* ( HX( I)* ( 2.D0* HX( I+1) + HX( I+2))
1 - HX( I+1)* ( HX( I+1) + HX( I+2)))
2 / HX( I-1)/( HX( I-1) + HX( I))
3 /( HX( I-1) + HX( I) + HX( I+1))
4 /( HX( I-1) + HX( I) + HX( I+1) + HX( I+2))
C
AX15 = 6.D0* HX( I+1)/ HX( I-1)/( HX( I-1) + HX( I))
1 /( HX( I-1) + HX( I) + HX( I+1))
2 /( HX( I-1) + HX( I) + HX( I+1) + HX( I+2))
C
AX16 = 24.D0* ( HX( I)- 2.D0* HX( I+1) - HX( I+2))
1 / HX( I-1)/( HX( I-1) + HX( I))
2 /( HX( I-1) + HX( I) + HX( I+1))
3 /( HX( I-1) + HX( I) + HX( I+1) + HX( I+2))
C
BX1 = ( - AX1
1 * ( HX( I-2) + HX( I-1))* ( HX( I-2) + HX( I-1) + HX( I))
2 * ( HX( I-2) + HX( I-1) + HX( I) + HX( I+1))
3 - HX( I)* ( HX( I) + HX( I+1)))
4 / HX( I-1)/( HX( I-1) + HX( I))
5 /( HX( I-1) + HX( I) + HX( I+1))
C
BX2 = ( - AX2
1 * ( HX( I-2) + HX( I-1))* ( HX( I-2) + HX( I-1) + HX( I))
2 * ( HX( I-2) + HX( I-1) + HX( I) + HX( I+1))
3 + 2.D0* ( 2.D0* HX( I) + HX( I+1)))
4 / HX( I-1)/( HX( I-1) + HX( I))
5 /( HX( I-1) + HX( I) + HX( I+1))
C
BX3 = ( - AX3
1 * ( HX( I-2) + HX( I-1))* ( HX( I-2) + HX( I-1) + HX( I))
2 * ( HX( I-2) + HX( I-1) + HX( I) + HX( I+1))
3 - 6.D0)/ HX( I-1)/( HX( I-1) + HX( I))
4 /( HX( I-1) + HX( I) + HX( I+1))
C
BX4 = ( - AX4
1 * ( HX( I-2) + HX( I-1))* ( HX( I-2) + HX( I-1) + HX( I))
2 * ( HX( I-2) + HX( I-1) + HX( I) + HX( I+1)))
3 / HX( I-1)/( HX( I-1) + HX( I))
4 /( HX( I-1) + HX( I) + HX( I+1))
C
BX5 = ( - AX5
1 * ( HX( I-3) + HX( I-2))* ( HX( I-3) + HX( I-2) + HX( I-1))
2 * ( HX( I-3) + HX( I-2) + HX( I-1) + HX( I))
3 - HX( I-1)* ( HX( I-1) + HX( I)))
4 / HX( I-2)/( HX( I-2) + HX( I-1))
5 /( HX( I-2) + HX( I-1) + HX( I))
C
BX6 = ( - AX6
1 * ( HX( I-3) + HX( I-2))* ( HX( I-3) + HX( I-2) + HX( I-1))
2 * ( HX( I-3) + HX( I-2) + HX( I-1) + HX( I))
3 + 2.D0* ( 2.D0* HX( I-1) + HX( I)))
4 / HX( I-2)/( HX( I-2) + HX( I-1))
5 /( HX( I-2) + HX( I-1) + HX( I))
C
BX7 = ( - AX7
1 * ( HX( I-3) + HX( I-2))* ( HX( I-3) + HX( I-2) + HX( I-1))
2 * ( HX( I-3) + HX( I-2) + HX( I-1) + HX( I)) - 6.D0)
3 / HX( I-2)/( HX( I-2) + HX( I-1))
4 /( HX( I-2) + HX( I-1) + HX( I))
C
BX8 = ( - AX8
1 * ( HX( I-3) + HX( I-2))* ( HX( I-3) + HX( I-2) + HX( I-1))
2 * ( HX( I-3) + HX( I-2) + HX( I-1) + HX( I)))/ HX( I-2)
3 /( HX( I-2) + HX( I-1))/( HX( I-2) + HX( I-1) + HX( I))
C
BX9 = ( - AX9
1 * ( HX( I-4) + HX( I-3))* ( HX( I-4) + HX( I-3) + HX( I-2))
2 * ( HX( I-4) + HX( I-3) + HX( I-2) + HX( I-1))
3 - HX( I-2)* ( HX( I-2) + HX( I-1)))/ HX( I-3)
4 /( HX( I-3) + HX( I-2))/( HX( I-3) + HX( I-2) + HX( I-1))
C
BX10= ( - AX10
1 * ( HX( I-4) + HX( I-3))* ( HX( I-4) + HX( I-3) + HX( I-2))
2 * ( HX( I-4) + HX( I-3) + HX( I-2) + HX( I-1))
3 + 2.D0* ( 2.D0* HX( I-2) + HX( I-1)))/ HX( I-3)
4 /( HX( I-3) + HX( I-2))/( HX( I-3) + HX( I-2) + HX( I-1))
C
BX11 = ( -AX11
1 * ( HX( I-4) + HX( I-3))* ( HX( I-4) + HX( I-3) + HX( I-2))
2 * ( HX( I-4) + HX( I-3) + HX( I-2) + HX( I-1))
3 - 6.D0)/ HX( I-3)/( HX( I-3) + HX( I-2))
4 /( HX( I-3) + HX( I-2) + HX( I-1))
C
BX12 = ( -AX12
1 * ( HX( I-4) + HX( I-3))* ( HX( I-4) + HX( I-3) + HX( I-2))
2 * ( HX( I-4) + HX( I-3) + HX( I-2) + HX( I-1)))
3 / HX( I-3)/( HX( I-3) + HX( I-2))
4 /( HX( I-3) + HX( I-2) + HX( I-1))
C
DX1 = ( AX1* ( HX( I-2) + HX( I-1))
1 * ( HX( I-2) + HX( I-1) + HX( I) + HX( I+1))
2 + BX1* HX( I-1)* ( HX( I-1) + HX( I) + HX( I+1))
3 + HX( I) + HX( I+1))/ HX( I)/ HX( I+1)
C
DX2 = ( AX2* ( HX( I-2) + HX( I-1))
1 * ( HX( I-2) + HX( I-1) + HX( I) + HX( I+1))
2 + BX2* HX( I-1)* ( HX( I-1) + HX( I) + HX( I+1))
3 - 2.D0)/ HX( I)/ HX( I+1)
C
DX3 = ( AX3* ( HX( I-2) + HX( I-1))
1 * ( HX( I-2) + HX( I-1) + HX( I) + HX( I+1))
2 + BX3* HX( I-1)* ( HX( I-1) + HX( I) + HX( I+1)))
3 / HX( I)/ HX( I+1)
C
DX4 = ( AX4* ( HX( I-2) + HX( I-1))
1 * ( HX( I-2) + HX( I-1) + HX( I) + HX( I+1))
2 + BX4* HX( I-1)* ( HX( I-1) + HX( I) + HX( I+1)))
3 / HX( I)/ HX( I+1)
C
DX5 = ( AX5* ( HX( I-3) + HX( I-2))
1 * ( HX( I-3) + HX( I-2) + HX( I-1) + HX( I))
2 + BX5* HX( I-2)* ( HX( I-2) + HX( I-1) + HX( I))
3 + HX( I-1) + HX( I))/ HX( I-1)/ HX( I)
C
DX6 = ( AX6* ( HX( I-3) + HX( I-2))
1 * ( HX( I-3) + HX( I-2) + HX( I-1) + HX( I))
2 + BX6* HX( I-2)* ( HX( I-2) + HX( I-1) + HX( I))
3 - 2.D0)/ HX( I-1)/ HX( I)
C
DX7 = ( AX7* ( HX( I-3) + HX( I-2))
1 * ( HX( I-3) + HX( I-2) + HX( I-1) + HX( I))
2 + BX7* HX( I-2)* ( HX( I-2) + HX( I-1) + HX( I)))
3 / HX( I-1)/ HX( I)
C
DX8 = ( AX8* ( HX( I-3) + HX( I-2))
1 * ( HX( I-3) + HX( I-2) + HX( I-1) + HX( I))
2 + BX8* HX( I-2)* ( HX( I-2) + HX( I-1) + HX( I)))
3 / HX( I-1)/ HX( I)
C
DX9 = ( AX9* ( HX( I-4) + HX( I-3))
1 * ( HX( I-4) + HX( I-3) + HX( I-2) + HX( I-1))
2 + BX9* HX( I-3)* ( HX( I-3) + HX( I-2) + HX( I-1))
3 + HX( I-2) + HX( I-1))/ HX( I-2)/ HX( I-1)
C
DX10= ( AX10* ( HX( I-4) + HX( I-3))
1 * ( HX( I-4) + HX( I-3) + HX( I-2) + HX( I-1))
2 + BX10* HX( I-3)* ( HX( I-3) + HX( I-2) + HX( I-1))
3 - 2.D0)/ HX( I-2)/ HX( I-1)
C
DX11 = ( AX11* ( HX( I-4) + HX( I-3))
1 * ( HX( I-4) + HX( I-3) + HX( I-2) + HX( I-1))
2 + BX11* HX( I-3)* ( HX( I-3) + HX( I-2) + HX( I-1)))
3 / HX( I-2)/ HX( I-1)
C
DX12 = ( AX12* ( HX( I-4) + HX( I-3))
1 * ( HX( I-4) + HX( I-3) + HX( I-2) + HX( I-1))
2 + BX12* HX( I-3)* ( HX( I-3) + HX( I-2) + HX( I-1)))
3 / HX( I-2)/ HX( I-1)
C
EX1 = ( AX1* ( HX( I-2) + HX( I-1)) + BX1* HX( I-1)
1 - DX1* HX( I) +1.D0)/( HX( I) + HX( I+1))
C
EX2 = ( AX2* ( HX( I-2) + HX( I-1)) + BX2* HX( I-1)
1 - DX2* HX( I))/( HX( I) + HX( I+1))
C
EX3 = ( AX3* ( HX( I-2) + HX( I-1)) + BX3* HX( I-1)
1 - DX3* HX( I))/( HX( I) + HX( I+1))
C
EX4 = ( AX4* ( HX( I-2) + HX( I-1)) + BX4* HX( I-1)
1 - DX4* HX( I))/( HX( I) + HX( I+1))
C
EX5 = ( AX5* ( HX( I-3) + HX( I-2)) + BX5* HX( I-2)
1 - DX5* HX( I-1) + 1.D0)/( HX( I-1) + HX( I))
C
EX6 = ( AX6* ( HX( I-3) + HX( I-2)) + BX6* HX( I-2)
1 - DX6* HX( I-1))/( HX( I-1) + HX( I))
C
EX7 = ( AX7* ( HX( I-3) + HX( I-2)) + BX7* HX( I-2)
1 - DX7* HX( I-1))/( HX( I-1) + HX( I))
C
EX8 = ( AX8* ( HX( I-3) + HX( I-2)) + BX8* HX( I-2)
1 - DX8* HX( I-1))/( HX( I-1) + HX( I))
C
EX9 = ( AX9* ( HX( I-4) + HX( I-3)) + BX9* HX( I-3)
1 - DX9* HX( I-2) +1.D0)/( HX( I-2) + HX( I-1))
C
EX10= ( AX10* ( HX( I-4) + HX( I-3)) + BX10* HX( I-3)
1 - DX10* HX( I-2))/( HX( I-2) + HX( I-1))
C
EX11 = ( AX11* ( HX( I-4) + HX( I-3)) + BX11* HX( I-3)
1 - DX11* HX( I-2))/( HX( I-2) + HX( I-1))
C
EX12 = ( AX12* ( HX( I-4) + HX( I-3)) + BX12* HX( I-3)
1 - DX12* HX( I-2))/( HX( I-2) + HX( I-1))
C
AGX1 = - AX1
AGX2 = AX2/ 2.D0
AGX3 = - AX3/ 6.D0
AGX4 = AX4/ 24.D0
C
BGX1 = - ( AX13+ BX1)
BGX2 = ( AX14+ BX2)/ 2.D0
BGX3 = - ( AX15+ BX3)/ 6.D0
BGX4 = ( AX16+ BX4)/ 24.D0
C
CGX1 = DX5 + EX9
CGX2 = - ( DX6 + EX10)/ 2.D0
CGX3 = ( DX7 + EX11)/ 6.D0
CGX4 = - ( DX8 + EX12)/ 24.D0
C
DGX1 = EX5
DGX2 = - EX6/ 2.D0
DGX3 = EX7/ 6.D0
DGX4 = - EX8/ 24.D0
C
C --------------------------------------------------------------------
DO 300 M = I, I + 1
C
C ----------------------------------------------------------
** FU(M) = ! FU( I-1/ 2) & FU( I+1/ 2)
C
FU(M) =
C
1 AX( I,J)* ( AGX1* U(M-2,J) + BGX1* U(M-1,J)
2 + CGX1* U(M,J) + DGX1* U(M+1,J))
3 + AX( I,J)**2 * ( AGX2* U(M-2,J) + BGX2* U(M-1,J)
4 + CGX2* U(M,J) + DGX2* U(M+1,J))
C
5 + AX( I,J)**3 * ( AGX3* U(M-2,J) + BGX3* U(M-1,J)
6 + CGX3* U(M,J) + DGX3* U(M+1,J))
7 + AX( I,J)**4 * ( AGX4* U(M-2,J) + BGX4* U(M-1,J)
8 + CGX4* U(M,J) + DGX4* U(M+1,J))
C
C ----------------------------------------------------------
C
** FV(M) = ! FV( I-1/ 2) & FV( I+1/ 2)
C
FV(M) =
C
1 AX( I,J)* ( AGX1* V(M-2,J) + BGX1* V(M-1,J)
2 + CGX1* V(M,J) + DGX1* V(M+1,J))
3 + AX( I,J)** 2 * ( AGX2* V(M-2,J) + BGX2* V(M-1,J)
4 + CGX2* V(M,J) + DGX2* V(M+1,J))
C
5 + AX( I,J)** 3 * ( AGX3* V(M-2,J) + BGX3* V(M-1,J)
6 + CGX3* V(M,J) + DGX3* V(M+1,J))
7 + AX( I,J)** 4 * ( AGX4* V(M-2,J) + BGX4* V(M-1,J)
8 + CGX4* V(M,J) + DGX4* V(M+1,J))
C
C ----------------------------------------------------------
C
300 CONTINUE
C --------------------------------------------------------------------
C
UA( I,J) = U( I,J) + WEI* ( FU( I) - FU( I+1))
VA( I,J) = V( I,J) + WEI* ( FV( I) - FV( I+1))
C
100 CONTINUE
C --------------------------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE VLSY ( U, V, UA, VA, HY, NX, NY, AY, BY, FU, FV )
C
IMPLICIT REAL*8( A-H,O-Z)
C
DIMENSION U ( NX,NY), V ( NX,NY), UA( NX,NY), VA( NX,NY),
1 HY( NY), FU( NY+3), FV( NY+3), BY( NX,NY),
2 AY( NX,NY)
C
COMMON /OUT/ DT, RE, WEI, NPR, LOPE, MSH
COMMON /TES/ A1, A2, A3, A4, A5, A6, A7, A8
C
C --------------------------------------------------------------------
DO 100 J = 2, NY-1
DO 100 I = 2, NX-1
C --------------------------------------------------------------------
C
AY( I,J) = V( I,J)* DT
C
IF ( J.EQ.2.OR.J.EQ.NY-1) GO TO 70
IF ( J.EQ.3) GO TO 70
IF ( J.EQ.4) GO TO 70
GO TO 80
C
C ----- Meshes with Equal length ---
C
** 50 CONTINUE
C
** UA( I,J) = UA( I,J)-BY( I,J)/ 2.D0* (U( I,J+1)- U( I,J-1))
** VA( I,J) = VA( I,J)-BY( I,J)/ 2.D0* (V( I,J+1)- V( I,J-1))
** GO TO 100
C
** 60 CONTINUE
C
** DO 200 M=J,J+1
C
** FU(M) =
** 1 BY( I,J) * ( A1* (U( I,M-1) +U( I,M)) + A5* (U( I,M-2) +U( I,M+1)))
** 2 + BY( I,J)**2* ( A2* (U( I,M-1)- U( I,M)) + A6* (U( I,M-2)- U( I,M+1)))
** 3 + BY( I,J)**3* ( A3* (U( I,M-1) +U( I,M)) + A7* (U( I,M-2) +U( I,M+1)))
** 4 + BY( I,J)**4* ( A4* (U( I,M-1)- U( I,M)) + A8* (U( I,M-2)- U( I,M+1)))
C
** FV(M) =
** 1 BY( I,J) * ( A1* (V( I,M-1) +V( I,M)) + A5* (V( I,M-2) +V( I,M+1)))
** 2 + BY( I,J)**2* ( A2* (V( I,M-1)- V( I,M)) + A6* (V( I,M-2)- V( I,M+1)))
** 3 + BY( I,J)**3* ( A3* (V( I,M-1) +V( I,M)) + A7* (V( I,M-2) +V( I,M+1)))
** 4 + BY( I,J)**4* ( A4* (V( I,M-1)- V( I,M)) + A8* (V( I,M-2)- V( I,M+1)))
** 200 CONTINUE
C
** UA( I,J) = UA( I,J) + WEI* ( FU( J)-FU( J+1))
** VA( I,J) = VA( I,J) + WEI* ( FV( J)-FV( J+1))
** GO TO 100
C
C
C ----- Meshes of Unequal length ---
C
70 CONTINUE
C
F4 = ( - HY( J)/ HY( J-1)/( HY( J-1) + HY( J)))
F5 = (( HY( J) - HY( J-1))/ HY( J-1)/ HY( J))
F6 = ( HY( J-1)/ HY( J)/( HY( J-1) + HY( J)))
C
UA( I,J) = UA( I,J)
1 - AY( I,J)* ( F4* U( I,J-1) + F5* U( I,J) + F6* U( I,J+1))
C
VA( I,J) = VA( I,J)
1 - AY( I,J)* ( F4* V( I,J-1) + F5* V( I,J) + F6* V( I,J+1))
C
GO TO 100
C
80 AY1 = HY( J-1)* HY( J)* ( HY( J) + HY( J+1))
1 / HY( J-2)/( HY( J-2) + HY( J-1))
2 /( HY( J-2) + HY( J-1) + HY( J))
3 /( HY( J-2) + HY( J-1) + HY( J) + HY( J+1))
C
AY2 = - 2.D0* ( HY( J-1)* ( 2.D0* HY( J) + HY( J+1))
1 - HY( J)* ( HY( J) + HY( J+1)))
2 / HY( J-2)/( HY( J-2) + HY( J-1))
3 /( HY( J-2) + HY( J-1) + HY( J))
4 /( HY( J-2) + HY( J-1) + HY( J) + HY( J+1))
C
AY3 = 6.D0* HY( J)/ HY( J-2)/( HY( J-2) + HY( J-1))
1 /( HY( J-2) + HY( J-1) + HY( J))
2 /( HY( J-2) + HY( J-1) + HY( J) + HY( J+1))
C
AY4 = 24.D0* ( HY( J-1)- 2.D0* HY( J) - HY( J+1))
1 / HY( J-2)/( HY( J-2) + HY( J-1))
2 /( HY( J-2) + HY( J-1) + HY( J))
3 /( HY( J-2) + HY( J-1) + HY( J) + HY( J+1))
C
AY5 = HY( J-2)* HY( J-1)* ( HY( J-1) + HY( J))
1 / HY( J-3)/( HY( J-3) + HY( J-2))
2 /( HY( J-3) + HY( J-2) + HY( J-1))
3 /( HY( J-3) + HY( J-2) + HY( J-1) + HY( J))
C
AY6 = - 2.D0* ( HY( J-2)* ( 2.D0* HY( J-1) + HY( J))
1 - HY( J-1)* ( HY( J-1) + HY( J)))
2 / HY( J-3)/( HY( J-3) + HY( J-2))
3 /( HY( J-3) + HY( J-2) + HY( J-1))
4 /( HY( J-3) + HY( J-2) + HY( J-1) + HY( J))
C
AY7 = 6.D0* HY( J-1)/ HY( J-3)/( HY( J-3) + HY( J-2))
1 /( HY( J-3) + HY( J-2) + HY( J-1))
2 /( HY( J-3) + HY( J-2) + HY( J-1) + HY( J))
C
AY8 = 24.D0* ( HY( J-2)- 2.D0* HY( J-1) - HY( J))
1 / HY( J-3)/( HY( J-3) + HY( J-2))
2 /( HY( J-3) + HY( J-2) + HY( J-1))
3 /( HY( J-3) + HY( J-2) + HY( J-1) + HY( J))
C
AY9 = HY( J-3)* HY( J-2)* ( HY( J-2) + HY( J-1))
1 / HY( J-4)/( HY( J-4) + HY( J-3))
2 /( HY( J-4) + HY( J-3) + HY( J-2))
3 /( HY( J-4) + HY( J-3) + HY( J-2) + HY( J-1))
C
AY10= - 2.D0* ( HY( J-3)* ( 2.D0* HY( J-2) + HY( J-1))
1 - HY( J-2)* ( HY( J-2) + HY( J-1)))
2 / HY( J-4)/( HY( J-4) + HY( J-3))
3 /( HY( J-4) + HY( J-3) + HY( J-2))
4 /( HY( J-4) + HY( J-3) + HY( J-2) + HY( J-1))
C
AY11 = 6.D0* HY( J-2)/ HY( J-4)/( HY( J-4) + HY( J-3))
1 /( HY( J-4) + HY( J-3) + HY( J-2))
2 /( HY( J-4) + HY( J-3) + HY( J-2) + HY( J-1))
C
AY12 = 24.D0* ( HY( J-3)- 2.D0* HY( J-2) - HY( J-1))
1 / HY( J-4)/( HY( J-4) + HY( J-3))
2 /( HY( J-4) + HY( J-3) + HY( J-2))
3 /( HY( J-4) + HY( J-3) + HY( J-2) + HY( J-1))
C
AY13 = HY( J)* HY( J+1)* ( HY( J+1) + HY( J+2))
1 / HY( J-1)/( HY( J-1) + HY( J))
2 /( HY( J-1) + HY( J) + HY( J+1))
3 /( HY( J-1) + HY( J) + HY( J+1) + HY( J+2))
C
AY14 = - 2.D0* ( HY( J)* ( 2.D0* HY( J+1) + HY( J+2))
1 - HY( J+1)* ( HY( J+1) + HY( J+2)))
2 / HY( J-1)/( HY( J-1) + HY( J))
3 /( HY( J-1) + HY( J) + HY( J+1))
4 /( HY( J-1) + HY( J) + HY( J+1) + HY( J+2))
C
AY15 = 6.D0* HY( J+1)/ HY( J-1)/( HY( J-1) + HY( J))
1 /( HY( J-1) + HY( J) + HY( J+1))
2 /( HY( J-1) + HY( J) + HY( J+1) + HY( J+2))
C
AY16 = 24.D0* ( HY( J)- 2.D0* HY( J+1) - HY( J+2))
1 / HY( J-1)/( HY( J-1) + HY( J))
2 /( HY( J-1) + HY( J) + HY( J+1))
3 /( HY( J-1) + HY( J) + HY( J+1) + HY( J+2))
C
BY1 = ( - AY1* ( HY( J-2) + HY( J-1))
1 * ( HY( J-2) + HY( J-1) + HY( J))
2 * ( HY( J-2) + HY( J-1) + HY( J) + HY( J+1))
3 - HY( J)* ( HY( J) + HY( J+1)))
4 / HY( J-1)/( HY( J-1) + HY( J))
5 /( HY( J-1) + HY( J) + HY( J+1))
C
BY2 = ( - AY2* ( HY( J-2) + HY( J-1))
1 * ( HY( J-2) + HY( J-1) + HY( J))
2 * ( HY( J-2) + HY( J-1) + HY( J) + HY( J+1))
3 + 2.D0* ( 2.D0* HY( J) + HY( J+1)))
4 / HY( J-1)/( HY( J-1) + HY( J))
5 / ( HY( J-1) + HY( J) + HY( J+1))
C
BY3 = ( - AY3* ( HY( J-2) + HY( J-1))
1 * ( HY( J-2) + HY( J-1) + HY( J))
2 * ( HY( J-2) + HY( J-1) + HY( J) + HY( J+1))
3 - 6.D0)/ HY( J-1)/( HY( J-1) + HY( J))
4 /( HY( J-1) + HY( J) + HY( J+1))
C
BY4 = ( - AY4* ( HY( J-2) + HY( J-1))
1 * ( HY( J-2) + HY( J-1) + HY( J))
2 * ( HY( J-2) + HY( J-1) + HY( J) + HY( J+1)))
3 / HY( J-1)/( HY( J-1) + HY( J))
4 /( HY( J-1) + HY( J) + HY( J+1))
C
BY5 = ( - AY5* ( HY( J-3) + HY( J-2))
1 * ( HY( J-3) + HY( J-2) + HY( J-1))
2 * ( HY( J-3) + HY( J-2) + HY( J-1) + HY( J))
3 - HY( J-1)* ( HY( J-1) + HY( J)))
4 / HY( J-2)/( HY( J-2) + HY( J-1))
5 /( HY( J-2) + HY( J-1) + HY( J))
C
BY6 = ( - AY6* ( HY( J-3) + HY( J-2))
1 * ( HY( J-3) + HY( J-2) + HY( J-1))
2 * ( HY( J-3) + HY( J-2) + HY( J-1) + HY( J))
3 + 2.D0* ( 2.D0* HY( J-1) + HY( J)))
4 / HY( J-2)/( HY( J-2) + HY( J-1))
5 /( HY( J-2) + HY( J-1) + HY( J))
C
BY7 = ( - AY7* ( HY( J-3) + HY( J-2))
1 * ( HY( J-3) + HY( J-2) + HY( J-1))
2 * ( HY( J-3) + HY( J-2) + HY( J-1) + HY( J))
3 - 6.D0)/ HY( J-2)/( HY( J-2) + HY( J-1))
4 /( HY( J-2) + HY( J-1) + HY( J))
C
BY8 = ( - AY8* ( HY( J-3) + HY( J-2))
1 * ( HY( J-3) + HY( J-2) + HY( J-1))
2 * ( HY( J-3) + HY( J-2) + HY( J-1) + HY( J)))
3 / HY( J-2)/( HY( J-2) + HY( J-1))
4 /( HY( J-2) + HY( J-1) + HY( J))
C
BY9 = ( - AY9* ( HY( J-4) + HY( J-3))
1 * ( HY( J-4) + HY( J-3) + HY( J-2))
2 * ( HY( J-4) + HY( J-3) + HY( J-2) + HY( J-1))
3 - HY( J-2)* ( HY( J-2) + HY( J-1)))
4 / HY( J-3)/( HY( J-3) + HY( J-2))
5 /( HY( J-3) + HY( J-2) + HY( J-1))
C
BY10= ( - AY10* ( HY( J-4) + HY( J-3))
1 * ( HY( J-4) + HY( J-3) + HY( J-2))
2 * ( HY( J-4) + HY( J-3) + HY( J-2) + HY( J-1))
3 + 2.D0* ( 2.D0* HY( J-2) + HY( J-1)))
4 / HY( J-3)/( HY( J-3) + HY( J-2))
5 /( HY( J-3) + HY( J-2) + HY( J-1))
C
BY11 = ( - AY11* ( HY( J-4) + HY( J-3))
1 * ( HY( J-4) + HY( J-3) + HY( J-2))
2 * ( HY( J-4) + HY( J-3) + HY( J-2) + HY( J-1))
3 - 6.D0)/ HY( J-3)/( HY( J-3) + HY( J-2))
4 /( HY( J-3) + HY( J-2) + HY( J-1))
C
BY12 = ( - AY12* ( HY( J-4) + HY( J-3))
1 * ( HY( J-4) + HY( J-3) + HY( J-2))
2 * ( HY( J-4) + HY( J-3) + HY( J-2) + HY( J-1)))
3 / HY( J-3)/( HY( J-3) + HY( J-2))
4 /( HY( J-3) + HY( J-2) + HY( J-1))
C
DY1 = ( AY1* ( HY( J-2) + HY( J-1))
1 * ( HY( J-2) + HY( J-1) + HY( J) + HY( J+1))
2 + BY1* HY( J-1)* ( HY( J-1) + HY( J) + HY( J+1))
3 + HY( J) + HY( J+1))/ HY( J)/ HY( J+1)
C
DY2 = ( AY2* ( HY( J-2) + HY( J-1))
1 * ( HY( J-2) + HY( J-1) + HY( J) + HY( J+1))
2 + BY2* HY( J-1)* ( HY( J-1) + HY( J) + HY( J+1))
3 - 2.D0)/ HY( J)/ HY( J+1)
C
DY3 = ( AY3* ( HY( J-2) + HY( J-1))
1 * ( HY( J-2) + HY( J-1) + HY( J) + HY( J+1))
2 + BY3* HY( J-1)* ( HY( J-1) + HY( J) + HY( J+1)))
3 / HY( J)/ HY( J+1)
C
DY4 = ( AY4* ( HY( J-2) + HY( J-1))
1 * ( HY( J-2) + HY( J-1) + HY( J) + HY( J+1))
2 + BY4* HY( J-1)* ( HY( J-1) + HY( J) + HY( J+1)))
3 / HY( J)/ HY( J+1)
C
DY5 = ( AY5* ( HY( J-3) + HY( J-2))
1 * ( HY( J-3) + HY( J-2) + HY( J-1) + HY( J))
2 + BY5* HY( J-2)* ( HY( J-2) + HY( J-1) + HY( J))
3 + HY( J-1) + HY( J)) / HY( J-1)/ HY( J)
C
DY6 = ( AY6* ( HY( J-3) + HY( J-2))
1 * ( HY( J-3) + HY( J-2) + HY( J-1) + HY( J))
2 + BY6* HY( J-2)* ( HY( J-2) + HY( J-1) + HY( J))
3 - 2.D0)/ HY( J-1)/ HY( J)
C
DY7 = ( AY7* ( HY( J-3) + HY( J-2))
1 * ( HY( J-3) + HY( J-2) + HY( J-1) + HY( J))
2 + BY7* HY( J-2)* ( HY( J-2) + HY( J-1) + HY( J)))
3 / HY( J-1)/ HY( J)
C
DY8 = ( AY8* ( HY( J-3) + HY( J-2))
1 * ( HY( J-3) + HY( J-2) + HY( J-1) + HY( J))
2 + BY8* HY( J-2)* ( HY( J-2) + HY( J-1) + HY( J)))
3 / HY( J-1)/ HY( J)
C
DY9 = ( AY9* ( HY( J-4) + HY( J-3))
1 * ( HY( J-4) + HY( J-3) + HY( J-2) + HY( J-1))
2 + BY9* HY( J-3)* ( HY( J-3) + HY( J-2) + HY( J-1))
3 + HY( J-2) + HY( J-1))/ HY( J-2)/ HY( J-1)
C
DY10= ( AY10* ( HY( J-4) + HY( J-3))
1 * ( HY( J-4) + HY( J-3) + HY( J-2) + HY( J-1))
2 + BY10* HY( J-3)* ( HY( J-3) + HY( J-2) + HY( J-1))
3 - 2.D0)/ HY( J-2)/ HY( J-1)
C
DY11 = ( AY11* ( HY( J-4) + HY( J-3))
1 * ( HY( J-4) + HY( J-3) + HY( J-2) + HY( J-1))
2 + BY11* HY( J-3)* ( HY( J-3) + HY( J-2) + HY( J-1)))
3 / HY( J-2)/ HY( J-1)
C
DY12 = ( AY12* ( HY( J-4) + HY( J-3))
1 * ( HY( J-4) + HY( J-3) + HY( J-2) + HY( J-1))
2 + BY12* HY( J-3)* ( HY( J-3) + HY( J-2) + HY( J-1)))
3 / HY( J-2)/ HY( J-1)
C
EY1 = ( AY1* ( HY( J-2) + HY( J-1)) + BY1* HY( J-1)
1 - DY1* HY( J) +1.D0)/( HY( J) + HY( J+1))
C
EY2 = ( AY2* ( HY( J-2) + HY( J-1)) + BY2* HY( J-1)
1 - DY2* HY( J))/( HY( J) + HY( J+1))
C
EY3 = ( AY3* ( HY( J-2) + HY( J-1)) + BY3* HY( J-1)
1 - DY3* HY( J))/( HY( J) + HY( J+1))
C
EY4 = ( AY4* ( HY( J-2) + HY( J-1)) + BY4* HY( J-1)
1 - DY4* HY( J))/( HY( J) + HY( J+1))
C
EY5 = ( AY5* ( HY( J-3) + HY( J-2)) + BY5* HY( J-2)
1 - DY5* HY( J-1) +1.D0)/( HY( J-1) + HY( J))
C
EY6 = ( AY6* ( HY( J-3) + HY( J-2)) + BY6* HY( J-2)
1 - DY6* HY( J-1))/( HY( J-1) + HY( J))
C
EY7 = ( AY7* ( HY( J-3) + HY( J-2)) + BY7* HY( J-2)
1 - DY7* HY( J-1))/( HY( J-1) + HY( J))
C
EY8 = ( AY8* ( HY( J-3) + HY( J-2)) + BY8* HY( J-2)
1 - DY8* HY( J-1))/( HY( J-1) + HY( J))
C
EY9 = ( AY9* ( HY( J-4) + HY( J-3)) + BY9* HY( J-3)
1 - DY9* HY( J-2) +1.D0)/( HY( J-2) + HY( J-1))
C
EY10= ( AY10* ( HY(J-4) + HY(J-3)) + BY10* HY(J-3)
1 - DY10* HY( J-2))/( HY( J-2) + HY( J-1))
C
EY11 = ( AY11* ( HY( J-4) + HY( J-3))
1 + BY11* HY( J-3) - DY11* HY( J-2))
2 /( HY( J-2) + HY( J-1))
C
EY12 = ( AY12* ( HY( J-4) + HY( J-3)) + BY12* HY( J-3)
1 - DY12* HY( J-2))/( HY( J-2) + HY( J-1))
C
AGY1 = -AY1
AGY2 = AY2/ 2.D0
AGY3 = -AY3/ 6.D0
AGY4 = AY4/ 24.D0
C
BGY1 = -( AY13 + BY1)
BGY2 = ( AY14 + BY2)/ 2.D0
BGY3 = -( AY15 + BY3)/ 6.D0
BGY4 = ( AY16 + BY4)/ 24.D0
C
CGY1 = DY5 + EY9
CGY2 = -( DY6 + EY10)/ 2.D0
CGY3 = ( DY7 + EY11)/ 6.D0
CGY4 = -( DY8 + EY12)/ 24.D0
C
DGY1 = EY5
DGY2 = - EY6/ 2.D0
DGY3 = EY7/ 6.D0
DGY4 = - EY8/ 24.D0
C
C ----------------------------------------------------------
DO 300 M = J, J + 1
C
FU(M) =
C
1 AY( I,J) * ( AGY1* U( I,M-2) + BGY1* U( I,M-1)
2 + CGY1* U( I,M) + DGY1* U( I,M+1))
3 + AY( I,J)** 2* ( AGY2* U( I,M-2) + BGY2* U( I,M-1)
4 + CGY2* U( I,M) + DGY2* U( I,M+1))
C
5 + AY( I,J)** 3* ( AGY3* U( I,M-2) + BGY3* U( I,M-1)
6 + CGY3* U( I,M) + DGY3* U( I,M+1))
7 + AY( I,J)** 4* ( AGY4* U( I,M-2) + BGY4* U( I,M-1)
8 + CGY4* U( I,M) + DGY4* U( I,M+1))
C
FV(M) =
C
1 AY( I,J) * ( AGY1* V( I,M-2) + BGY1* V( I,M-1)
2 + CGY1* V( I,M) + DGY1* V( I,M+1))
3 + AY( I,J)** 2* ( AGY2* V( I,M-2) + BGY2* V( I,M-1)
4 + CGY2* V( I,M) + DGY2* V( I,M+1))
C
5 + AY( I,J)** 3* ( AGY3* V( I,M-2) + BGY3* V( I,M-1)
6 + CGY3* V( I,M) + DGY3* V( I,M+1))
7 + AY( I,J)** 4* ( AGY4* V( I,M-2) + BGY4* V( I,M-1)
8 + CGY4* V( I,M) + DGY4* V( I,M+1))
C
300 CONTINUE
C ----------------------------------------------------------
C
UA( I,J) = UA( I,J) + WEI* ( FU( J) - FU( J+1))
VA( I,J) = VA( I,J) + WEI* ( FV( J) - FV( J+1))
C
100 CONTINUE
C --------------------------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE VLSF ( U, V, UA, VA, PA, HX, HY, NX, NY )
C
IMPLICIT REAL*8( A-H,O-Z)
C
DIMENSION U ( NX,NY), V ( NX,NY), UA( NX,NY), VA( NX,NY),
1 PA( NX,NY), HX( NX), HY( NY)
C
COMMON /OUT/ DT, RE, WEI, NPR, LOPE, MSH
C
C --------------------------------------------------------------------
DO 100 J = 2,NY-1
DO 100 I = 2,NX-1
C --------------------------------------------------------------------
C
F = 1.D0 + HX( I)* HY( J)
1 /( HX( I)** 2 + HY( J)** 2 )* U( I,J)* V( I,J)* RE* DT
C
C --------------------------------------------------------------------
C
E1 = ( 2.D0/ HX( I-1)/( HX( I-1) + HX( I)))
E2 = ( - 2.D0/ HX( I-1)/ HX( I))
E3 = ( 2.D0/ HX( I)/( HX( I-1) + HX( I)))
E4 = ( 2.D0/ HY( J-1)/( HY( J-1) + HY( J)))
E5 = ( - 2.D0/ HY( J-1)/ HY( J))
E6 = ( 2.D0/ HY( J)/( HY( J-1) + HY( J)))
C
F1 = ( - HX( I)/ HX( I-1)/( HX( I-1) + HX( I)))
F2 = (( HX( I) - HX( I-1))/ HX( I-1)/ HX( I))
F3 = ( HX( I-1)/ HX( I)/( HX( I-1) + HX( I)))
F4 = ( - HY( J)/ HY( J-1)/( HY( J-1) + HY( J)))
F5 = (( HY( J) - HY( J-1))/ HY( J-1)/ HY( J))
F6 = ( HY( J-1)/ HY( J)/( HY( J-1) + HY( J)))
C
C --------------------------------------------------------------------
UA( I,J) = UA( I,J)
C
1 + F* DT/ RE* ( E1* U( I-1,J) + E2* U ( I,J) + E3* U ( I+1,J)
2 + E4* U( I,J-1) + E5* U ( I,J) + E6* U ( I,J+1))
3 - DT* ( F1* PA( I-1,J) + F2* PA( I,J) + F3* PA( I+1,J))
C
C --------------------------------------------------------------------
VA( I,J) = VA( I,J)
C
1 + F* DT/ RE* ( E1* V ( I-1,J) + E2* V ( I,J) + E3* V ( I+1,J)
2 + E4* V ( I,J-1) + E5* V ( I,J) + E6* V ( I,J+1))
3 - DT* ( F4* PA( I,J-1) + F5* PA( I,J) + F6* PA( I,J+1))
C
C --------------------------------------------------------------------
C
100 CONTINUE
C --------------------------------------------------------------------
C
RETURN
END
C
C **********************************************************************
C
SUBROUTINE PRVL ( U, V, UA, VA, P, PA, DP, OU, OV,
1 HX, HY, NX, NY )
C
IMPLICIT REAL*8( A-H,O-Z)
C
DIMENSION U( NX,NY), V ( NX,NY), UA( NX,NY), VA( NX,NY),
1 P( NX,NY), PA( NX,NY), DP( NX,NY),
2 OU( NX,NY), OV( NX,NY), HX( NX), HY( NY)
C
COMMON /OUT/ DT, RE, WEI, NPR, LOPE, MSH
C
C --------------------------------------
DO 100 J = 1, NY
DO 100 I = 1, NX
C
P( I,J) = PA( I,J) + DP( I,J)
C
100 CONTINUE
C ---------------------------------------------------------------
DO 200 J = 2,NY-1
C
DO 200 I = 2,NX-1
C
OU( I,J) = U( I,J)
OV( I,J) = V( I,J)
C
F1 = ( - HX( I)/ HX( I-1)/( HX( I-1) + HX( I)))
F2 = (( HX( I) - HX( I-1))/ HX( I-1)/ HX( I))
F3 = ( HX( I-1)/ HX( I)/( HX( I-1) + HX( I)))
F4 = ( - HY( J)/ HY( J-1)/( HY( J-1) + HY( J)))
F5 = (( HY( J) - HY( J-1))/ HY( J-1)/ HY( J))
F6 = ( HY( J-1)/ HY( J)/( HY( J-1) + HY( J)))
C
U( I,J) = UA( I,J)
1 - DT* ( F1* DP( I-1,J) + F2* DP( I,J) + F3* DP( I+1,J))
C
V( I,J) = VA( I,J)
1 - DT* ( F4* DP( I,J-1) + F5* DP( I,J) + F6* DP( I,J+1))
C
200 CONTINUE
C ---------------------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE FVLB ( U, V, UB, VB, HX, HY, NX, NY)
C
IMPLICIT REAL*8 ( A-H,O-Z)
C
DIMENSION U ( NX,NY), V ( NX,NY), UB( NX,NY), VB( NX,NY),
1 HX( NX), HY( NY)
C
COMMON /OUT/ DT, RE, WEI, NPR, LOPE, MSH
C
C --------------------------------------------------------------------
DO 100 J = 2, NY-1
DO 100 I = 2, NX-1
C --------------------------------------------------------------------
C
C UB( I,J) = U( I,J)
C 1 + DT* (1.D0/RE/ HX( I)**2* (U( I+1,J)- 2.D0* U( I,J) +U( I-1,J))
C
C 2 +1.D0/RE/ HY( J)**2* (U( I,J+1)- 2.D0* U( I,J) +U( I,J-1))
C 3 -0.5D0/ HX( I)* U( I,J)* (U( I+1,J)- U( I-1,J))
C 4 -0.5D0/ HY( J)* V( I,J)* (U( I,J+1)- U( I,J-1)))
C
C VB( I,J) = V( I,J)
C 1 + DT* (1.D0/RE/ HX( I)**2* (V( I+1,J)- 2.D0* V( I,J) +V( I-1,J))
C
C 2 +1.D0/RE/ HY( J)**2* (V( I,J+1)- 2.D0* V( I,J) +V( I,J-1))
C 3 -0.5D0/ HX( I)* U( I,J)* (V( I+1,J)- V( I-1,J))
C 4 -0.5D0/ HY( J)* V( I,J)* (V( I,J+1)- V( I,J-1)))
C
C
E1 = ( 2.D0/ HX( I-1)/( HX( I-1) + HX( I)))
E2 = ( - 2.D0/ HX( I-1)/ HX( I))
E3 = ( 2.D0/ HX( I)/( HX( I-1) + HX( I)))
E4 = ( 2.D0/ HY( J-1)/( HY( J-1) + HY( J)))
E5 = ( - 2.D0/ HY( J-1)/ HY( J))
E6 = ( 2.D0/ HY( J)/( HY( J-1) + HY( J)))
C
F1 = ( -HX( I)/ HX( I-1)/( HX( I-1) + HX( I)))
F2 = (( HX( I) - HX( I-1))/ HX( I-1)/ HX( I))
F3 = ( HX( I-1)/ HX( I)/( HX( I-1) + HX( I)))
F4 = ( -HY( J)/ HY( J-1)/( HY( J-1) + HY( J)))
F5 = (( HY( J) - HY( J-1))/ HY( J-1)/ HY( J))
F6 = ( HY( J-1)/ HY( J)/( HY( J-1) + HY( J)))
C
C --------------------------------------------------------------------
UB( I,J) = U( I,J)
C
1 + DT* (( E1* U( I-1,J) + E2* U( I,J) + E3* U( I+1,J)
2 + E4* U( I,J-1) + E5* U( I,J) + E6* U( I,J+1))/ RE
C
3 - U( I,J)* ( F1* U( I-1,J) + F2* U( I,J) + F3* U( I+1,J))
4 - V( I,J)* ( F4* U( I,J-1) + F5* U( I,J) + F6* U( I,J+1)))
C
C --------------------------------------------------------------------
VB( I,J) = V( I,J)
C
1 + DT* (( E1* V( I-1,J) + E2* V( I,J) + E3* V( I+1,J)
2 + E4* V( I,J-1) + E5* V( I,J) + E6* V( I,J+1))/ RE
C
3 - U( I,J)* ( F1* V( I-1,J) + F2* V( I,J) + F3* V( I+1,J))
4 - V( I,J)* ( F4* V( I,J-1) + F5* V( I,J) + F6* V( I,J+1)))
C
C ---------------------------------------------------------------
C
100 CONTINUE
C --------------------------------------------------------------------
C
RETURN
END
C*********************************************************************
C
SUBROUTINE SINP
C
IMPLICIT REAL*8( A-H,O-Z)
C
COMMON /VRD/ AP
COMMON /OUT/ DT, RE, WEI, NPR, LOPE, MSH
C
C ------------------------------------------------
WRITE (6,2000)
2000 FORMAT(/' ',' DT = ?')
READ(5,*) DT
C
* DT = 0.014D0 ! 20x20 Re = 100 ( Diverg. for DT = 0.015D0
C
* DT = 0.007D0 ! 40x40 Re = 400
* DT = 0.009D0 ! 40x40 Re = 400
C
* DT = 0.02D0 ! UE10 Re=100 XXX because of 4th-ordered FDM
C =============================================
WRITE (6,2100)
2100 FORMAT(' ',' Reynolds Number = ? ')
READ (5,*) RE
C
* RE = 100.D0
C
* RE = 400.D0
C
C ------------------------------------------------
WRITE (6,2200)
2200 FORMAT(' ',' Timing of Print out = ? ')
*** READ (5,*) NPR
C
NPR = 100
C
C ------------------------------------------------
WRITE (6,2300)
2300 FORMAT(' ',' Loop = ? ')
*** READ (5,*) LOPE
C
LOPE = 50
LOPE = 1000
LOPE = 4000
* LOPE = 1000
C
C ------------------------------------------------
WRITE(6,*) ' AP = ?'
** READ(5,*) AP
AP = 1.D0
AP = 1.8D0
AP = 1.4D0 ! Optimised Value
C ==================================
C
C ----- Weighting Parameter ----------------------
C
WEI = 1.D0
C ------------------------------------------------
C
WRITE(6,2400) RE, DT, LOPE, WEI, AP, NPR
2400 FORMAT(/,' *** Re=',F6.1,' DT=',F6.4,' LOPE=',I5,
1 ' WEI=',F5.2,' Accel.Par.=',F7.3,' NPR=',I3,/)
C
RETURN
END
C **********************************************************************
C
SUBROUTINE SORP ( SU, SV, SP, HX, HY, NX, NY, ITT, B,
1 DX, DY, SSP, LOP )
C
IMPLICIT REAL*8( A-H,O-Z)
C
DIMENSION SU( NX,NY), SV( NX,NY), SP( NX,NY), SSP( NX,NY),
1 EU( 4), EP( 4), HX( NX), HY ( NY),
2 DX( NX), DY( NY), B ( NX,NY)
C
COMMON /VRD/ AP
COMMON /OUT/ DT, RE, WEI, NPR, LOPE, MSH
C
ITT = 0
ORM = 0.D0
C
C ----- Boundary condition ---
C
C --------------------------------------
DO 100 I = 2, NX-1
C
SP( I, 1) = SP( I,2)
SP( I,NY) = SP( I,NY-1)
C
100 CONTINUE
C --------------------------------------
DO 200 J = 2, NY-1
C
SP( 1, J) = SP( 2,J)
SP( NX,J) = SP( NX-1,J)
C
200 CONTINUE
C --------------------------------------
SP( 1, 1) = SP( 2,2)
SP( NX, 1) = SP( NX-1,2)
SP( 1, NY) = SP( 2,NY-1)
SP( NX,NY) = SP( NX-1, NY-1)
C
C ----------------------------------------------------------
999 ITT = ITT + 1
ITN = 0
C
C ------------------------------------------------
DO 300 J = 1, NY
DO 300 I = 1, NX
C ------------------------------------------------
C
IF ( I.EQ.1) THEN
EU(1) = SU( I,J)
EP(1) = SP( I,J)
ELSE
EU(1) = SU( I-1,J)
EP(1) = SP( I-1,J)
C
END IF
C --------------------------------------
IF ( I.EQ.NX) THEN
EU( 2) = SU( I,J)
EP( 2) = SP( I,J)
ELSE
EU( 2) = SU( I+1,J)
EP( 2) = SP( I+1,J)
C
END IF
C --------------------------------------
C
C ----- Down ---
C
IF ( J.EQ.1) THEN
EU(3) = SV( I,J)
EP(3) = SP( I,J)
ELSE
EU(3) = SV( I,J-1)
EP(3) = SP( I,J-1)
END IF
C ----- Upper ---
C
IF ( J.EQ.NY) THEN
EU(4) = SV( I,J)
EP(4) = SP( I,J)
ELSE
EU(4) = SV( I,J+1)
EP(4) = SP( I,J+1)
END IF
C --------------------------------------
C
DX( I) =1.D0/ HX( I)** 2
DY( J) =1.D0/ HY( J)** 2
C
C --------------------------------------
IF ( I.EQ.1) THEN
II = I + 1
ELSE
II = I
END IF
C --------------------------------------
IF ( J.EQ.1) THEN
JJ = J + 1
ELSE
JJ = J
END IF
C --------------------------------------
C
E1 = ( 2.D0/ HX( II-1)/( HX( II-1) + HX( I)))
E2 = ( - 2.D0/ HX( II-1)/ HX( I))
E3 = ( 2.D0/ HX( I)/( HX( II-1) + HX( I)))
E4 = ( 2.D0/ HY( JJ-1)/( HY( JJ-1) + HY( J)))
E5 = ( - 2.D0/ HY( JJ-1)/ HY( J))
E6 = ( 2.D0/ HY( J)/( HY( JJ-1) + HY( J)))
C
F1 = ( -HX( I)/ HX( II-1)/( HX( II-1) + HX( I)))
F2 = (( HX( I) - HX( II-1))/ HX( II-1)/ HX( I))
F3 = ( HX( II-1)/ HX( I)/( HX( II-1) + HX( I)))
F4 = ( -HY( J)/ HY( JJ-1)/( HY( JJ-1) + HY( J)))
F5 = (( HY( J) - HY( JJ-1))/ HY( JJ-1)/ HY( J))
F6 = ( HY( JJ-1)/ HY( J)/( HY( JJ-1) + HY( J)))
C
B( I,J) = ( F1* EU( 1) + F2* SU( I,J) + F3* EU( 2)
1 + F4* EU( 3) + F5* SV( I,J) + F6* EU( 4))/ DT
C
SSP( I,J) = ( - E1* EP( 1) - E3* EP( 2) - E4* EP( 3)
1 - E6* EP( 4) + B( I,J))/( E2 + E5)
C
C--------------------------------------------------------
C
ORM = DABS( SSP( I,J) - SP( I,J))
C
SP( I,J) = SP( I,J) + AP* (SSP( I,J) - SP( I,J))
C =====================================================
C
IF ( ORM.GT.1.D-4) ITN =1
B( I,J) = B( I,J)/ DX( I)
C
300 CONTINUE
C ------------------------------------------------
C
IF ( ITN.EQ.0) RETURN
IF ( ITT.LT.( NX-1)*50) GO TO 999
C
C ----------------------------------------------------------
WRITE(6,2000) ITT, ORM, LOP
2000 FORMAT(' * Divergence !',' ITT at SORP =',I5,', ORM=',F10.6,
1 ' Loop=',I5,' * STOP *')
STOP
C
END
C **********************************************************************