戻る

C
PROGRAM MAIN
C **************************************************************** C
C C
C F 2 D - D 4 2 C
C C
C Cavity Fluid Flow Analysis C
C C
C by 4th order & 2nd order Finite Diference Method C
C C
C ( Equal length meshes : DX = DY = DS ) C
C =============== C
C C
C ( Version 4.0 ) C
C C
C Copyright : Yasuhiro MATSUDA C
C C
C **************************************************************** C
C
IMPLICIT REAL*8 (A-H,O-Z)
C
CHARACTER*8 FSM
C
C ----- ( 1) ---
C
PARAMETER ( IX= 11, IY= 11, IX1=IX-1, IY1=IY-1 ) ! 10 X 10
*** PARAMETER ( IX= 21, IY= 21, IX1=IX-1, IY1=IY-1 ) ! 20 X 20
*** PARAMETER ( IX= 41, IY= 41, IX1=IX-1, IY1=IY-1 ) ! 40 X 40
C --------------
C
LOGICAL*1 WPT( IX,IY,4), PPT( IX,IY,4), EPT( IX,IY,4)
C
DIMENSION PSF( IX, IY), WP1( IX, IY), WB ( IX, IY),
1 PS ( IX, IY), S ( IX, IY), U ( IX, IY),
2 V ( IX, IY), UU ( IX*IY), VV ( IX*IY),
3 UB ( IX, IY), VB ( IX, IY), PE ( IX*IY),
C
4 XX ( IX*IY), YY ( IX*IY), XY ( 2,IX*IY),
5 RNU( IX, IY), UN ( IX1*IY1), VN ( IX1*IY1),
6 XN ( IX1*IY1), YN ( IX1*IY1), NDE( IX1*IY1,4)
C
COMMON /HIST/ WMX( 5000), DMX( 5000)
COMMON /DAT1/ LOP, RE, NLM, IPR, NNP, NEL
COMMON /DAT2/ DT, DS, DTS, COF, UL, UV, VIS, WEI
C
CALL TTLE
C
C ----- ( 2) ---
C
OPEN ( 7, FILE='F2DD42_10_FIG.DAT' )
*** OPEN ( 7, FILE='F2DD42_20_FIG.DAT' )
*** OPEN ( 7, FILE='F2DD42_40_FIG.DAT' )
C --------------
C
9999 CONTINUE
C
C ----- ( 3) ---
C
OPEN ( 8, FILE='F2DD42_10.DAT' )
*** OPEN ( 8, FILE='F2DD42_40.DAT' )
*** OPEN ( 8, FILE='F2DD42_20.DAT' )
C
LOP = 0
C
CALL BEGN ( WPT, PPT, EPT, WP1, PSF, S, PS, U, V, UB,
1 VB, WB, FSM, ISC, ICR, IX, IY, NBC )
C
CALL GEOI ( WP1, PSF, WPT, PPT, EPT, IX, IY )
C
CALL STND ( XX, YY, XN, YN, NDE, IX, IY, IX1, IY1 )
C
C --------------------------------------------------------------------
900 LOP = LOP + 1
C
CALL SFCL ( PSF, PS, S, WP1, PPT, ADF, IX, IY )
C
IF ( NBC.EQ.0) CALL VORB0 ( WP1, PSF, WPT, IX, IY )
C
IF ( NBC.EQ.1) CALL VORB1 ( WP1, PSF, WPT, IX, IY )
C
CALL CMSB ( S, WP1, PSF, U, V, RNU, PS, WPT, EPT,
1 IX, IY, ISC, ICR )
C
CALL VRVL ( U, V, UB, VB, UMX, VMX, IX, IY )
C
CALL VRVR ( WP1, WB, PSF, PL, WL1, WL2, VRV, FVV, IX, IY )
C
CALL CHCK ( U, V, UU, VV, PSF, PE, XX, YY, UN, VN,
1 XN, YN, XY, NDE, FSM, ADF, IX, IY, IX1, IY1,
2 *333 )
C
IF ( LOP.LT.NLM) GO TO 900
C
333 CONTINUE
C
CALL RPLT ( U, V, UU, VV, UN, VN, PSF, PE, XX, YY,
1 NDE, XY, XN, YN, FSM, IX, IY, IX1, IY1,
2 NBC )
C
CLOSE ( 8)
C
WRITE(6, 2000)
2000 FORMAT(/,' *** Simulation to be continued ? ( 1: Yes, 0: No)')
** READ(5,*) ICO
ICO = 0
IF ( ICO.EQ.1) GO TO 9999
C
CLOSE ( 7)
C
STOP
END
C **********************************************************************
C
SUBROUTINE DINP ( FSM, IX, IY, ISC, ICR, NBC )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
CHARACTER*8 FSM
C
COMMON /DAT1/ LOP, RE, NLM, IPR, NNP, NEL
COMMON /DAT2/ DT, DS, DTS, COF, UL, UV, VIS, WEI
C
WRITE(6,2000)
2000 FORMAT(/,' * Please input data ')
WRITE(6,2100)
2100 FORMAT('1: FWA(C) 2: SWA(C)'/,' * Please select the scheme'/)
C
** READ(5,*) IFSC
IFSC = 1
C
C ---------------------------------
C
IF ( IFSC.EQ.1) GO TO 22
IF ( IFSC.EQ.2) GO TO 44
C
C ---------------------------------
22 ISC = 0
ICR = 1 ! ICR : Key for Correction Method
FSM = '* FWA(C)'
GO TO 20
C
44 ISC = 1
ICR = 1
FSM = '* SWA(C)'
GO TO 20
C
20 CONTINUE
WRITE(6,2200)
2200 FORMAT(' *** Please input a weighting parameter (0〜1)')
*** READ(5,*) WEI
C
C ----- Attention ---
C
WEI = 0.D0
C ================
WRITE(6,2300)
2300 format(' *** Please input time step')
C
** READ(5,*) DT
DT = 0.1D0
C
WRITE(6,2400)
2400 FORMAT(' *** How many Loops do you want ?')
** READ(5,*) NLM
NLM = 150
C
WRITE(6,2500)
2500 format(/' *** Please input Reynolds number')
** READ(5,*) RE
RE = 100.D0
C
VIS = 1.D0/ RE
C ===================
C
UL = 0.D0
C
UV = 1.D0
C
C ***** Cavity Fluid Flow Analysis **********
C
DS = 1.D0/ DFLOAT( IX-1) ! DX = DY = DS
C ===== Attention ===============================
C
DTS = DT/ DS
COF = VIS* DT/ ( DS* DS )
C
WRITE(6,2600) DS, DTS, COF
2600 FORMAT(/,' * DS =', F7.3,' DTS =', F7.3,' COF =', F8.4)
C
WRITE(6,2700)
2700 FORMAT(/,' *** Approx.of Boundary Condition : = 1 ( By Ghia ), = 0
1 ( By Wood )')
READ(5,*) NBC
WRITE(6,2800) NBC
2800 FORMAT(' ** NBC =',I3)
C
RETURN
END
C **********************************************************************
C
SUBROUTINE TTLE
C
WRITE(6,20)
WRITE(6,21)'*****************************************************'
WRITE(6,21)' '
WRITE(6,21)' ##### ### #### ### # ### '
WRITE(6,21)' # # # # # # # ## # # '
WRITE(6,21)' #### # # # #### # # # # # '
WRITE(6,21)' # ## # # # # ###### ## '
WRITE(6,21)' ## ##### #### ### # ##### '
WRITE(6,21)' '
WRITE(6,21)' ( Version 4.0 ) '
WRITE(6,21)' '
WRITE(6,21)' for Two Dimensional Viscous Fluid Flow '
WRITE(6,21)' '
WRITE(6,21)' by the 4th-ordered and the 2nd-ordered FDM '
WRITE(6,21)' '
WRITE(6,21)' ( ** Equal mesh size : DX = DY = DS ** ) '
WRITE(6,21)' '
WRITE(6,21)' Copyright 2011 : Yasuhiro MATSUDA '
WRITE(6,21)' '
WRITE(6,21)'*****************************************************'
C
20 FORMAT(/)
21 FORMAT(A70)
C
RETURN
END
C ***********************************************************************
C
SUBROUTINE RPLT ( U, V, UU, VV, UN, VN, PSF, PE, XX,
1 YY, NDE, XY, XN, YN, FSM, IX, IY, IX1,
2 IY1, NBC )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
CHARACTER*8 FSM
C
COMMON /HIST/ WMX( 5000), DMX( 5000)
COMMON /DAT1/ LOP, RE, NLM, IPR, NNP, NEL
COMMON /DAT2/ DT, DS, DTS, COF, UL, UV, VIS, WEI
C
DIMENSION U ( IX, IY), V ( IX, IY), UU( IX*IY),
1 VV( IX*IY), PSF( IX, IY), PE( IX*IY),
2 XX( IX*IY), YY ( IX*IY), XY( 2,IX*IY),
3 UN( IX1*IY1), VN ( IX1*IY1), XN( IX1*IY1),
4 YN( IX1*IY1), NDE( IX1*IY1,4)
C
C ----------------------------------------------------------
C
CALL STMP ( U, UU, IX, IY )
C
CALL STMP ( V, VV, IX, IY )
C
CALL STMP ( PSF, PE, IX, IY )
C
CALL STND ( XX, YY, XN, YN, NDE, IX, IY, IX1, IY1 )
C
C ----------------------------------------------------------
DO 100 I = 1, NEL
C
UN( I) = ( UU( NDE( I,1)) + UU( NDE( I,2)))/ 2.D0
VN( I) = ( VV( NDE( I,1)) + VV( NDE( I,2)))/ 2.D0
C
100 CONTINUE
C ----------------------------------------------------------
C
K = 0
C ----------------------------------
DO 200 J = 1, IY
DO 200 I = 1, IX
C
K = K + 1
II = IX* (IY-J) + I
XY( 1,K) = XX( II)
XY( 2,K) = YY( II)
C
200 CONTINUE
C ----------------------------------
C
WRITE(7,2000) FSM, WEI, RE, DT, LOP, NEL, NNP, NBC
2000 FORMAT(A8,3F9.3,3I4,I2)
WRITE(7,2100) ( XN( I), YN( I), I = 1, NEL )
WRITE(7,2100) ( UN( I), VN( I), I = 1, NEL )
2100 FORMAT(10F8.4)
C
VLMX = 0.D0
DO 300 I = 1, NEL
C
IF ( DABS( UN( I)).GE.VLMX ) VLMX = DABS( UN( I))
300 CONTINUE
C
WRITE(6,2200) VLMX
2200 FORMAT(/,' *** Umax.=',F7.3)
WRITE(7,2300) ( WMX( I)*100.D0, I = 2, LOP )
WRITE(7,2300) ( DMX( I)*100.D0, I = 2, LOP )
2300 FORMAT(10F8.2)
C
WRITE(7,2300) ( ( XY( I,J), I = 1, 2 ), J = 1, NNP )
WRITE(7,2300) ( PE( I), I = 1, NNP )
WRITE(7,2300) ( ( U( I,J), I = 1, IX ), J = 1, IY )
WRITE(7,2300) ( ( V( I,J), I = 1, IX ), J = 1, IY )
C
C ---------------------------------------
C
CALL STMP ( PSF, PE, IX, IY )
C
CALL STMP ( U, UU, IX, IY )
C
CALL STMP ( V, VV, IX, IY )
C
C ---------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE CHCK ( U, V, UU, VV, PSF, PE, XX, YY, UN,
1 VN, XN, YN, XY, NDE, FSM, ADF, IX, IY,
2 IX1, IY1, * )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
CHARACTER*8 FSM
C
COMMON /HIST/ WMX( 5000), DMX( 5000)
COMMON /DAT1/ LOP, RE, NLM, IPR, NNP, NEL
COMMON /DAT2/ DT, DS, DTS, COF, UL, UV, VIS, WEI
C
DIMENSION U ( IX,IY), V ( IX, IY), UU ( IX*IY),
1 VV( IX*IY), PSF( IX, IY), PE ( IX*IY),
2 XX( IX*IY), YY ( IX*IY),
3 UN( IX1*IY1), VN ( IX1*IY1), XN ( IX1*IY1),
4 YN( IX1*IY1), XY ( 2,IX*IY), NDE( IX1*IY1,4)
C
RFR = WMX( LOP)* 100.D0
RFV = DMX( LOP)* 100.D0
C
WRITE(6,2000) LOP, ADF, RFR, RFV
2000 FORMAT(' ** LOP=',I4,'* ADF =',E12.5, ' * Vr_VOR=',F8.3,' (%)',
1 ' * Vr_UV=',F8.3,' (%)')
C
IF ( RFR.LT.0.01D0) WRITE(6,*) ' '
IF ( RFR.LT.0.01D0)
1WRITE(6,*) ' *** Stationary condition reached at L00p =', LOP
IF ( RFR.LT.0.1D0) RETURN 1 ! .LT. 0.1 (%)
IF ( IPR.LT.0) GO TO 50
IPR = IPR + 1
IF ( IPR.LT.NLM) GO TO 999
50 CONTINUE
C
C ----------------------------------------------------------
C
CALL STMP ( U, UU, IX, IY )
C
CALL STMP ( V, VV, IX, IY )
C
CALL STMP ( PSF, PE, IX, IY )
C
CALL STND ( XX, YY, XN, YN, NDE, IX, IY, IX1, IY1 )
C
C ----------------------------------------------------------
DO 100 I = 1, NEL
C
UN( I) = ( UU( NDE( I,1)) + UU( NDE( I,2)))/ 2.D0
VN( I) = ( VV( NDE( I,1)) + VV( NDE( I,2)))/ 2.D0
C
100 CONTINUE
C ----------------------------------------------------------
C
IPR = 0
999 CONTINUE
C ----------------------------------------------------------
DO 200 I = 1, IX
DO 200 J = 1, IY
C
IF ( DABS( U( I,J)).GE.1.1D0 ) GO TO 123
200 CONTINUE
C ----------------------------------------------------------
RETURN
123 WRITE(6,2100)
2100 FORMAT(/,' *** Diverged ! ***')
RETURN 1
C
END
C **********************************************************************
C
SUBROUTINE GEOI ( WP1, PSF, WPT, PPT, EPT, IX, IY )
C
IMPLICIT REAL*8 (A-H,O-Z)

DIMENSION INDX( 5), INDY( 5), INDE( 5)
DIMENSION WP1( IX,IY), PSF( IX,IY)
C
LOGICAL*1 WPT( IX,IY,4), PPT( IX,IY,4), EPT( IX,IY,4)
C
LOGICAL*1 OMG, PSI, SLP, MOV, CNT
NAMELIST /IDAT/ OMG, PSI, SLP, MOV, CNT
C
ILP = 0
C
C ---------------------------------
DO 100 I = 1, IX
DO 100 J = 1, IY
C
WPT( I,J,1) = .TRUE.
PPT( I,J,1) = .TRUE.
C
100 CONTINUE
C ---------------------------------
DO 150 I = 1, IX
C
EPT( I,1, 2) = .TRUE.
EPT( I,IY,1) = .TRUE.
C
150 CONTINUE
C ---------------------------------
DO 300 I = 1, IY
C
EPT( 1, I,4) = .TRUE.
EPT( IX,I,3) = .TRUE.
C
300 CONTINUE
C ---------------------------------
C
C ----- Edge points ---
C
900 CONTINUE
OMG = .FALSE.
PSI = .FALSE.
SLP = .FALSE.
MOV = .FALSE.
CNT = .FALSE.
C
READ(8,IDAT,ERR=202)
C
ILP = ILP + 1
C
IF ( OMG) GO TO 205
IF ( PSI) GO TO 205
IF ( CNT) GO TO 205
C
202 WRITE(6,2100)
2100 FORMAT (' * Error in DATA')
WRITE(6,2200) OMG, PSI, SLP, MOV, CNT
2200 FORMAT(' *** OMG, PSI, SLP, MOV, CNT',5L2)
GO TO 900
C
205 CONTINUE
IICT = 1
ISLP = -10
C
C ----------------------------
DO 400 L = 1, 5
C
INDX( L) = 0
INDY( L) = 0
INDE( L) = 0
C
400 CONTINUE
C ----------------------------
C
WVAL = 0.D0
PVAL = 0.D0
READ(8,1000,ERR=225) ISLP, WVAL, PVAL,
1 ( INDX( L), INDY( L), INDE( L), L = 1, 5 )
1000 FORMAT ( I2,1X,2(F6.3,1X),5(3I3,1X))
C
IKP = 0
C
C ISLP = -1 : Boundary slopes down ISLP = 0 : Horizontal
C ISLP = 1 : Boundary slopes upward ISLP = 10 : Vertical
C
IF ( ISLP+ 1) 225, 233, 221
221 IF ( ISLP ) 225, 232, 222
222 IF ( ISLP- 1) 225, 231, 223
223 IF ( ISLP-10) 225, 230, 225
C
225 WRITE(6,2300)
2300 FORMAT(' *** Error in DATA ***')
WRITE(6,2500) ISLP, WVAL, PVAL, ( INDX( L),INDY( L),INDE( L),
1 L = 1, 5 )
2500 FORMAT(5X,I2,1X,2(F6.3,1X),5(3I3,1X))
GO TO 333
C
230 IKP = IKP + 1
231 IKP = IKP + 1
232 IKP = IKP + 1
233 IKP = IKP + 1
C
C IKP = 1 : Slopes Down IKP = 2: Horizontal
C IKP = 3 : Slopes Upper IKP = 4, Vertical
C
C ----------------------------------------------------------
DO 500 L = 1, 5
C
IUPR = INDE( L)
GO TO ( 240, 240, 240, 235 ), IKP
C
235 ILOW = INDY( L)
GO TO 245
C
240 ILOW = INDX( L)
245 I = INDX( L)
J = INDY( L)
C
C ------------------------------------------------
DO 550 IDX = ILOW, IUPR
C
IF ( I) 255, 500, 250
250 IF ( J) 255, 500, 251
251 IF ( I-IX) 252, 252, 255
252 IF ( J-IY) 260, 260, 255
C
C ----- IF I or J = 0 : No bopundary secified ---
C
C ----- Check to see this point is within array boundaries ---
C
255 WRITE (6,2400) L, IX, IY
2400 FORMAT ('* Error in DATA', I1,' IX =', I3,' IY =', I3 )
WRITE(6,2500) ISLP, WVAL, PVAL,
1 ( INDX( M),INDY( M),INDE( M), M = 1, 5 )
GO TO 500
C
260 IF ( .NOT. OMG) GO TO 265
WPT( I,J,1) = .FALSE.
WPT( I,J,2) = SLP
WPT( I,J,3) = MOV
WP1( I,J) = WVAL
C
265 IF ( .NOT. PSI) GO TO 270
PPT( I,J,1) = .FALSE.
PPT( I,J,2) = SLP
PPT( I,J,3) = MOV
PSF( I,J) = PVAL
C
270 IF ( .NOT. CNT) GO TO 280
WPT( I,J,4) =.TRUE.
PPT( I,J,4) =.TRUE.
280 GO TO ( 286, 288, 287, 285 ), IKP
C
285 J = J + 1
GO TO 550
C
286 J = J - 1
GO TO 288
C
287 J = J + 1
288 I = I + 1
550 CONTINUE
C ------------------------------------------------
C
C ----- From ILOW to INDE( L) : Assigned ---
C
500 CONTINUE
C ----------------------------------------------------------
C
333 CONTINUE
C
IF ( ILP.EQ.3) GO TO 310
C ----- Attention ---
C
GO TO 900
C
310 CONTINUE
C
C -----------------------------------------
C
CALL IFLW ( PSF, PPT, IX, IY )
C
C -----------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE IFLW ( PSF, PPT, IX, IY )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION PSF( IX, IY)
C
LOGICAL *1 PPT ( IX, IY,4)
C
COMMON /DAT2/ DT, DS, DTS, COF, UL, UV, VIS, WEI
C
N = IY-1
FLOW = ( PSF(IX,IY) - PSF(IX,1))/ (IY-1)
C
C ------------------------------------------------
DO 100 J = 2, N
C
PSF( IX,J) = PSF( IX,J-1) + FLOW
100 CONTINUE
C ------------------------------------------------
DO 200 I = 1, IX
DO 200 J = 1, IY
C
IF ( .NOT.PPT( I,J,1)) GO TO 200
PSF( I,J) = PSF( IX,J)
C
200 CONTINUE
C ------------------------------------------------
IF (.NOT.PPT( 1,1,3)) UL = 0.D0 ! Not moving
C
RETURN
END
C **********************************************************************
C
SUBROUTINE CMSB ( S, WP1, PSF, U, V, RNU, PS, WPT, EPT,
1 IX, IY, ISC, ICR )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
LOGICAL*1 WPT( IX,IY,4), EPT( IX,IY,4)
C
COMMON /DAT1/ LOP, RE, NLM, IPR, NNP, NEL
COMMON /DAT2/ DT, DS, DTS, COF, UL, UV, VIS, WEI
C
DIMENSION S( IX,IY), WP1( IX,IY), PSF( IX,IY), U( IX,IY),
1 V( IX,IY), RNU( IX,IY), PS ( IX,IY)
C
C ---------------------------------
DO 100 I = 1, IX
DO 100 J = 1, IY
C
S( I,J) = WP1( I,J)
C
100 CONTINUE
C ---------------------------------
C
IF ( ISC) 120, 120, 130
C ----------------------------------------------------------
C
120 CALL WUCL ( WP1, S, PSF, WPT, IX, IY, U )
C
C ----------------------------------------------------------
C
CALL WVCL ( WP1, S, PSF, WPT, EPT, IX, IY, V )
C
C ----------------------------------------------------------
GO TO 155
130 CONTINUE
C ----------------------------------------------------------
C
CALL WUC2 ( WP1, S, PSF, WPT, IX, IY, U )
C
C ----------------------------------------------------------
C
CALL WVC2 ( WP1, S, PSF, WPT, EPT, IX, IY, V )
C
C ----------------------------------------------------------
155 CONTINUE
C ----------------------------------------------------------
DO 200 I = 1, IX
DO 200 J = 1, IY
C
IF ( .NOT.WPT( I,J,1).AND.EPT( I,J,3)) U( I,J) = 0.D0
IF ( .NOT.WPT( I,J,1).AND.EPT( I,J,4)) U( I,J) = 0.D0
C
IF ( .NOT.WPT( I,J,1).AND.EPT( I,J,1)) V( I,J) = 0.D0
IF ( .NOT.WPT( I,J,1).AND.EPT( I,J,2)) V( I,J) = 0.D0
C
200 CONTINUE
C ---------------------------------------------------------------------
C
IF (ICR.GT.0) CALL WCOR ( WP1, S, WPT, EPT, IX, IY, U, V, RNU )
C
C ---------------------------------------------------------------------
C
IF (ICR.LE.0) CALL WCLC ( WP1, S, WPT, EPT, IX, IY )
C
C ---------------------------------------------------------------------
DO 300 I = 1, IX
DO 300 J = 1, IY
C
S ( I,J) = PS ( I,J)
PS( I,J) = PSF( I,J)
C
300 CONTINUE
C ----------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE BEGN ( WPT, PPT, EPT, WP1, PSF, S, PS, U, V, UB, VB,
1 WB, FSM, ISC, ICR, IX, IY, NBC )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
CHARACTER*8 FSM
C
LOGICAL*1 WPT( IX,IY,4), PPT( IX,IY,4), EPT( IX,IY,4)
C
COMMON /KKY / ICT
COMMON /HIST/ WMX( 5000), DMX( 5000)
COMMON /DAT1/ LOP, RE, NLM, IPR, NNP, NEL
C
DIMENSION WP1( IX,IY), PSF( IX,IY), S ( IX,IY), PS( IX,IY),
1 U ( IX,IY), V ( IX,IY), UB( IX,IY), VB( IX,IY),
2 WB ( IX,IY)
C
IPR = - 1
ICT = 0
C
NNP = IX*IY
NEL = (IX-1)* (IY-1)
C
C --------------------------------------
DO 100 J = 1, IY
DO 100 I = 1, IX
C
C ----------------------------
DO 150 K = 1,4
C
WPT( I,J,K) = .FALSE.
PPT( I,J,K) = .FALSE.
EPT( I,J,K) = .FALSE.
C
150 CONTINUE
C ----------------------------
C
WP1( I,J) = 0.D0
PSF( I,J) = 0.D0
S ( I,J) = 0.D0
PS ( I,J) = 0.D0
C
U ( I,J) = 0.D0
V ( I,J) = 0.D0
UB ( I,J) = 0.D0
VB ( I,J) = 0.D0
WB ( I,J) = 0.D0
C
100 CONTINUE
C ---------------------------------------
DO 200 I = 1, 5000
C
WMX( I) = 0.D0
DMX( I) = 0.D0
C
200 CONTINUE
C --------------------------------------
C
C ----- Input Data --------------------------
C
CALL DINP ( FSM, IX, IY, ISC, ICR, NBC )
C
C -------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE SFCL ( PSF, PS, S, WP1, PPT, ADF, IX, IY )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION PSF( IX,IY), PS( IX,IY), S( IX,IY), WP1( IX,IY)
C
LOGICAL*1 PPT( IX,IY,4), DRT
C
COMMON /KKY / ICT
COMMON /HIST/ WMX( 5000), DMX( 5000)
COMMON /DAT1/ LOP, RE, NLM, IPR, NNP, NEL
COMMON /DAT2/ DT, DS, DTS, COF, UL, UV, VIS, WEI
C
ACR = 1.5D0
EPS = 5.D-5
CFI = 0.25D0
MIN = 10
C
ITT = 0
IF ( ICT.GT.1) GO TO 50
ICT = ICT + 1
GO TO 104
C
50 CONTINUE
C ------------------------------------------------
DO 100 J = 1, IY
DO 100 I = 1, IX
C
PSF( I,J) = 2.D0 * PS( I,J) - S( I,J)
C
100 CONTINUE
C ------------------------------------------------
C
104 DRT = .FALSE.
JDM = IY + 1
105 IDFC = 0
ITT = ITT + 1
IF ( ITT.GT.100) GO TO 999
C
C ----------------------------------------------------------
DO 200 JJ = 1, IY
C
10 DRT = .NOT. DRT
IF ( DRT ) J = JJ
IF ( .NOT. DRT ) J = JDM - JJ
C
C ------------------------------------------------
DO 250 I = 1, IX
C
IF ( PPT( I,J,1)) GO TO 1
IF ( PPT( I,J,2)) GO TO 250
C
C ----- Rigid, No-SLP, Not-moving
C : Keep PSF as given also for Free-Slip ---
GO TO 250
C
1 CONTINUE
IF ( I .EQ. IX ) GO TO 2
IF ( I .EQ. 1 ) GO TO 12
C
C ----- Normal interior to fluid point ---
C
DLP = ( PSF( I+1,J) + PSF( I-1,J)+ ( PSF( I,J+1)
1 + PSF( I,J-1) + DS* DS* WP1( I,J)))/ 4.D0
ADF = DABS( DLP - PSF( I,J))
C
PSF( I,J) = PSF( I,J) + ACR* ( DLP - PSF( I,J))
C
IF ( ADF .GT.EPS ) GO TO 5
IF ( ITT .LT. MIN) GO TO 5
GO TO 250
C
5 IDFC = 1
IXX = I
IYY = J
GO TO 250
C
2 IF ( PPT( I,J,4)) GO TO 8
C
C ----- Periodic Right ---
C
PSF( I,J) = PSF( I-(IX-1), J)
GO TO 250
C
C ----- Continuative Right ---
C
8 PSF( I,J) = ( PSF( I,J+1) + PSF( I,J-1)
1 + DS* DS* WP1( I,J))* 0.5D0
GO TO 250
C
12 IF ( PPT( I,J,4)) GO TO 18
C
C ----- Periodic Left ---
C
DLP = ( PSF(I+1,J) + PSF(I-1+(IX-1),J) +
1 ( PSF(I,J+1) + PSF(I,J-1) + DS* DS* WP1(I,J)))* 0.25
C
PSF( I,J) = PSF( I,J) + ACR* ( DLP - PSF( I,J))
C *****
GO TO 250
C
C ----- Continuative Left ---
C
18 PSF( I,J) = ( PSF(I,J) + PSF(I,J) + ( PSF(I,J+1)
1 + PSF(I,J-1) + DS* DS* WP1( I,J)))* 0.25
C
250 CONTINUE
C ------------------------------------------------
IF ( DRT) GO TO 10
C
200 CONTINUE
C ----------------------------------------------------------
C
IF ( IDFC.GT.0) GO TO 105
RETURN
C
999 CONTINUE
WRITE(6,*) ' *** Attention: ** ITT.GT.100 ***'
RETURN
END
C **********************************************************************
C
SUBROUTINE STMP ( DM1, DM2, IX, IY )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
COMMON /DAT1/ LOP, RE, NLM, IPR, NNP, NEL
C
DIMENSION DM1( IX, IY), DM2( NNP)
C
NN = 0
C ------------------------------
DO 100 J = 1, IY
DO 100 I = 1, IX
C
NN = NN + 1
DM2( NN) = DM1( I,J)
C
100 CONTINUE
C ------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE STND ( XX, YY, XN, YN, NDE,IX, IY, IX1, IY1 )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
COMMON /DAT1/ LOP, RE, NLM, IPR, NNP, NEL
COMMON /DAT2/ DT, DS, DTS, COF, UL, UV, VIS, WEI
C
DIMENSION XX( NNP), YY( NNP), NDE( NEL,4), XN( IX1*IY1),
1 YN( IX1*IY1)
C
IF ( IX.EQ.IY) GO TO 900
STOP
C
900 CONTINUE
C
C -----------------------------------------------------
NN = 0
DO 100 J = 1, IY
DO 100 I = 1, IX
C
NN = NN + 1
XX( NN ) = DS* DFLOAT( I-1)
YY( NN ) = DS* DFLOAT( J-1)
C
100 CONTINUE
C -----------------------------------------------------
DO 200 J = 1, IY-1
DO 200 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
200 CONTINUE
C -----------------------------------------------------
DO 300 I = 1, NEL
C
XN( I) = ( XX( NDE( I,1)) + XX( NDE( I,2))
1 + XX( NDE( I,3)) + XX( NDE( I,4)))/ 4.D0
C
YN( I) = ( YY( NDE( I,1)) + YY( NDE( I,2))
1 + YY( NDE( I,3)) + YY( NDE( I,4)))/ 4.D0
C
300 CONTINUE
C -----------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE VRVL ( U, V, UB, VB, UMX, VMX, IX, IY )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
COMMON /HIST/ WMX( 5000), DMX( 5000)
COMMON /DAT1/ LOP, RE, NLM, IPR, NNP, NEL
C
DIMENSION U( IX,IY), V( IX,IY), UB( IX,IY), VB( IX,IY)

UMX = 0.D0
VMX = 0.D0
C
C ----------------------------------------------------------
DO 100 I = 1, IX
DO 100 J = 1, IY
C
IF ( DABS( U( I,J)).LT.0.01D0 ) GO TO 140
UMX = DMAX1( UMX, DABS( ( U( I,J) - UB( I,J))/ U( I,J)))
140 IF ( DABS( V( I,J)).LT.0.01D0 ) GO TO 150
VMX = DMAX1( VMX, DABS( ( V( I,J) - VB( I,J))/ V( I,J)))
150 CONTINUE
C
100 CONTINUE
C ----------------------------------------------------------
DMX( LOP) = DMAX1( UMX, VMX )
C
C ----------------------------------------------------------
DO 200 I = 1, IX
DO 200 J = 1, IY
C
UB( I,J) = U( I,J)
VB( I,J) = V( I,J)
C
200 CONTINUE
C ----------------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE VRVR ( WP1, WB, PSF, PL, WL1, WL2, VRV, FVV, IX,IY )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
COMMON /HIST/ WMX( 5000), DMX( 5000)
COMMON /DAT1/ LOP, RE, NLM, IPR, NNP, NEL
C
DIMENSION WP1( IX,IY), WB( IX,IY), PSF( IX,IY)
C
WL1 = 0.01D0
VRV = 0.D0
FVV = 0.D0
C ------------------------------------------------
DO 100 I = 1, IX
DO 100 J = 1, IY
C
WL1 = DMAX1( WL1, DABS( WP1( I,J)))
C
100 CONTINUE
C ------------------------------------------------
WL2 = 0.001D0* WL1
C
C ------------------------------------------------
DO 200 I = 1, IX
DO 200 J = 1, IY
C
IF ( DABS( WP1( I,J)).LT.WL2) GO TO 200
VRV = DMAX1( VRV, DABS( ( WP1( I,J) - WB( I,J))/ WP1( I,J)))
C
200 CONTINUE
C ------------------------------------------------
WMX( LOP) = VRV
FVV = 100.D0* VRV
C ------------------------------------------------
DO 300 I = 1, IX
DO 300 J = 1, IY
C
WB( I,J) = WP1( I,J)
C
300 CONTINUE
C ------------------------------------------------
C
PL = PSF( 1,1)
II = 1
JJ = 1
C
C --------------------------------------
DO 400 I = 1, IX
DO 400 J = 1, IY
C
IF ( PSF( I,J).GT.PL) THEN
PL = PSF( I,J)
II = I
JJ = J
END IF
C
400 CONTINUE
C --------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE WCLC ( WP1, W, WPT, EPT, IX, IY )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION WP1( IX,IY), W( IX,IY)
C
LOGICAL*1 WPT( IX,IY,4), EPT( IX,IY,4)
C
COMMON /DAT2/ DT, DS, DTS, COF, UL, UV, VIS, WEI
C
C ------------------------------------------------
DO 100 I = 1, IX
DO 100 J = 1, IY
C
IF ( WPT( I,J,1)) GO TO 100
WP1( I,J) = W( I,J)
C
100 CONTINUE
C ------------------------------------------------
DO 200 I = 1, IX
DO 200 J = 1, IY
C ------------------------------------------------
C
IF ( WPT( I,J,1)) GO TO 1
IF ( WPT( I,J,2)) GO TO 101
IF ( WPT( I,J,3)) GO TO 200
C
C ----- Rigid, No-Slip, Not moving ---
GO TO 101
C
1 CONTINUE
IF ( EPT( I,J,3)) GO TO 2
IF ( EPT( I,J,4)) GO TO 12
IF ( EPT( I,J,1)) GO TO 250
IF ( EPT( I,J,2)) GO TO 250
C
WP1( I,J) = WP1( I,J) + COF* ( W( I,J+1) + W( I,J-1) + ( W( I+1,J)
1 + W( I-1,J) - 4.D0* W( I,J)))
GO TO 200
C
2 IF ( WPT( I,J,4)) GO TO 8
C
C ----- Periodic Right ---
C
WP1( I,J) = WP1( I,J) + COF* ( W( I,J+1) + W( I,J-1)
1 + ( W( I+1-(IX-1),J) + W( I-1,J) - 4.D0* W( I,J)))
GO TO 200
C
C ----- Continuative Right ---
C
8 WP1( I,J) = WP1( I,J) + COF*( W( I,J+1) + W( I,J-1)-2.D0*W( I,J))
GO TO 200
C
12 IF ( WPT( I,J,4)) GO TO 18
C
C ----- Periodic Left ---
C
WP1( I,J) = WP1( I,J) + COF* ( W( I,J+1) + W( I,J-1)
1 + ( W(I+1,J) + W(I-1+(IX-1),J) - 4.D0 * W(I,J)))
GO TO 200
C
C ----- Continuative Left ---
C
18 WP1( I,J) = WP1( I,J) + COF*( W( I,J+1) + W( I,J-1)
1 + ( W( I,J) + W( I,J) - 4.D0 * W( I,J)))
GO TO 200
C
C ----- Free Slip - No Diffusion ---
C
101 GO TO 200
200 CONTINUE
C ------------------------------------------------
RETURN
C
250 CONTINUE
C
WRITE (6,2000)
2000 FORMAT (' Moving Wall : Not Defined.')
C
RETURN
END
C **********************************************************************
C
SUBROUTINE WCOR ( WP1, W, WPT, EPT, IX, IY, U, V, RNU )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION WP1(IX,IY), W(IX,IY), U(IX,IY), V(IX,IY), RNU(IX,IY)
C
LOGICAL*1 WPT( IX,IY,4), EPT( IX,IY,4)
C
COMMON /DAT2/ DT, DS, DTS, COF, UL, UV, VIS, WEI
C
C ----- Correction Method ---
C
C ------------------------------------------------
DO 100 I = 1, IX
DO 100 J = 1, IY
C
RNU( I,J) = 0.5D0* U( I,J)* V( I,J)* DT
C
100 CONTINUE
C ------------------------------------------------
DO 200 I = 1, IX
DO 200 J = 1, IY
C
IF ( WPT( I,J,1)) GO TO 200
WP1( I,J) = W( I,J)
C
200 CONTINUE
C ------------------------------------------------
DO 300 I = 1, IX
DO 300 J = 1, IY
C
AVIS = ( VIS + RNU( I,J))* DT/( DS* DS)
IF ( WPT( I,J,1)) GO TO 10
IF ( WPT( I,J,2)) GO TO 300
IF ( WPT( I,J,3)) GO TO 300
C
C ----- Rigid, No-Slip, Not moving ---
GO TO 300

10 CONTINUE
IF ( EPT( I,J,3)) GO TO 20
IF ( EPT( I,J,4)) GO TO 12
IF ( EPT( I,J,1)) GO TO 900
IF ( EPT( I,J,2)) GO TO 900
C
WP1( I,J) = WP1( I,J) + AVIS* ( W(I,J+1) + W(I,J-1)
1 + ( W(I+1,J) + W(I-1,J) - 4.D0* W(I,J)))
GO TO 300
C
20 IF ( WPT( I,J,4)) GO TO 30
C
C ----- Periodic Right ---
C
WP1( I,J) = WP1( I,J) + AVIS* ( W(I,J+1) + W(I,J-1)
1 + ( W(I+1-(IX-1),J) + W(I-1,J) - 4.D0* W( I,J)))
GO TO 300
C
C ----- Continuative Right ---
C
30 WP1( I,J) = WP1( I,J)+ AVIS* ( W(I,J+1) + W(I,J-1) - 2.D0*W(I,J))
GO TO 300
C
12 IF ( WPT( I,J,4)) GO TO 40
C
C ----- Periodic Left ---
C
WP1( I,J) = WP1( I,J) + AVIS* ( W(I,J+1) + W(I,J-1)
1 + ( W(I+1,J) + W(I-1+(IX-1),J) - 4.D0* W(I,J)))
GO TO 300
C
C ----- Continuative Left ---
C
40 WP1( I,J) = WP1( I,J) + AVIS* ( W(I,J+1) + W(I,J-1)
1 + ( W(I,J) + W(I,J) - 4.D0* W(I,J)))
GO TO 300
C
C ----- Free Slip - No Diffusion ---
C
300 CONTINUE
C ------------------------------------------------
RETURN
C
900 CONTINUE
RETURN
END
C **********************************************************************
C
SUBROUTINE VORB0 ( WP1, PSF, WPT, IX, IY )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION WP1( IX,IY), PSF( IX,IY)
C
LOGICAL*1 WPT( IX,IY,4)
C
COMMON /DAT2/ DT, DS, DTS, COF, UL, UV, VIS, WEI
C
C ----- Vorticity Boundary Conditions ---
C
C Attention : Moving condition
C
C ----- No-Slip Surfaces updated ---
C
DSQ = DS* DS
C ==================
C
C ------------------------------------------------
DO 100 I = 1, IX
DO 100 J = 1, IY
C ------------------------------------------------
C
IF ( WPT( I,J,1)) GO TO 100 ! Fluid
IF ( WPT( I,J,2)) GO TO 50 ! Slip
IF ( WPT( I,J,3)) GO TO 200 ! Moving
IF ( J.EQ.IY) GO TO 10
IF ( J.EQ. 1) GO TO 20
IF ( I.EQ.IX) GO TO 30
IF ( I.EQ.1 ) GO TO 40
C
C ----- No slip surfaces ---
C
IF ( .NOT. WPT( I+1,J,1)) GO TO 1
WA = ( PSF( I,J) - PSF( I+1,J))/ DSQ
C
13 IF ( .NOT. WPT( I-1,J,1)) GO TO 3
WB = ( PSF( I,J) - PSF( I-1,J))/ DSQ
C
12 IF ( .NOT. WPT( I,J+1,1)) GO TO 2
WC = ( PSF( I,J) - PSF( I,J+1))/ DSQ
C
14 IF ( .NOT. WPT( I,J-1,1)) GO TO 4
WD = ( PSF( I,J) - PSF( I,J-1))/ DSQ
C
110 WP1( I,J) = WA + WB + WC + WD
GO TO 100
C
1 WA = ( PSF( I,J) - PSF( I-1,J))/ DSQ
GO TO 13
C
3 WB = ( PSF( I,J) - PSF( I+1,J))/ DSQ
GO TO 12
C
2 WC = ( PSF( I,J) - PSF( I,J-1))/ DSQ
GO TO 14
C
4 WD = ( PSF( I,J) - PSF( I,J+1))/ DSQ
GO TO 110
C
C ----- No slip surfaces at edges ---
C
50 WP1( I,J) = 0.D0
GO TO 100
C
C ----- Attention : BC : Wood's Approximation ------------------------
C
10 CONTINUE
WP1( I,J) = 3.D0* ( PSF( I,J) - PSF( I,J-1) - DS )/ DSQ
1 - 0.5D0* WP1( I,J-1)
C
IF ( I.EQ.1.OR.I.EQ.IX ) WP1( I,J) = 0.D0
GO TO 100
C
20 CONTINUE
WP1( I,J) = 3.D0*( PSF(I,J) - PSF(I,J+1))/ DSQ - 0.5D0* WP1(I,J+1)
GO TO 100
C
30 CONTINUE
WP1( I,J) = 3.D0*( PSF(I,J) - PSF(I-1,J))/ DSQ - 0.5D0* WP1(I-1,J)
GO TO 100
C
40 CONTINUE
WP1( I,J) = 3.D0*( PSF(I,J) - PSF(I+1,J))/ DSQ - 0.5D0* WP1(I+1,J)
GO TO 100
C
200 CONTINUE
C
IF ( J .EQ.IY) GO TO 10
IF ( J .EQ. 1) GO TO 20
C
100 CONTINUE
C ------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE VORB1 ( WP1, PSF, WPT, IX, IY )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION WP1( IX,IY), PSF( IX,IY)
LOGICAL*1 WPT( IX,IY,4)
C
COMMON /DAT2/ DT, DS, DTS, COF, UL, UV, VIS, WEI
C
C ----- Vorticity Boundary Conditions --- Ghia's Approximation ---
C
C Attention : Moving condition
C
C ----- No-Slip Surfaces updated ---
C
DSQ = DS* DS
C ==================
C
C ------------------------------------------------
DO 100 I = 1, IX
DO 100 J = 1, IY
C ------------------------------------------------
C
IF ( WPT( I,J,1)) GO TO 100 ! Fluid
IF ( WPT( I,J,2)) GO TO 50 ! Slip
IF ( WPT( I,J,3)) GO TO 200 ! Moving
IF ( J.EQ.IY) GO TO 10
IF ( J.EQ. 1) GO TO 20
IF ( I.EQ.IX) GO TO 30
IF ( I.EQ.1 ) GO TO 40
C
C ----- No slip surfaces ---
C
IF ( .NOT. WPT( I+1,J,1)) GO TO 1
WA = ( PSF( I,J) - PSF( I+1,J))/ DSQ
C
13 IF ( .NOT. WPT( I-1,J,1)) GO TO 3
WB = ( PSF( I,J) - PSF( I-1,J))/ DSQ
C
12 IF ( .NOT. WPT( I,J+1,1)) GO TO 2
WC = ( PSF( I,J) - PSF( I,J+1))/ DSQ
C
14 IF ( .NOT. WPT( I,J-1,1)) GO TO 4
WD = ( PSF( I,J) - PSF( I,J-1))/ DSQ
C
110 WP1( I,J) = WA + WB + WC + WD
GO TO 100
C
1 WA = ( PSF( I,J) - PSF( I-1,J))/ DSQ
GO TO 13
C
3 WB = ( PSF( I,J) - PSF( I+1,J))/ DSQ
GO TO 12
C
2 WC = ( PSF( I,J) - PSF( I,J-1))/ DSQ
GO TO 14
C
4 WD = ( PSF( I,J) - PSF( I,J+1))/ DSQ
GO TO 110
C
C ----- No slip surfaces at edges ---
C
50 WP1( I,J) = 0.D0
GO TO 100
C
C ----- Attention : BC : Ghia's Approximation ---
C
10 CONTINUE
WP1( I,J) = ( PSF( I,J-2) - 8.D0* PSF( I,J-1) + 7.D0* PSF( I,J) )
1 / ( 2.D0* DSQ ) - 3.D0/ DS
C
IF ( I.EQ.1.OR.I.EQ.IX ) WP1( I,J) = 0.D0
GO TO 100
C
20 CONTINUE
WP1( I,J) = ( 7.D0* PSF(I,J) - 8.D0* PSF(I,J+1) + PSF(I,J+2) )
1 /( 2.D0* DSQ )
GO TO 100
C
30 CONTINUE
WP1( I,J) = ( 7.D0* PSF(I,J) - 8.D0* PSF(I-1,J) + PSF(I-2,J) )
1 /( 2.D0* DSQ )
GO TO 100
C
40 CONTINUE
WP1( I,J) = ( 7.D0* PSF(I,J) - 8.0D0* PSF(I+1,J) + PSF(I+2,J) )
1 /( 2.D0* DSQ )
GO TO 100
C
200 CONTINUE
IF ( J .EQ.IY) GO TO 10
IF ( J .EQ. 1) GO TO 20
C
100 CONTINUE
C ------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE WUC2 ( WP1, S , PSF, WPT, IX, IY, UF )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
C ----- Attention --- "S" ---
C
DIMENSION WP1( IX,IY), S( IX,IY), PSF( IX,IY), UF( IX,IY)
C
LOGICAL*1 WPT( IX,IY,4)
C
COMMON /DAT2/ DT, DS, DTS, COF, UL, UV, VIS, WEI
C
C ----------------------------------------------------------
DO 100 I = 1, IX
DO 100 J = 1, IY
C
IF ( I.EQ. IX) GO TO 100
IF ( J.EQ. IY) GO TO 120
IF ( J.EQ. 1 ) GO TO 140
C
UF( I,J) = ( PSF( I,J+1) - PSF( I,J-1)
1 + PSF( I+1,J+1)- PSF( I+1,J-1))/( 4.D0* DS )
GO TO 160
C
120 UF( I,J) = UV
C ===== Upper wall ===
C
GO TO 160

140 UF( I,J) = UL
160 ALF = UF( I,J)* DTS
C
IF ( ALF.LT. 0.D0) GO TO 130
IF ( .NOT. WPT( I,J,1)) GO TO 150
IF ( WPT( I,J,4)) GO TO 150
C
A1 = 0.5D0* ALF
A2 = A1* ALF
B2 = A2 - 2.D0* A1
TRP = ( S( I-1,J) + S( I,J))* ( 1.D0 - WEI)
TRM = ( S( I-1,J) - S( I,J))* ( 1.D0 - WEI)
C
TPL = ( S( I,J) + S( I+1,J))* WEI
TML = ( S( I,J) - S( I+1,J))* WEI
C
TF = A1* TPL + A1* TRP + A2* TML + B2* TRM
WP1( I+1,J) = WP1( I+1,J) + TF
WP1( I, J) = WP1( I,J ) - TF
GO TO 100
C
130 IF ( WPT( I,J,4)) GO TO 170
C
K = I + 1
A1 = 0.5D0* ALF
C =====================
A2 = A1* ALF
B2 = A2 + 2.D0* A1
TRP = ( S( K,J) + S( K+1,J))* ( 1.D0 - WEI)
TRM = ( S( K,J) - S( K+1,J))* ( 1.D0 - WEI)
TPL = ( S( K-1,J) + S( K, J) )* WEI
TML = ( S( K-1,J) - S( K, J) )* WEI
C
TF = A1* TPL + A1* TRP + A2* TML + B2* TRM
C
WP1( K, J) = WP1( K,J) + TF
WP1( K-1,J) = WP1( K-1,J) - TF
GO TO 100
C
150 TF = ALF* S( I,J)
GO TO 180
C
170 TF = ALF* S( I+1,J)
180 WP1( I,J) = WP1( I,J) - TF
WP1( I+1,J) = WP1( I+1,J) + TF
IF ( WPT( I,J,4)) WP1( I,J) = WP1( I,J) + TF
C
100 CONTINUE
C ----------------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE WUCL ( WP1, S , PSF, WPT, IX, IY, UF )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION WP1( IX,IY), S( IX,IY), PSF( IX,IY), UF( IX,IY)
C
LOGICAL*1 WPT ( IX,IY,4)
C
COMMON /DAT2/ DT, DS, DTS, COF, UL, UV, VIS, WEI
C
A11 = 7.D0/ 12.D0
A12 = 5.D0/ 8.D0
A13 = 1.D0/ 2.D0 - A11
A14 = 1.D0/ 2.D0 - A12
C
A21 = A13
A22 = A14/ 3.D0
A23 = - A21
A24 = - A22
C
C ----------------------------------------------------------
DO 100 I = 1, IX
DO 100 J = 1, IY
C
IF ( I .EQ. IX) GO TO 100
IF ( J .EQ. IY) GO TO 10
IF ( J .EQ. 1) GO TO 20
C
UF( I,J) = ( PSF( I,J+1) - PSF( I,J-1)
1 + PSF( I+1,J+1) - PSF( I+1,J-1))/( 4.D0* DS)
C ============================================================
GO TO 30
C
C ----- Attention *** Moving Wall ---
C
10 UF( I,J) = UV
C ===================================
C
C ----- Upper wall ---
C
GO TO 30
C
20 UF( I,J) = UL
C ====================
C
30 ALF = UF( I,J)* DTS
IF ( ALF .LT. 0.D0) GO TO 250
IF ( .NOT. WPT( I,J,1)) GO TO 600
IF ( WPT( I,J,4)) GO TO 600
IF ( WPT( I+1,J,4)) GO TO 600
IF ( I.EQ.1) GO TO 400
C
IF ( .NOT. WPT( I-1,J,1)) GO TO 200
IF ( WPT( I-1,J,4)) GO TO 200
IF ( .NOT. WPT( I+1,J,1)) GO TO 200
IF ( WPT( I+1,J,4)) GO TO 200
C
400 CONTINUE
AA = ALF*ALF
AAA = ALF*AA
AAAA = ALF*AAA
C
A1 = A11*ALF
A2 = A12*AA
A3 = A13*AAA
A4 = A14*AAAA
C
B1 = A21*ALF
B2 = A22*AA
B3 = A23*AAA
B4 = A24*AAAA
C
B = ALF - 1.D0
BB = B*B
BBB = B*BB
BBBB = B*BBB
C
C1 = A11*B + A11
C2 = A12*BB - A12
C3 = A13*BBB + A13
C4 = A14*BBBB - A14
C
D1 = A21*B
D2 = A22*BB
D3 = A23*BBB
D4 = A24*BBBB
C
IF ( I .EQ. 1) GO TO 120
IF ( I-1.EQ. 1) GO TO 140
IF ( I+1.EQ.IX) GO TO 160
C
FDP = ( S( I-2,J) + S( I+1,J))* ( 1.D0 - WEI)
FDM = ( S( I-2,J) - S( I+1,J))* ( 1.D0 - WEI)
185 FBP = ( S( I-1,J) + S( I+2,J))* WEI
FBM = ( S( I-1,J) - S( I+2,J))* WEI
C
190 FCP = ( S( I-1,J) + S( I,J)) * ( 1.D0 - WEI)
FCM = ( S( I-1,J) - S( I,J)) * ( 1.D0 - WEI)
195 FAP = ( S( I, J) + S( I+1,J))* WEI
FAM = ( S( I, J) - S( I+1,J))* WEI
C
TF = A1* FAP + A2* FAM + A3* FAP + A4* FAM
1 + B1* FBP + B2* FBM + B3* FBP + B4* FBM
2 + C1* FCP + C2* FCM + C3* FCP + C4* FCM
3 + D1* FDP + D2* FDM + D3* FDP + D4* FDM
C
WP1( I,J) = WP1( I,J) - TF
WP1( I+1,J) = WP1( I+1,J) + TF
C
IF ( I+1.LT.IX) GO TO 180
IF ( .NOT. WPT( I+1,J,1)) GO TO 100
C
C ----- Periodic Right ---
C
WP1( I+1-(IX-1),J) = WP1( I+1-(IX-1),J) + TF
GO TO 100
C
180 IF ( I.GT.1) GO TO 100
IF ( .NOT.WPT( I,J,1)) GO TO 100
C
C ----- Periodic Left ---
C
WP1( I+(IX-1),J) = WP1( I+(IX-1),J) - TF
GO TO 100
C
120 FDP = ( S(I-2+(IX-1),J) + S( I+1,J))* ( 1.D0 - WEI)
FDM = ( S(I-2+(IX-1),J) - S( I+1,J))* ( 1.D0 - WEI)
FBP = ( S(I-1+(IX-1),J) + S( I+2,J))* WEI
C
FBM = ( S(I-1+(IX-1),J) - S( I+2,J))* WEI
FCP = ( S(I-1+(IX-1),J) + S( I,J)) * ( 1.D0 - WEI)
FCM = ( S(I-1+(IX-1),J) - S( I,J)) * ( 1.D0 - WEI)
GO TO 195
C
140 FDP = ( S(I-2+(IX-1), J) + S(I+1,J))* ( 1.D0 - WEI)
FDM = ( S(I-2+(IX-1), J) - S(I+1,J))* ( 1.D0 - WEI)
GO TO 185
C
160 FBP = ( S( I-1,J) + S( I+2-(IX-1), J))* WEI
FBM = ( S( I-1,J) - S( I+2-(IX-1), J))* WEI
FDP = ( S( I-2,J) + S( I+1,J))* ( 1.D0 - WEI)
FDM = ( S( I-2,J) - S( I+1,J))* ( 1.D0 - WEI)
GO TO 190
C
250 IF ( .NOT. WPT( I+1,J,1)) GO TO 601
IF ( WPT( I+1,J,4)) GO TO 601
IF ( WPT( I,J,4) ) GO TO 601
C
K = I+1
IF ( K .EQ.IX) GO TO 401
IF ( .NOT. WPT( K+1,J,1)) GO TO 201
IF ( WPT( K+1,J,4)) GO TO 201
IF ( .NOT. WPT( K-1,J,1)) GO TO 201
IF ( WPT( K-1,J,4)) GO TO 201
C
401 AA = ALF*ALF
AAA = ALF*AA
AAAA = ALF*AAA
C
A1 = A11*ALF
A2 = A12*AA
A3 = A13*AAA
A4 = A14*AAAA
C
B1 = A21*ALF
B2 = A22*AA
B3 = A23*AAA
B4 = A24*AAAA
C
B = ALF + 1.D0
BB = B*B
BBB = B*BB
BBBB = B*BBB
C
C1 = A11*B - A11
C2 = A12*BB - A12
C3 = A13*BBB - A13
C4 = A14*BBBB - A14
C
D1 = A21*B
D2 = A22*BB
D3 = A23*BBB
D4 = A24*BBBB
C
IF ( K. EQ. IX ) GO TO 314
IF ( K+1 .EQ.IX) GO TO 315
IF ( K-1 .EQ.1 ) GO TO 319
C
FDP = ( S( K-1,J) + S( K + 2,J))* ( 1.D0 - WEI)
FDM = ( S( K-1,J) - S( K + 2,J))* ( 1.D0 - WEI)
318 FBP = ( S( K-2,J) + S( K+1,J)) * WEI
FBM = ( S( K-2,J) - S( K+1,J)) * WEI
C
317 FCP = ( S( K,J) + S( K+1, J))* ( 1.D0 - WEI)
FCM = ( S( K,J) - S( K+1, J))* ( 1.D0 - WEI)
316 FAP = ( S( K-1,J) + S( K, J))* WEI
FAM = ( S( K-1,J) - S( K, J))* WEI
C
TF = A1* FAP + A2* FAM + A3* FAP + A4* FAM
1 + B1* FBP + B2* FBM + B3* FBP + B4* FBM
2 + C1* FCP + C2* FCM + C3* FCP + C4* FCM
3 + D1* FDP + D2* FDM + D3* FDP + D4* FDM

WP1( K-1,J) = WP1( K-1,J) - TF
WP1( K, J) = WP1( K, J) + TF
C
IF ( K .LT. IX ) GO TO 312
IF ( .NOT. WPT( K,J,1)) GO TO 100
C
C ----- Periodic Right ---
C
WP1( K-(IX-1),J) = WP1( K-(IX-1),J) + TF
GO TO 100
C
312 CONTINUE
C
IF ( K-1.GT. 1) GO TO 100
IF ( .NOT. WPT( K-1,J,1)) GO TO 100
C
C ----- Periodic Left ---
C
WP1( K-1+(IX-1),J) = WP1( K-1+(IX-1),J) - TF
GO TO 100
C
314 FDP = ( S( K-1,J) + S( K+2-(IX-1), J))* ( 1.D0 - WEI)
FDM = ( S( K-1,J) - S( K+2-(IX-1), J))* ( 1.D0 - WEI)
FBP = ( S( K-2,J) + S( K+1-(IX-1), J))* WEI
C
FBM = ( S( K-2,J) - S( K+1-( IX-1), J))* WEI
FCP = ( S( K, J ) + S( K+1-( IX-1),J))* ( 1.D0 - WEI)
FCM = ( S( K, J ) - S( K+1-( IX-1),J))* ( 1.D0 - WEI)
GO TO 316
C
315 FDP = ( S( K-1,J) + S( K+2-(IX-1), J))* ( 1.D0 - WEI)
FDM = ( S( K-1,J) - S( K+2-(IX-1), J))* ( 1.D0 - WEI)
GO TO 318
C
319 FBP = ( S( K-2+(IX-1),J) + S( K+1,J))* WEI
FBM = ( S( K-2+(IX-1),J) - S( K+1,J))* WEI
FDP = ( S( K-1,J) + S( K+2,J))* ( 1.D0 - WEI)
FDM = ( S( K-1,J) - S( K+2,J))* ( 1.D0 - WEI)
GO TO 317
C
600 TF = ALF* S( I,J)
GO TO 602
C
601 TF = ALF* S( I+1,J)
602 WP1( I,J) = WP1( I,J) - TF
WP1( I+1,J) = WP1( I+1,J) + TF
C
IF ( WPT( I,J,4)) WP1( I,J) = WP1( I,J) + TF
IF ( WPT( I+1,J,4))
1 WP1( I+1,J) = S( I+1,J) + ALF* ( S( I,J) - S( I+1,J))
GO TO 100
C
200 A1 = 0.5D0 * ALF
A2 = A1 * ALF
B2 = A2 - 2.D0 * A1
C
TRP = ( S( I-1,J) + S( I,J)) * ( 1.D0 - WEI)
TRM = ( S( I-1,J) - S( I,J)) * ( 1.D0 - WEI)
TPL = ( S( I,J ) + S( I+1,J))* WEI
TML = ( S( I,J ) - S( I+1,J))* WEI
C
TF = A1* TPL + A1* TRP + A2* TML + B2* TRM
C
WP1( I+1,J) = WP1( I+1,J) + TF
WP1( I,J) = WP1( I,J) - TF
GO TO 100
C
201 A1 = 0.5D0* ALF
A2 = A1* ALF
B2 = A2 + 2.D0* A1
C
TRP = ( S( K, J) + S( K+1,J))* ( 1.D0 - WEI)
TRM = ( S( K, J) - S( K+1,J))* ( 1.D0 - WEI)
TPL = ( S( K-1,J) + S( K, J)) * WEI
TML = ( S( K-1,J) - S( K, J)) * WEI

TF = A1* TPL + A1* TRP + A2* TML + B2* TRM

WP1( K, J) = WP1( K,J) + TF
WP1( K-1,J) = WP1( K-1,J) - TF
C
100 CONTINUE
C **********************************************
C
RETURN
END
C *********************************************************************
C
SUBROUTINE WVC2 ( WP1, S , PSF, WPT, EPT, IX, IY, VF )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION WP1( IX,IY), S( IX,IY), PSF( IX,IY), VF( IX,IY)

LOGICAL*1 WPT( IX,IY,4), EPT( IX,IY,4)

COMMON /DAT2/ DT, DS, DTS, COF, UL, UV, VIS, WEI
C
C ----------------------------------------------------------
C
DO 100 I = 1, IX
DO 100 J = 1, IY
C
IF ( EPT( I,J,1)) GO TO 100
IF ( EPT( I,J,3)) GO TO 100
IF ( EPT( I,J,4)) GO TO 100
C
VF( I,J) = ( PSF( I-1,J) - PSF( I+1,J)
1 + PSF( I-1,J+1) - PSF( I+1,J+1))/ ( 4.D0* DS )
C
BET = VF( I,J)* DTS
C ========================
C
IF ( BET .LT. 0.D0) GO TO 201
IF ( .NOT. WPT( I,J,1)) GO TO 600
IF ( EPT( I,J,2)) GO TO 200
C
200 A1 = 0.5D0 * BET
A2 = A1 * BET
B2 = A2 - 2.D0 * A1
TRP = ( S( I,J-1) + S( I,J)) * ( 1.D0 - WEI)
TRM = ( S( I,J-1) - S( I,J)) * ( 1.D0 - WEI)
TPL = ( S( I,J ) + S( I,J+1))* WEI
TML = ( S( I,J ) - S( I,J+1))* WEI
C
TF = A1* TPL + A1* TRP + A2* TML + B2* TRM
C
WP1( I,J+1) = WP1( I,J+1) + TF
WP1( I,J) = WP1( I,J) - TF
GO TO 100
C
201 IF ( .NOT.WPT( I,J+1,1)) GO TO 601
K = J + 1
A1 = 0.5D0* BET
A2 = A1* BET
B2 = A2 + 2.D0 * A1
TRP = ( S( I,K ) + S( I,K+1))* ( 1.D0 - WEI)
TRM = ( S( I,K ) - S( I,K+1))* ( 1.D0 - WEI)
TPL = ( S( I,K-1) + S( I,K ))* WEI
TML = ( S( I,K-1) - S( I,K ))* WEI
C
TF = A1* TPL + A1* TRP + A2* TML + B2* TRM
C
WP1( I,K) = WP1( I,K) + TF
WP1( I,K-1) = WP1( I,K-1) - TF
GO TO 100
C
600 TF = BET* S( I,J)
C
GO TO 602
C
601 TF = BET* S( I,J+1)
C
602 WP1( I,J) = WP1( I,J) - TF
WP1( I,J+1) = WP1( I,J+1) + TF
C
100 CONTINUE
C ----------------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE WVCL ( WP1, S , PSF, WPT, EPT, IX, IY, VF )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION WP1( IX,IY), S( IX,IY), PSF( IX,IY), VF( IX,IY)
C
LOGICAL*1 WPT( IX,IY,4), EPT( IX,IY,4)
C
COMMON /DAT2/ DT, DS, DTS, COF, UL, UV, VIS, WEI
C
DATA A11 /0.58333333333333D0/, A12 / 0.625D0 /,
1 A13 /-0.08333333333333D0/, A14 /-0.125D0 /,
2 A21 /-0.08333333333333D0/, A22 /-0.0416666666666667D0/,
2 A23 / 0.08333333333333D0/, A24 / 0.0416666666666667D0/
C
C A11 = 7/12, A12 = 5/8, A13 = 1/2-A11, A14 = 1/2-A12,
C
C A21 = A13, A22 = A14/3, A23 = - A21, A24 = - A22
C
C ----------------------------------------------------------
DO 100 I = 1, IX
DO 100 J = 1, IY
C ----------------------------------------------------------
C
IF ( EPT( I,J,1)) GO TO 100
IF ( EPT( I,J,3)) GO TO 502
IF ( EPT( I,J,4)) GO TO 512
C ----------------------------------------------------------
C
VF( I,J)= ( PSF( I-1,J) - PSF( I+1,J) + PSF( I-1,J+1)
1 - PSF( I+1,J+1))/ ( 4.D0* DS)
C
C ----------------------------------------------------------
C
505 BET = VF( I,J)* DTS
C
IF ( BET .LT. 0.D0) GO TO 250
IF ( .NOT. WPT( I,J,1) ) GO TO 600
IF ( EPT( I,J,2)) GO TO 400
IF ( .NOT. WPT( I,J-1,1)) GO TO 200
IF ( .NOT. WPT( I,J+1,1)) GO TO 200
C
400 AA = BET*BET
AAA = BET*AA
AAAA = BET*AAA
C
A1 = A11*BET
A2 = A12*AA
A3 = A13*AAA
A4 = A14*AAAA
C
B1 = A21*BET
B2 = A22*AA
B3 = A23*AAA
B4 = A24*AAAA
C
B = BET - 1.D0
BB = B*B
BBB = B*BB
BBBB = B*BBB
C
C1 = A11*B + A11
C2 = A12*BB - A12
C3 = A13*BBB + A13
C4 = A14*BBBB - A14
C
D1 = A21*B
D2 = A22*BB
D3 = A23*BBB
D4 = A24*BBBB
C
FDP = ( S( I,J-2) + S( I,J+1))* ( 1.D0 - WEI)
FDM = ( S( I,J-2) - S( I,J+1))* ( 1.D0 - WEI)
FBP = ( S( I,J-1) + S( I,J+2))* WEI
FBM = ( S( I,J-1) - S( I,J+2))* WEI
C
FCP = ( S( I,J-1) + S( I,J)) * ( 1.D0 - WEI)
FCM = ( S( I,J-1) - S( I,J)) * ( 1.D0 - WEI)
FAP = ( S( I,J) + S( I,J+1))* WEI
FAM = ( S( I,J) - S( I,J+1))* WEI
C
TF = A1* FAP + A2* FAM + A3* FAP + A4* FAM
1 + B1* FBP + B2* FBM + B3* FBP + B4* FBM
2 + C1* FCP + C2* FCM + C3* FCP + C4* FCM
3 + D1* FDP + D2* FDM + D3* FDP + D4* FDM
C
WP1( I,J) = WP1( I,J) - TF
WP1( I,J+1) = WP1( I,J+1) + TF
GO TO 100
C
502 IF ( WPT( I,J+1,4) .OR. WPT( I,J,4)) GO TO 202
IF ( .NOT.WPT( I,J+1,1).AND. .NOT.WPT( I,J,1)) GO TO 100
C
C ----- Periodic Right ---
C
VF( I,J) = ( PSF( I-1,J) - PSF( I+1-(IX-1),J) + PSF( I-1,J+1)
1 - PSF( I+1-(IX-1), J+1))/ ( 4.D0 * DS )
GO TO 505
C
512 IF ( WPT( I,J+1,4) .OR. WPT( I,J,4)) GO TO 100
IF ( .NOT.WPT( I,J+1,1).AND. .NOT. WPT( I,J,1))
1 GO TO 100
C
C ----- Periodic Left ---
C
VF( I,J) = ( PSF( I-1+(IX-1),J) - PSF( I+1,J)
1 + PSF( I-1+(IX-1),J+1) - PSF( I+1,J+1))/( 4.D0* DS)
C
GO TO 505
C
202 VF( I,J) = ( PSF( I-1,J) - PSF( I,J) + PSF( I-1,J+1)
1 - PSF( I,J+1))/( 2.D0* DS)
C
BET = VF( I,J)* DTS
IF ( BET .LT. 0.D0) GO TO 601
GO TO 600
C
250 IF ( .NOT. WPT( I,J+1,1)) GO TO 601
C
K = J + 1
IF ( EPT( I,K,1)) GO TO 401
IF ( .NOT.WPT( I,K+1,1)) GO TO 201
IF ( .NOT.WPT( I,K-1,1)) GO TO 201
C
401 AA = BET*BET
AAA = BET*AA
AAAA = BET*AAA
C
A1 = A11*BET
A2 = A12*AA
A3 = A13*AAA
A4 = A14*AAAA
C
B1 = A21*BET
B2 = A22*AA
B3 = A23*AAA
B4 = A24*AAAA
C
B = BET + 1.D0
BB = B*B
BBB = B*BB
BBBB = B*BBB
C
C1 = A11*B - A11
C2 = A12*BB - A12
C3 = A13*BBB - A13
C4 = A14*BBBB - A14
C
D1 = A21*B
D2 = A22*BB
D3 = A23*BBB
D4 = A24*BBBB
C
FDP = ( S( I,K-1) + S( I,K+2))* ( 1.D0 - WEI)
FDM = ( S( I,K-1) - S( I,K+2))* ( 1.D0 - WEI)
FBP = ( S( I,K-2) + S( I,K+1))* WEI
FBM = ( S( I,K-2) - S( I,K+1))* WEI
C
FCP = ( S( I,K) + S( I,K+1))* ( 1.D0 - WEI)
FCM = ( S( I,K) - S( I,K+1))* ( 1.D0 - WEI)
FAP = ( S( I,K-1) + S( I,K)) * WEI
FAM = ( S( I,K-1) - S( I,K)) * WEI
C
TF = A1* FAP + A2* FAM + A3* FAP + A4* FAM
1 + B1* FBP + B2* FBM + B3* FBP + B4* FBM
2 + C1* FCP + C2* FCM + C3* FCP + C4* FCM
3 + D1* FDP + D2* FDM + D3* FDP + D4* FDM
C
WP1( I,K-1) = WP1( I,K-1) - TF
WP1( I,K ) = WP1( I,K) + TF
GO TO 100
C
600 TF = BET* S( I,J)
GO TO 602
C
601 TF = BET* S( I,J+1)
C
602 WP1( I,J) = WP1( I,J) - TF
WP1( I,J+1) = WP1( I,J+1) + TF
C
GO TO 100
C
200 A1 = 0.5D0* BET
A2 = A1* BET
B2 = A2 - 2.D0* A1
C
TRP = ( S( I,J-1) + S( I,J)) * ( 1.D0 - WEI)
TRM = ( S( I,J-1) - S( I,J)) * ( 1.D0 - WEI)
TPL = ( S( I,J) + S( I,J+1))* WEI
TML = ( S( I,J) - S( I,J+1))* WEI
C
TF = A1* TPL + A1* TRP + A2* TML + B2* TRM
C
WP1( I,J+1) = WP1( I,J+1) + TF
WP1( I,J) = WP1( I,J) - TF
GO TO 100
C
201 CONTINUE
C
A1 = 0.5D0* BET
A2 = A1* BET
B2 = A2 + 2.D0* A1
C
TRP = ( S( I,K) + S( I,K+1))* ( 1.D0 - WEI)
TRM = ( S( I,K) - S( I,K+1))* ( 1.D0 - WEI)
TPL = ( S( I,K-1) + S( I,K)) * WEI
TML = ( S( I,K-1) - S( I,K)) * WEI
C
TF = A1* TPL + A1* TRP + A2* TML + B2* TRM
C
WP1( I,K) = WP1( I,K) + TF
WP1( I,K-1) = WP1( I,K-1) - TF
C
100 CONTINUE
C
RETURN
END
C **********************************************************************