戻る

PROGRAM MAIN

C *********************************************************** C
C C
C F 2 D - G F D C
C C
C ( Fourth or Second-ordered FDM ) C
C C
C ( V.5.3 ) C
C C
C Copyright : Yasuhiro MATSUDA C
C C
C ( Seven Flow Problems / Mesh of equal length ) C
C C
C V-Shape / Square / Circle / Backstep( 1) C
C C
C Backstep( 2) / Cavity( S) / Cavity( L) C
C C
C *********************************************************** C
C
IMPLICIT REAL*8 (A-H,O-Z )
C
INTEGER X, Y
C
PARAMETER ( X = 60, Y = 60, ! Max. Computational Area
C
1 NNP= X* Y, NEL= ( X-1)*( Y-1)) ! Temp. Max. Definition
C
LOGICAL*1 WPT( X,Y,4), PPT( X,Y,4), EPT( X,Y,4), FLV( NEL)
C
DIMENSION PS1( X,Y), WP1( X,Y), PS ( X,Y), WM1( X,Y),
1 U ( X,Y), V ( X,Y), WS ( X,Y), XN ( NNP),
2 YN ( NNP), PS2( NNP), WP2( NNP),
3 NDE( NEL,4)
C
DIMENSION XYZ( 2,NNP), UU( NNP), VV( NNP), UE( NEL),
1 VE ( NEL), XE( NEL), YE( NEL)
C
COMMON /HST/ WMX( 5000), PMA( 5000)
COMMON /CNS/ NPB, WEI, BV, RV, KCR, OMW
COMMON /IDT/ DT, LMT, UL, VU, RD, VIS, REN, LOP
COMMON /UDT/ DX, DY, DX2, DY2, RS, RS2, CF, DTX, DTY
C
CHARACTER*14 PRB
C
OPEN (12, FILE ='F2D_GFD_VRW.DAT')
OPEN (15, FILE ='F2D_GFD_PLT.DAT' )
C
C ----- Attention : IX & IY ---
C
IX = X
IY = Y
C
CALL INPT ( WP1, PS1, S, PS, WPT, PPT, EPT, IX, IY,
1 WS, PRB, XN, YN, NDE, KFS, NNP, NEL )
C
CALL GMIN ( WP1, PS1, WPT, PPT, EPT, IX, IY, IS, JS, IE )
C
LOP = 0
TIM = 0.D0
C
999 CONTINUE
C
LOP = LOP + 1
TIM = TIM + DT
C
CALL PSCL ( PS1, WP1, PPT, IX, IY, PS1, WP1, ITT, NNP )
C
CALL WFVR ( WP1, PS1, WPT, IX, IY, IS, JS, IE, WM1,
1 PS2, PS, WP2, EPT, NNP )
C
C -------------------------
IF ( KFS.EQ.4) THEN
C
CALL WUC4 ( U, WP1, WM1, PS1, WPT, IX, IY )
C
CALL WVC4 ( V, WP1, WM1, PS1, WPT, EPT, IX, IY )
C
END IF
C -------------------------
IF ( KFS.EQ.2) THEN
C
CALL WUC2 ( U, WP1, WM1, PS1, WPT, IX, IY, NNP, NEL )
C
CALL WVC2 ( V, WP1, WM1, PS1, WPT, EPT, IX, IY, NNP, NEL )
C
END IF
C -------------------------
C
CALL FSUB ( U, V, WP1, WM1, WPT, EPT, IX, IY,
1 PS1, PS, WS, XN, YN, IE, JS, X1R,
2 TIM, ITT, NNP, *333 )
C
IF ( LOP.LT.LMT) GO TO 999
C
333 CONTINUE
C
CALL RPLT ( U, V, UU, VV, UE, VE, PS, PS2,
1 XN, YN, NDE, XYZ, XE, YE, PRB, IX,
2 IY, PPT, WPT, FLV, NNP, NEL )
C
CLOSE ( 10)
CLOSE ( 12)
CLOSE ( 15)
C
STOP
END
C **********************************************************************
C
SUBROUTINE BGIN ( PRB, IY, IX, KFS, NNP, NEL )
C
IMPLICIT REAL*8 (A-H,O-Z )
C
CHARACTER*14 PRB
C
COMMON /CNS/ NPB, WEI, BV, RV, KCR, OMW
C
WRITE(6,*) ' '
WRITE(6,2000)
2000 FORMAT(
1 ' ',' C **************************************************** C'/
2 ' ',' C C'/
3 ' ',' C F2D - GFD C'/
4 ' ',' C C'/
5 ' ',' C Fourth-ordered ( or Second-ordered ) C'/
5 ' ',' C C'/
6 ' ',' C Finite Difference Method C'/
6 ' ',' C C'/
7 ' ',' C ( 7 - Flow Problems ) C'/
8 ' ',' C C'/
9 ' ',' C ( V.5.3 ) C'/
A ' ',' C C'/
B ' ',' C Copyright 2011 : Yasuhiro MATSUDA C'/
C ' ',' C C'/
D ' ',' C **************************************************** C'/)
C
WRITE(6,2100)
2100 FORMAT(' ',' ** Which Problem ? ')
WRITE(6,*) ' '
WRITE(6,*) '1: V-Shape 2: Square, 3: Circle'
WRITE(6,*) '4: Backstep (1) 5: Backstep (2) 6: Cavity(S)'
WRITE(6,*) '7: Cavity( L)'
C
READ(5,*) NPB
C
IF ( NPB.EQ.1) PRB =' V-Shape '
IF ( NPB.EQ.2) PRB =' Square '
IF ( NPB.EQ.3) PRB =' Circle '
C
IF ( NPB.EQ.4) PRB =' Backstep (1)'
IF ( NPB.EQ.5) PRB =' Backstep (2)'
C
IF ( NPB.EQ.6) PRB =' Cavity (S)'
IF ( NPB.EQ.7) PRB =' Cavity (L)'
C
C ----- Attention ---
C
IX = 53
IY = 27
C
C ----- Back step ---
C
IF ( NPB.EQ.4 .OR. NPB.EQ.5) THEN
IX = 45
IY = 15
END IF
C
C ----- Cavity Flow with many meshes ---
C
IF ( NPB.EQ.7) THEN
IX = 51
IY = 51
END IF
C
C ----- Cavity Flow with lesser meshes ---
C
IF ( NPB.EQ.6) THEN
IX = 11
IY = 11
END IF
C
C ----- Attention ----- Max. Region : X * Y ---------------
C
C --- To Define New " NNP & NEL " for each problem ---
C
NNP = IX* IY
NEL = ( IX-1)*( IY-1)
C
C ---------------------------------------------------------
C
WRITE(6,*) ' ** Solution Order = ? 4 or 2'
** READ(5,*) KFS
KFS = 4
C =============
C
WRITE(6,2200) NPB, PRB, IY, IX, KFS
2200 FORMAT(/,' ** NPB =',I3,' * Problem :',A14,' ( IY=',I3,' IX =',
1 I3,' )',' * KFS =', I2/)
C
C ---------------------------------------------------
C
IF ( NPB.EQ.1) OPEN ( 10, FILE='PFDM1.DAT')
IF ( NPB.EQ.2) OPEN ( 10, FILE='PFDM2.DAT')
IF ( NPB.EQ.3) OPEN ( 10, FILE='PFDM3.DAT')
C
IF ( NPB.EQ.4) OPEN ( 10, FILE='PFDM7.DAT')
IF ( NPB.EQ.5) OPEN ( 10, FILE='PFDM4.DAT')
IF ( NPB.EQ.6) OPEN ( 10, FILE='PFDM6.DAT')
IF ( NPB.EQ.7) OPEN ( 10, FILE='PFDM5.DAT')
C
C ---------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE CFSB ( ALF, KY )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
COMMON /ABCD/ A1, A2, A3, A4, B1, B2, B3, B4, C1, C2, C3, C4,
1 D1, D2, D3, D4
C
DATA A11/ 0.58333 33333 33333D0/, ! 7/12
1 A12/ 0.625D0/, ! 5/8
2 A13/ -0.08333 33333 33333D0/, ! -1/12
3 A14/ -0.125D0/, ! -1/8
4 A21/ -0.08333 33333 33333D0/, ! -1/12
5 A22/ -0.04166 66666 66667D0/, ! -1/24
6 A23/ 0.08333 33333 33333D0/, ! 1/12
7 A24/ 0.04166 66666 66667D0/ ! 1/24
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
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
IF ( KY.EQ.0) THEN
B = ALF - 1.D0
A11M = A11
A13M = A13
END IF
C
IF ( KY.EQ.1) THEN
B = ALF + 1.D0
A11M = - A11
A13M = - A13
END IF
C
BB = B * B
BBB = B * BB
BBBB = B * BBB
C
C1 = A11 * B + A11M
C2 = A12 * BB - A12
C3 = A13 * BBB + A13M
C4 = A14 * BBBB - A14
C
D1 = A21 * B
D2 = A22 * BB
D3 = A23 * BBB
D4 = A24 * BBBB
C
RETURN
END
C **********************************************************************
C
SUBROUTINE DIST ( PS1, X, Y, IX, IY, IE, JS, X1R )
C
IMPLICIT REAL*8 (A-H,O-Z )
C
COMMON /CNS/ NPB, WEI, BV, RV, KCR, OMW
C
DIMENSION PS1( IX,IY), X( IX,IY), Y( IX,IY), DIS( 10)
C
NI = 0
C
IF ( NPB.EQ.5) IE = 1
C ----- Attention -----------
C
C --------------------------------------------------------------------
DO 100 I = IE, IX-1
C
IF (( PS1( I,2)+1.D0)* ( PS1( I+1,2)+1.D0).GT.1.D-7) GO TO 100
C
NI = NI + 1
DIS( NI) = X( I,2) + ( -1.D0 - PS1( I,2))
1 /( PS1( I+1,2) - PS1( I,2))* ( X( I+1,2) - X( I,2))
C
100 CONTINUE
C --------------------------------------------------------------------
C
II = NI
C
DO 200 I = 1, NI
C
DIS( I) = ( DIS( I) - X( IE,2))/( Y( 1,JS) - Y( 1,1))
IF ( DIS( I).GT.2.D0 .AND. DIS( I).LT.12.D0) II = I
C
200 CONTINUE
C --------------------------------------------------------------------
C
X1R = DIS( II)
C
RETURN
END
C **********************************************************************
C
SUBROUTINE FSUB ( U, V, WP1, WM1, WPT, EPT, X,
1 Y, PS1, PS, WS, XN, YN, IE,
2 JS, X1R, TIM, ITT, NNP, * )
C
IMPLICIT REAL*8 (A-H,O-Z )
C
COMMON /CNS/ NPB, WEI, BV, RV, KCR, OMW
COMMON /IDT/ DT, LMT, UL, VU, RD, VIS, REN, LOP
COMMON /UDT/ DX, DY, DX2, DY2, RS, RS2, CF, DTX, DTY
C
COMMON /HST/ WMX( 5000), PMA( 5000)
C
INTEGER X, Y
C
DIMENSION WP1( X,Y), WM1( X,Y), PS1( X,Y), PS( X,Y),
1 U ( X,Y), V ( X,Y), WS ( X,Y), XN( NNP),
2 YN ( NNP)
C
LOGICAL*1 WPT( X,Y,4), EPT( X,Y,4)
C
IF ( KCR.EQ.1) CALL WDCR ( WP1, WM1, WPT, EPT, X, Y, PS1, PS,
1 U, V )
C
IF ( KCR.EQ.1) GO TO 123
C
IF ( KCR.NE.1) CALL WDNC ( WP1, WM1, WPT, EPT, X, Y, PS1, PS )
C
123 CONTINUE
C
IF ( LOP.EQ.1) GO TO 999
C
C ----- Var. of Vorticity ---
C
WL1 = 0.D0
WDEL = 0.D0
C
C ----------------------------------------------------------
DO 100 J = 1,Y
DO 100 I = 1,X
C
WL1 = DMAX1( WL1, DABS( WP1( I,J)))
100 CONTINUE
C ----------------------------------------------------------
WL2 = 0.01D0* WL1
C
C ----------------------------------------------------------
DO 200 J = 1, Y
DO 200 I = 1, X
C
IF ( DABS( WP1( I,J)).LT.WL2) GO TO 200
WDEL = DMAX1( WDEL, DABS(( WP1( I,J) - WS( I,J))/ WP1( I,J)))
C
200 CONTINUE
C ----------------------------------------------------------
WMX( LOP) = 100.D0* WDEL
C
C ----- Attention --- PS1: All Zero ---
C
PMX = 0.D0
PDL = 0.D0
C
C ------------------------------------------------
DO 300 J = 1,Y
DO 300 I = 1,X
C
PMX = DMAX1( PMX, DABS( PS1( I,J)))
300 CONTINUE
C ------------------------------------------------
PL2 = 0.01D0* PMX
C
C ------------------------------------------------
DO 400 J = 1, Y
DO 400 I = 1, X
C
IF ( DABS( PS1( I,J)).LT.PL2) GO TO 400
PDL = DMAX1( PDL, DABS(( PS1( I,J) - PS( I,J))/ PS1( I,J)))
C
400 CONTINUE
C ------------------------------------------------
PMA( LOP) = 100.D0* PDL
C
C ----------------------------------------------------------
IF ( NPB.EQ.5.OR.NPB.EQ.4) THEN
C
C ----- To estimate an attachmennt length ---
C
VXL = 0.D0
SXR = X1R
C
C -----------------------------------------------------
C
CALL DIST ( PS1, XN, YN, X, Y, IE, JS, X1R )
C
C -----------------------------------------------------
IF ( X1R.LE.0.0001D0) GO TO 1122
C
VXL = DABS((X1R - SXR)/ X1R)* 100.D0
C
C ----- VXL > = 0.5 (%) ==> Stationary condition ---
C
C --------------------------------------------------
IF ( VXL.LT.0.5D0) THEN
WRITE(6,2000) LOP
2000 FORMAT(/,' *** Stationary condition reached : Loop =',I4,/)
WRITE(6,2100) LOP, TIM, ITT, WMX( LOP), PMA( LOP)
RETURN 1
END IF
C --------------------------------------------------
C
1122 CONTINUE
C
END IF
C ----------------------------------------------------------
IF ( PMA( LOP).LE.0.1) THEN ! Var_Str_Function < = 0.1 (%)
C
WRITE(6,2000) LOP
WRITE(6,2100) LOP, TIM, ITT, WMX( LOP), PMA( LOP)
RETURN 1
C
END IF
C ----------------------------------------------------------
IF ( NPB.EQ.5.OR.NPB.EQ.4) THEN
C
WRITE(6,2200) LOP, TIM, ITT, WMX( LOP), PMA( LOP), X1R, VXL
2200 FORMAT(' *Lp=',I4,' T=',F7.3,' IT=',I3,' V_W=',F7.2,' %',
1 ' V_P=',F7.2,' %',' Xr=',F6.2,' V_X=',F7.2)
GO TO 999
C
END IF
C ----------------------------------------------------------
C
WRITE(6,2100) LOP, TIM, ITT, WMX( LOP), PMA( LOP)
2100 FORMAT(' * Loop =',I4,' TIME=',F8.4,' ITT=', I4,' Var.W=',F8.3,
1 '(%)',' Var.P=',F8.3,'(%)' )
C
999 CONTINUE
C
C ----- To store present PS1 & WP1 -----
C
DO 500 J = 1, Y
DO 500 I = 1, X
C
PS( I,J) = PS1( I,J)
WS( I,J) = WP1( I,J)
C
500 CONTINUE
C --------------------------------------
C
IF ( LOP.EQ.LMT) WRITE(6,*) ' '
IF ( LOP.EQ.LMT) WRITE(6,*) ' * Normal end *'
C
RETURN
END
C **********************************************************************
C
SUBROUTINE GMIN ( WP1, PS1, WPT, PPT, EPT, IX, IY, IS, JS, IE )
C
IMPLICIT REAL*8 (A-H,O-Z )
C
LOGICAL*1 OMG, PSI, CNT, SLP, MVE, WPT, PPT, EPT
C
NAMELIST /SPCS/ OMG, PSI, SLP, MVE, CNT
C
COMMON /CNS/ NPB, WEI, BV, RV, KCR, OMW
C
DIMENSION WP1( IX,IY), PS1( IX,IY), WPT( IX,IY,4),
1 PPT( IX,IY,4), EPT( IX,IY,4), IDX( 5),
2 IDY( 5), IDE( 5)
C
C ----------------------------
C
IER = 0
ICT = 0
C
C ----------------------------
DO 100 J = 1, IY
DO 100 I = 1, IX
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 170 I = 1, IY
C
EPT( 1, I,4) = .TRUE.
EPT( IX,I,3) = .TRUE.
C
170 CONTINUE
C ----------------------------
C
C --------------------------------------
900 CONTINUE
C
ICRD = 1
C
OMG = .FALSE.
PSI = .FALSE.
SLP = .FALSE.
MVE = .FALSE.
CNT = .FALSE.
C
READ(10,SPCS,END=310,ERR=202)
*** WRITE(6,2100) OMG, PSI, SLP, MVE, CNT
C
ICT = ICT + 1
C
IF ( OMG ) GO TO 205
IF ( PSI ) GO TO 205
IF ( CNT) GO TO 205
C
202 IER = 1
WRITE(6,2000)
2000 FORMAT(' * Error in SPEC DATA')
WRITE(6,2100) OMG, PSI, SLP, MVE, CNT
2100 FORMAT(' *** OMG, PSI, SLP, MVE, CNT ', 5L2)
GO TO 315
C
205 CONTINUE
C
C ------------------------------------------------
DO 300 ICNT = 1, ICRD
C
ISLP = - 10
C
C --------------------------------------
DO 210 L = 1, 5
C
IDX( L) = 0
IDY( L) = 0
IDE( L) = 0
C
210 CONTINUE
C --------------------------------------
C
WVL = 0.D0
PVL = 0.D0
C
READ(10,1000,END=310,ERR=225) ISLP, WVL, PVL,
1 ( IDX( L), IDY( L), IDE( L), L =1, 5 )
1000 FORMAT( I2,1X,F6.3,1X,F6.3,5(3I3,1X ))
C
IF ( ICT.EQ.2) THEN
IS = IDX( 1)
JS = IDE( 1)
IE = IDX( 2)
END IF
C
IPT = 0
C
C ISLP = -1 : Boundary slopes downward
C ISLP = 0 : Boundary is horiz.
C ISLP = 1 : Boundary slopes upward
C ISLP = 10 : Boundary is vert.
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

225 IER = 1
WRITE(11,2200)
2200 FORMAT(' * Read Error or Error in SLOPE ')
WRITE(11,2300) ISLP, WVL, PVL,
1 ( IDX( L), IDY( L), IDE( L), L=1, 5 )
2300 FORMAT( 5X, I2, 1X, 2( F6.3,1X ), 5( 3I3,1X ))
GO TO 300
C
230 IPT = IPT + 1
231 IPT = IPT + 1
232 IPT = IPT + 1
233 IPT = IPT + 1
C
C IPT = 1 : Slopes downward. IPT = 2 : Horizontal
C IPT = 3 : Slopes upward. IPT = 4 : Vetical
C
C --------------------------------------
DO 400 L = 1, 5
C
IUPR = IDE( L)
GO TO ( 240, 240, 240, 235), IPT
C
235 ILOW = IDY( L)
GO TO 245
C
240 ILOW = IDX( L)
245 I = IDX( L)
J = IDY( L)
C
C ----- Initial X & Y coordinates ---
C
C ----------------------------
DO 450 INX = ILOW, IUPR
C
IF ( I) 255, 400, 250
C
250 IF ( J) 255, 400, 251
251 IF ( I-IX) 252, 252, 255
252 IF ( J-IY) 260, 260, 255
C
C ----- To check to see if this point is within aRVay boundaries
C
255 IER = 1
C
WRITE(11,2400) L, IX,IY
2400 FORMAT(' * Error in coordinate spec. number', I1,
1 '. IX =', I3, ' IY =', I3)
WRITE(11,2300) ISLP, WVL, PVL, ( IDX(M), IDY(M), IDE(M), M=1,5)
GO TO 400
C
260 CONTINUE
C
IF ( .NOT. OMG ) GO TO 265
C --------------------------------
C
WPT( I,J,1) = .FALSE.
WPT( I,J,2) = SLP
C
WPT( I,J,3) = MVE
WP1( I,J) = WVL
C
265 IF ( .NOT. PSI) GO TO 270
C
PPT( I,J,1) = .FALSE.
PPT( I,J,2) = SLP
PPT( I,J,3) = MVE
PS1( I,J) = PVL
C
270 IF ( .NOT. CNT ) GO TO 280
C
WPT( I,J,4) =.TRUE.
PPT( I,J,4) =.TRUE.
C
280 GO TO ( 286, 288, 287, 285 ), IPT
C
285 J = J + 1
GO TO 450
C
286 J = J - 1
GO TO 288
C
287 J = J + 1
288 I = I + 1
C
450 CONTINUE
C ----------------------------
C
400 CONTINUE
C --------------------------------------
C
300 CONTINUE
C ------------------------------------------------
C
GO TO 900
C
310 CONTINUE
WRITE(6,*) ' '
C
C ------------------------------------------------
C
CALL GPRN ( WPT, EPT, IX, IY )
C
C ------------------------------------------------
C
CALL IGSS ( PS1, PPT, WP1, IX, IY )
C
C ------------------------------------------------
C
RETURN
C
315 CONTINUE
C
WRITE(6,*) ' *** Error ***'
C
STOP
END
C **********************************************************************
C
SUBROUTINE GPRN ( ARY, ARY2, IX, IY )
C
IMPLICIT REAL*8 (A-H,O-Z )
C
COMMON /CNS/ NPB, WEI, BV, RV, KCR, OMW
C
LOGICAL*1 ARY( IX,IY,4), ARY2( IX,IY,4)
C
CHARACTER LBF( 53), LTC, LTS, LTM, LTP, LTR, LOO
C
C ***** Attention : SIZE = 53 ***
C
DATA LTC/'C'/, LTS/'S'/, LTM/'M'/, LTP/'P'/, LTR/'R'/, LOO/' '/
C
C ----------------------------
DO 100 I = 1, 53
C
LBF( I) = LOO
100 CONTINUE
C ----------------------------
C
N = 1
NE = 53
IF ( IX.LT.53) NE = IX
C
99 CONTINUE
C
C ----------------------------------------------------------
DO 150 JJ = 1, IY
C
J = IY+1-JJ
C
C ------------------------------------------------
DO 200 I = N, NE
C
IF ( ARY( I,J,1)) GO TO 125
IF ( ARY( I,J,2)) GO TO 120
IF ( ARY( I,J,3)) GO TO 115
C
LBF( I+1-N) = LTR
C
C ------ Rigid, No SLP boundary ---
C
GO TO 200
C
115 LBF( I+1-N) = LTM
C
C ----- Moving, No SLP boundary ---
C
GO TO 200
C
120 LBF( I+1-N) = LTS
C
C ----- Interior to a boundary or a SLP boundary ---
C
GO TO 200
C
125 CONTINUE
C
DO 130 K = 1, 4
C
IF ( ARY2( I,J,K)) GO TO 135
130 CONTINUE
C
LBF( I+1-N) = LOO
GO TO 200
C
135 CONTINUE
C
IF ( ARY( I,J,4)) GO TO 140
C
LBF( I+1-N) = LTP
C
C ----- Edge point with periodic boundary conditions ---
C
GO TO 200
C
140 LBF( I+1-N) = LTC
C
C ----- Edge point with contnuative boundary conditions ---
C
200 CONTINUE
C ------------------------------------------------
C
WRITE(6,2000) ( LBF( I), I = 1, 53)
2000 FORMAT(' ',8X,53A1)
C
C ----- Print 1 line of alphabetic characters ---
C
150 CONTINUE
C ----------------------------------------------------------
C
WRITE(6,*) ' '
IF ( N + 52. GE.IX) RETURN
IF ( N + 104.LE.IX) GO TO 300
NCVL = N + 52 - IX
N = IX - 52
NE = N + 52
C
WRITE(6,2100) NCVL
2100 FORMAT(/,' Overlay point array with ',I2,' columns')
GO TO 99
C
300 N = N + 52
NE = N + 52
GO TO 99
C
END
C **********************************************************************
C
SUBROUTINE IGSS ( PS1, PPT, WP1, IX, IY )
C
IMPLICIT REAL*8 (A-H,O-Z )
C
DIMENSION PS1( IX,IY), WP1( IX,IY)
C
LOGICAL*1 PPT( IX,IY,4)
C
COMMON /CNS/ NPB, WEI, BV, RV, KCR, OMW
COMMON /IDT/ DT, LMT, UL, VU, RD, VIS, REN, LOP
COMMON /UDT/ DX, DY, DX2, DY2, RS, RS2, CF, DTX, DTY
C
N = IY - 1
C
FLOW = ( PS1( IX,IY) - PS1( IX,1))/( IY-1)
C
C --------------------------------------
DO 100 J = 2, N
C
PS1( IX,J) = PS1( IX,J-1) + FLOW
100 CONTINUE
C --------------------------------------
DO 200 J = 1, IY
DO 200 I = 1, IX
C
IF (.NOT.PPT( I,J,1)) GO TO 200
PS1( I,J) = PS1( IX,J) ! Fluid point
C
200 CONTINUE
C --------------------------------------
C
IF ( .NOT.PPT( 1,1,3)) UL = 0.D0 ! Not moving
C
RETURN
END
C **********************************************************************
C
SUBROUTINE INIT ( WP1, PS1, S, PS, WPT, PPT, EPT, IX, IY, WS )
C
IMPLICIT REAL*8 (A-H,O-Z )
C
LOGICAL*1 WPT( IX,IY,4), PPT( IX,IY,4), EPT( IX,IY,4)
C
DIMENSION PS1( IX,IY), WP1( IX,IY), PS( IX,IY), S( IX,IY),
1 WS ( IX,IY)
C
COMMON /CNS/ NPB, WEI, BV, RV, KCR, OMW
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
PS1( I,J) = 0.D0
S ( I,J) = 0.D0
PS ( I,J) = 0.D0
WS ( I,J) = 0.D0
C
100 CONTINUE
C ------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE INPT ( WP1, PS1, S, PS, WPT, PPT, EPT, IX, IY,
1 WS, PRB, XN, YN, NDE, KFS, NNP, NEL )
C
IMPLICIT REAL*8 (A-H,O-Z )
C
LOGICAL*1 WPT( IX,IY,4), PPT( IX,IY,4), EPT( IX,IY,4)
CHARACTER*14 PRB
C
COMMON /CNS/ NPB, WEI, BV, RV, KCR, OMW
COMMON /IDT/ DT, LMT, UL, VU, RD, VIS, REN, LOP
COMMON /UDT/ DX, DY, DX2, DY2, RS, RS2, CF, DTX, DTY
C
DIMENSION WP1( IX,IY), PS1( IX,IY), S ( IX,IY), PS ( IX,IY),
1 WS ( IX,IY), XN ( NNP), YN( NNP), NDE( NEL,4)
C
NAMELIST /IDAT/ RS, UL, VU, D
C
C ***** Attention : To Define New " NNP & NEL " for each problem ***
C
CALL BGIN ( PRB, IY, IX, KFS, NNP, NEL )
C
C -------------------------------------------------------------------
C
CALL INIT ( WP1, PS1, S, PS, WPT, PPT, EPT, IX, IY, WS )
C
C -------------------------------------------------------------------
C
DX = 2.D0/ DFLOAT( IX-1)
DY = 1.D0/ DFLOAT( IY-1)
C
C -------------------------------
C
IF ( NPB.EQ.7.OR.NPB.EQ.6) DX = 1.D0/ DFLOAT( IX-1)
C
C ----- Default value ---
C
RD = 1.D0
C
READ(10,IDAT)
C ==================
C
C ----- DX : A* RS, DY : A ---
C
WRITE(6,*) ' * Reynolds Number = ?'
** READ(5,*) REN
REN = 100.D0
C =================
C
WRITE(6,*) ' * Weighting Factor = ? ( 0 or 0.5 or 1)'
** READ(5,*) WEI
C
WEI = 0.5D0
C ================
OMW = 1.D0 - WEI
C =====================
C
WRITE(6,*) ' * Correction = ? ( 1: YES, 2: NO )'
** READ(5,*) KCR
C
KCR = 1 ! Correction Method for Viscousity
C ============
C
DX = 2.D0/ DFLOAT( IX-1)
DY = 1.D0/ DFLOAT( IY-1)
C
C ----- Cavity problem ---
C
IF ( NPB.EQ.7 .OR. NPB.EQ.6 ) THEN
DX = 1.D0/ DFLOAT( IX-1)
DY = DX
RS = DX/ DY
END IF
C
C ----- Back step Problem ------------------------
C
IF ( NPB.EQ.5.OR.NPB.EQ.4) THEN
C
WRITE(6,2000)
2000 FORMAT(' ',' * RL = ? ')
*** READ(5,*) RLL
C
RLL = 30.D0
C ================
C
C ----- Attention --- Area : 2 X RLL ---
C
DX = RLL / DFLOAT( IX-1)
DY = 2.D0/ DFLOAT( IY-1)
C --------------------------------------
C
RS = DX/ DY
C
END IF
C --------------------------------------------------------
C
CALL STND ( XN, YN, NDE, IX, IY, NNP, NEL )
C
C --------------------------------------------------------
C
DX2 = DX* DX
DY2 = DY* DY
C
RS2 = RS* RS
CF = 2.D0* ( RS2 + 1.D0 )
C
VIS = UL* RD/ REN
C
IF ( NPB.EQ.5.OR.NPB.EQ.4) VIS = 1.D0/ REN
IF ( NPB.EQ.7.OR.NPB.EQ.6) VIS = 1.D0/ REN
C
WRITE(6,*) ' * Number of Loop = ?'
** READ (5,*) LMT
C
LMT = 200
C
IF ( NPB.LE.3) LMT = 1000
IF ( NPB.EQ.7) LMT = 1000
C ------------------------------
C
WRITE(6,2100) LMT
2100 FORMAT(' ** Loop=',I4)
C
WRITE(6,*) ' * DT = ?'
** READ(5,*) DT
C
** IF ( NPB.EQ.7 ) READ(5,*) DT
C
C --------------------------------------
C
IF ( NPB.LE.3) DT = 0.005D0
IF ( NPB.EQ.4) DT = 0.02D0
C
IF ( NPB.EQ.5) DT = 0.02D0
IF ( NPB.EQ.6) DT = 0.1D0 ! S.Cav.
IF ( NPB.EQ.7) DT = 0.005D0 ! L.Cav.
C
C --------------------------------------
C
BV = VU* DT/ DY
C
C ----- DX : A* RS ---
C
IF ( NPB.LE.5) BV = UL* DT/ DX
C
RV = VIS* DT/( DX* DX)
C
DTX = DT/ DX
DTY = DT/ DY
C
XMX = 0.D0
YMX = 0.D0
XMN = 1000.D0
YMN = 1000.D0
C
C ------------------------------------------------
DO 100 I = 1, NNP
C
IF ( XN( I).GE.XMX) XMX = XN( I)
IF ( YN( I).GE.YMX) YMX = YN( I)
C
IF ( XN( I).LE.XMN) XMN = XN( I)
IF ( YN( I).LE.YMN) YMN = YN( I)
C
100 CONTINUE
C ------------------------------------------------
WRITE(6,*) ' '
WRITE(6,2200) XMX, XMN, YMX, YMN
2200 FORMAT(' * Xmax.=', F6.2,' Xmin.=', F6.2,' Ymax.=', F6.2,
1 ' Ymin.=',F6.2)
C
WRITE(6,2300) DX, DY, REN, DT, LMT
2300 FORMAT(/,' * DX =',F7.4,' DY =',F7.4,' Re =',F7.2,' DT =',
1 F7.4,' Loop=',I4)
C ------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE PSCL ( PS1, WP1, PPT, IX, IY, PS1L, WP1L, ITT, NNP )
C
IMPLICIT REAL*8 (A-H,O-Z )
C
DIMENSION PS1( IX,IY), WP1( IX,IY), PS1L( NNP), WP1L( NNP)
LOGICAL*1 PPT( IX,IY,4), DKY
C
COMMON /CNS/ NPB, WEI, BV, RV, KCR, OMW
COMMON /IDT/ DT, LMT, UL, VU, RD, VIS, REN, LOP
COMMON /UDT/ DX, DY, DX2, DY2, RS, RS2, CF, DTX, DTY
C
EPS = 0.00005D0
C
ACL = 1.5D0
C
IUU = 0
IUD = 0
ITT = 0
C
DKY = .FALSE.
C===================
C
JDMY = IY + 1
C
C ----------------------------------------------------------
900 IDFC = 0
ITT = ITT + 1
C
IF ( ITT.GT.1000 ) GO TO 98
C
C -------------------------------------------------
DO 100 JJ = 1, IY
C
10 CONTINUE
DKY = .NOT.DKY
C
J = JJ
IF ( DKY) J = JJ
IF (.NOT.DKY) J = JDMY - JJ
C
C --------------------------------------
DO 150 I = 1, IX
C
IF ( PPT( I,J,1)) GO TO 91 ! Fluid point
C
C ----- Except Fluid points ---
C
GO TO 150
C
91 CONTINUE
C
IF ( I.EQ.IX) GO TO 2
IF ( I.EQ.1 ) GO TO 12
C
C ----- Normal interior to fluid point ---
C
IP = ( J-1)*IX + I
C
DELP = ( PS1L( IP+1) + PS1L( IP-1)
1 + RS2* ( PS1L( IP+IX) + PS1L( IP-IX) + DY2* WP1L( IP)))
2 / CF
C
DDIF = DABS( DELP - PS1L( IP))
C
PS1L( IP) = PS1L( IP) + ACL* ( DELP - PS1L( IP))
IF ( DDIF.GT.EPS) GO TO 5
C
GO TO 150
C
5 IDFC = 1
IXX = I
IYY = J
GO TO 150
C
2 CONTINUE
C
IF ( PPT( I,J,4)) GO TO 8 ! Continuative
C
C ----- Periodic Right ---
C
PS1( I,J) = PS1( I-( IX-1),J)
GO TO 150
C
C ----- Continuative Right ---
C
8 PS1( I,J) = ( PS1( I,J+1) + PS1( I,J-1) + DY2* WP1( I,J))* 0.5D0
GO TO 150
C
12 IF ( PPT( I,J,4)) GO TO 18 ! Continuative
C
C ------ Perodic Left ---
C
DELP = ( PS1( I+1,J) + PS1( I-1+( IX-1),J) + RS2* ( PS1( I,J+1)
1 + PS1( I,J-1) + DY2* WP1( I,J)))/ CF
C
PS1( I,J) = PS1( I,J) + ACL* ( DELP - PS1( I,J))
GO TO 150
C
C ----- Continuative Left ---
C
18 PS1( I,J) = ( PS1( I,J) + PS1( I,J) + RS2* ( PS1( I,J+1)
1 + PS1( I,J-1) + DY2* WP1( I,J)))/ CF
C
150 CONTINUE
C --------------------------------------
C
IF ( DKY) GO TO 10
C
100 CONTINUE
C ------------------------------------------------
C
IF ( IUD.EQ.0) GO TO 101
C
PSS = 0.D0
C
C ------------------------------------------------
DO 300 I = 1, IX
C
PSS = PSS + PS1( I,2) + 0.5D0* DY2* WP1( I,1)
300 CONTINUE
C ------------------------------------------------
C
PSS = ( PSS/IX ) - DY* UL
C
C ------------------------------------------------
DO 400 I = 1, IX
C
PS1( I,1) = PSS
400 CONTINUE
C ------------------------------------------------
C
101 IF ( IUU.EQ.0) GO TO 102
C
PSS = 0.D0
C ------------------------------------------------
DO 500 I = 1, IX
C
PSS = PSS + PS1( I,IY-1) + 0.5D0* DY2* WP1( I,IY)
500 CONTINUE
C ------------------------------------------------
C
PSS = ( PSS/ IX) + DY* VU
C
C ------------------------------------------------
DO 600 I = 1, IX
C
PS1( I,IY) = PSS
600 CONTINUE
C ------------------------------------------------
C
IYY = ( IY/ 2) + 1
TMP = PS1( 1,IYY)
C
C ------------------------------------------------
DO 700 J = 1, IY
DO 700 I = 1, IX
C
PS1( I,J) = PS1( I,J) - TMP
700 CONTINUE
C ------------------------------------------------
C
102 IF ( IDFC.GT.0 ) GO TO 900
C ----------------------------------------------------------
RETURN
C
98 CONTINUE
WRITE(6,2000) IXX, IYY
2000 FORMAT(' ** ITT : Over 1000, ==> Check ', I4, 1X, I4,
1 ' * PSI - Solution incompl1ete : STOP ')
C
STOP
END
C **********************************************************************
C
SUBROUTINE RPLT ( U, V, UU, VV, UE, VE, PS,
1 PS2, XX, YY, NDE, XYZ, XE, YE,
2 PRB, IX, IY, PPT, WPT, FLV, NNP,
3 NEL )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
CHARACTER*14 PRB
C
COMMON /HST / WMX( 5000), PMA( 5000)
COMMON /CNS / NPB, WEI, BV, RV, KCR, OMW
COMMON /IDT / DT, LMT, UL, VU, RD, VIS, REN, LOP
COMMON /UDT / DX, DY, DX2, DY2, RS, RS2, CF, DTX, DTY
C
DIMENSION U ( IX,IY), V ( IX,IY), UU( NNP), VV( NNP),
1 PS ( IX,IY), PS2( NNP), XX( NNP), YY( NNP),
2 UE ( NEL), VE ( NEL), XE( NEL), YE( NEL),
3 XYZ( 2,NNP), NDE( NEL,4)
C
LOGICAL*1 PPT( IX,IY,4), WPT( IX,IY,4), FLV( NEL)
C
C ----------------------------------------------------------
C
CALL STMP ( U, UU, IX, IY, NNP )
CALL STMP ( V, VV, IX, IY, NNP )
CALL STMP ( PS, PS2, IX, IY, NNP )
C
CALL STNN ( XX, YY, XE, YE, NDE, IX, IY, NNP, NEL )
C
C ----------------------------------------------------------
DO 100 I = 1, NEL
C
UE( I) = ( UU( NDE( I,1)) + UU( NDE( I,2)))/ 2.D0
VE( I) = ( VV( NDE( I,1)) + VV( NDE( I,2)))/ 2.D0
C
100 CONTINUE
C ----------------------------------------------------------
C
IF ( NPB.EQ.7 .OR. NPB.EQ.6) THEN
C
C ----------------------------------------------------------
DO 200 I = 1, NEL
C
IF ( DABS( UE( I)).GT.1.1D0) THEN
WRITE(6,*) ' *** U > = 1.1 *** STOP *** U =', UE( I),' * NEL=', I
STOP
END IF
C
200 CONTINUE
C ----------------------------------------------------------
C
END IF
C
K = 0
C
C ----- Attention --------------------------------
C
DO 300 J = 1, IY
DO 300 I = 1, IX
C
K = K + 1
II = IX* ( IY-J) + I
C
XYZ( 1,K) = XX( II)
XYZ( 2,K) = YY( II)
C
300 CONTINUE
C ------------------------------------------------
C
UEMX = 0.D0
C
C ----------------------------------------------------------
DO 400 I = 1, NEL
C
IF ( DABS( UE( I)).GE.UEMX) UEMX = DABS( UE( I))
400 CONTINUE
C ----------------------------------------------------------
C
BV = UEMX* DT/ DX
C ======================
C
RV = VIS* DT/( DX* DX)
C
IF ( NPB.LE.3) REN = UEMX* RD/ VIS
IF ( NPB.EQ.5. OR. NPB.EQ.4) REN = UEMX* RD/ VIS
IF ( NPB.EQ.7 .OR. NPB.EQ.6) REN = VU * RD/ VIS
C =========================================================
C
IF ( UEMX.GE.1000.D0) WRITE(6,*) ' *** Umax. > = 1000, STOP'
IF ( UEMX.GE.1000.D0) STOP
C
XEMX = 0.D0
YEMX = 0.D0
C
C ----------------------------------------------------------
DO 500 I = 1, NEL
C
IF ( XE( I).GE.XEMX) XEMX = XE( I)
IF ( YE( I).GE.YEMX) YEMX = YE( I)
C
500 CONTINUE
C ----------------------------------------------------------
C
WRITE(6,2000) UEMX, BV, RV, REN
2000 FORMAT(/,' * Umax.=',F7.3,' * BV =',F8.4,' * RV =',F8.4,' * Re =',
1 F7.2)
C
WRITE(15,2100) NPB, PRB, WEI, REN, DT, LOP, NEL, NNP, IX, IY,
1 DX, DY
2100 FORMAT(I2,A14,3F9.3,5I4,2F9.3)
C
WRITE(15,2200) ( XE( I), YE( I), I = 1, NEL)
WRITE(15,2200) ( UE( I), VE( I), I = 1, NEL)
2200 FORMAT(10F8.4)
C
WRITE(15,2300) ( WMX( I), I = 1, LOP)
WRITE(15,2300) ( PMA( I), I = 1, LOP)
2300 FORMAT(10F8.2)
C
WRITE(15,2200) (( XYZ( I,J), I = 1, 2), J = 1, NNP)
WRITE(15,2200) ( PS2( I), I = 1, NNP)
WRITE(15,2400) (( NDE( I,J), J = 1, 4), I = 1, NEL)
2400 FORMAT(20I4)
C
C ----------------------------------------------------------
DO 600 I = 1, NEL
C
FLV( I) = .FALSE. ! Not Fluid points ( Not to Plot )
600 CONTINUE
C ----------------------------------------------------------
DO 700 J = 1, IY-1
DO 700 I = 1, IX-1
C
NEE = ( J-1)* ( IX-1) + I
IF ( WPT( I,J,2)) GO TO 700
FLV( NEE) = .TRUE. ! Fluid points ( To Plot )
C
700 CONTINUE
C ----------------------------------------------------------
C
WRITE(15,2500) ( FLV( I), I=1, NEL)
2500 FORMAT(80L1)
C
RETURN
END
C **********************************************************************
C
SUBROUTINE STAB ( PS1, EPT, IX, IY )
C
C ----- To check "DT" ---
C
IMPLICIT REAL*8 (A-H,O-Z )
C
DIMENSION PS1( IX,IY)
C
LOGICAL*1 EPT( IX,IY,4)
C
COMMON /CNS/ NPB, WEI, BV, RV, KCR, OMW
COMMON /IDT/ DT, LMT, UL, VU, RD, VIS, REN, LOP
COMMON /UDT/ DX, DY, DX2, DY2, RS, RS2, CF, DTX, DTY
C
SMAX = 0.D0
C
C ----------------------------------------------------------
DO 100 J = 1, IY
DO 100 I = 1, IX
C
IF ( EPT( I,J,2)) GO TO 100 ! Lower edge
UST = DABS ( PS1( I,J) - PS1( I,J-1)) ! No-SLP
IF ( UST.GT.SMAX) SMAX = UST
C
100 CONTINUE
C ----------------------------------------------------------
DO 200 J = 1, IY
DO 200 I = 1, IX
C
IF ( EPT( I,J,4)) GO TO 200 ! Left edge
UST = DABS ( PS1( I,J) - PS1( I-1,J)) ! Periodic
IF ( UST.GT.SMAX) SMAX = UST
C
200 CONTINUE
C ----------------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE STND ( X, Y, NDE, IX, IY, NNP, NEL )
C
IMPLICIT REAL*8 (A-H,O-Z )
C
COMMON /CNS/ NPB, WEI, BV, RV, KCR, OMW
COMMON /UDT/ DX, DY, DX2, DY2, RS, RS2, CF, DTX, DTY
C
DIMENSION X( NNP), Y( NNP), NDE( NEL,4)
C
NN = 0
C
C --------------------------------------
DO 100 J = 1, IY
DO 100 I = 1, IX
C
NN = NN + 1
X( NN) = DX* DFLOAT( I-1)
Y( NN) = DY* 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 --------------------------------------
C
RETURN
END
C ***********************************************************************
C
SUBROUTINE STMP ( DM1, DM2, IX, IY, NNP )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
COMMON /CNS/ NPB, WEI, BV, RV, KCR, OMW
C
DIMENSION DM1( IX,IY), DM2( NNP)
C
NN = 0
C
C -----Attention -----------------------
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 STNN ( XX, YY, XE, YE, NDE, IX, IY, NNP, NEL )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
COMMON /CNS/ NPB, WEI, BV, RV, KCR, OMW
COMMON /UDT/ DX, DY, DX2, DY2, RS, RS2, CF, DTX, DTY
C
DIMENSION XX( NNP), YY( NNP), NDE( NEL,4), XE( NEL), YE( NEL)
C
NN = 0
C
C --------------------------------------
DO 100 J = 1, IY
DO 100 I = 1, IX
C
NN = NN + 1
XX( NN) = DX* DFLOAT( I-1)
YY( NN) = DY* 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
C
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
XE( I) = ( XX( NDE( I,1)) + XX( NDE( I,2))
1 + XX( NDE( I,3)) + XX( NDE( I,4)))/ 4.D0
C
YE( 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 VORB ( WP1, PS1, WPT, IX, IY, NNP, NEL )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION WP1( IX,IY), PS1( IX,IY)
LOGICAL*1 WPT( IX,IY,4)
C
COMMON /CNS/ NPB, WEI, BV, RV, KCR, OMW
COMMON /IDT/ DT, LMT, UL, VU, RD, VIS, REN, LOP
COMMON /UDT/ DX, DY, DX2, DY2, RS, RS2, CF, DTX, DTY
C
C ----- Vorticity Boundary Conditions ---
C
C *** Attention : Moving condition ==> Check
C
C No-SLP surfaces updated ( Ref.: BS-S/W )
C
C ------------------------------------------------
DO 100 I = 1, IX
DO 100 J = 1, IY
C ------------------------------------------------
C
IF ( WPT( I,J,1)) GO TO 100 ! Fluid point
IF ( WPT( I,J,2)) GO TO 50 ! SLP
IF ( WPT( I,J,3)) GO TO 200 ! Moving
C
C ----- WPT( I,J,4) : Continuative ---
C
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 SLP surfaces ---
C
IF ( .NOT. WPT( I+1,J,1)) GO TO 1
WA = ( PS1( I,J) - PS1( I+1,J))/ DX2 ! Fluid point

13 IF ( .NOT. WPT( I-1,J,1)) GO TO 3
WB = ( PS1( I,J) - PS1( I-1,J))/ DX2

12 IF ( .NOT. WPT( I,J+1,1)) GO TO 2
WC = ( PS1( I,J) - PS1( I,J+1))/ DY2

14 IF ( .NOT. WPT( I,J-1,1)) GO TO 4
WD = ( PS1( I,J) - PS1( I,J-1))/ DY2
C
110 WP1( I,J) = WA + WB + WC + WD
GO TO 100
C
1 WA = ( PS1( I,J) - PS1( I-1,J))/ DX2 ! Boundary point
GO TO 13
C
3 WB = ( PS1( I,J) - PS1( I+1,J))/ DX2
GO TO 12
C
2 WC = ( PS1( I,J) - PS1( I,J-1))/ DY2
GO TO 14
C
4 WD = ( PS1( I,J) - PS1( I,J+1))/ DY2
GO TO 110
C
C ----- No SLP surfaces at edges ---
C
50 WP1( I,J) = 0.D0
GO TO 100
C
C ----- Moving boundary with Wood's Approx. ---
C
10 CONTINUE
C
WP1( I,J) = 3.D0* ( PS1( I,J) - PS1( I,J-1) - DY)
C *** Moving direction ***
1 / DY2 - 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 WP1( I,J) = 3.D0* ( PS1( I,J) - PS1( I,J+1))
1 / DY2 - 0.5D0* WP1( I,J+1)
GO TO 100
C
30 WP1( I,J) = 3.D0* ( PS1( I,J) - PS1( I-1,J))
1 / DX2 - 0.5D0* WP1( I-1,J)
GO TO 100
C
40 WP1( I,J) = 3.D0* ( PS1( I,J) - PS1( I+1,J))
1 / DX2 - 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 WDCR ( WP1, WM1, WPT, EPT, IX, IY, PS1, PS, U, V )
C
IMPLICIT REAL*8 (A-H,O-Z )
C
DIMENSION WP1( IX,IY), WM1( IX,IY), U( IX,IY), V( IX,IY),
1 PS1( IX,IY), PS ( IX,IY)
C
LOGICAL*1 WPT( IX,IY,4), EPT( IX,IY,4)
C
COMMON /CNS/ NPB, WEI, BV, RV, KCR, OMW
COMMON /IDT/ DT, LMT, UL, VU, RD, VIS, REN, LOP
COMMON /UDT/ DX, DY, DX2, DY2, RS, RS2, CF, DTX, DTY
C
C ------------------------------------------------
DO 200 J = 1, IY
DO 200 I = 1, IX
C
IF ( WPT( I,J,1)) GO TO 200
WP1( I,J) = WM1( I,J)
C
200 CONTINUE
C ------------------------------------------------
DO 100 J = 1, IY
DO 100 I = 1, IX
C ------------------------------------------------
C
C ----- Artificial Viscousity --- ( RS = DX/ DY ) -----------
C
AVSC = VIS + U( I,J)* V( I,J)* DT* RS/( RS* RS + 1.D0)
C -----------------------------------------------------------
C
C ----- Art. Viscousity Term ---
C
AVT = AVSC* DT/ DY2
C ------------------------------
C
IF ( WPT( I,J,1)) GO TO 1
IF ( WPT( I,J,2)) GO TO 150
IF ( WPT( I,J,3)) GO TO 100
C
C ----- Rigid, No-SLP, Not moving ---
C
GO TO 150
C
1 IF ( EPT( I,J,3)) GO TO 2
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
C ----- Fluid point ---
C
C ----- To compute Diffusion term --- dW/dT = Art_V* Nbr-W ------------
C
WP1( I,J) = WP1( I,J) + AVT* ( WM1( I,J+1) + WM1( I,J-1)
1 + ( WM1( I+1,J) + WM1( I-1,J) - CF* WM1( I,J))/ RS2 )
C
C ----------------------------------------------------------------------
C
GO TO 100
C
2 CONTINUE
C
IF ( WPT( I,J,4)) GO TO 8
C
C ------ Periodic Right ---
C
WP1( I,J) = WP1( I,J) + AVT* ( WM1( I,J+1) + WM1( I,J-1)
1 + ( WM1( I+1-( IX-1),J) + WM1( I-1,J) - CF* WM1( I,J))
2 / RS2 )
GO TO 100
C
C ----- Continuative Right ---
C
8 IF ( NPB.NE.7.OR.NPB.NE.6) WP1( I,J) = WP1( I-1,J)
C
IF ( NPB.EQ.7.OR.NPB.EQ.6)
1 WP1( I,J) = WP1( I,J)
2 + AVT* ( WM1( I,J+1) + WM1( I,J-1) - 2.D0* WM1( I,J))
GO TO 100
C
12 IF ( WPT( I,J,4)) GO TO 18
C
C ----- Periodic Left ---
C
WP1( I,J) = WP1( I,J) + AVT* ( WM1( I,J+1) + WM1( I,J-1)
1 + ( WM1( I+1,J) + WM1( I-1+( IX-1),J) - CF* WM1( I,J))/ RS2 )
GO TO 100
C
C ----- Continuative Left ---
C
18 WP1( I,J) = WP1( I,J) + AVT* ( WM1( I,J+1) + WM1( I,J-1)
1 + ( WM1( I,J) + WM1( I,J) - CF* WM1( I,J))/ RS2 )
GO TO 100
C
C ----- Free SLP - No diffusion ---
C
150 GO TO 100
C
100 CONTINUE
C ------------------------------------------------
C
C --------------------------------------
DO 400 J = 1, IY
DO 400 I = 1, IX
C
WM1( I,J) = PS( I,J)
400 CONTINUE
C --------------------------------------
RETURN
C
900 CONTINUE
C
RETURN
END
C **********************************************************************
C
SUBROUTINE WDNC ( WP1, WM1, WPT, EPT, IX, IY, PS1, PS )
C
IMPLICIT REAL*8 (A-H,O-Z )
C
DIMENSION WP1( IX,IY), WM1( IX,IY), PS1( IX,IY), PS ( IX,IY)
C
LOGICAL*1 WPT( IX,IY,4), EPT( IX,IY,4)
C
COMMON /CNS/ NPB, WEI, BV, RV, KCR, OMW
COMMON /IDT/ DT, LMT, UL, VU, RD, VIS, REN, LOP
COMMON /UDT/ DX, DY, DX2, DY2, RS, RS2, CF, DTX, DTY
C
C ------------------------------------------------
DO 100 J = 1, IY
DO 100 I = 1, IX
C
IF ( WPT( I,J,1)) GO TO 100
WP1( I,J) = WM1( I,J)
C
100 CONTINUE
C ------------------------------------------------
C
C ------------------------------------------------
DO 200 J = 1, IY
DO 200 I = 1, IX
C ------------------------------------------------
C
C ----- In case of No-Correction for Viscousity ---
C
AVT = VIS* DT/ DY2
C -------------------------------------------------
IF ( WPT( I,J,1)) GO TO 1
IF ( WPT( I,J,2)) GO TO 200
IF ( WPT( I,J,3)) GO TO 200
C
C ----- Rigid, No-SLP, Not moving ---
C
GO TO 200
C
1 IF ( EPT( I,J,3)) GO TO 2
IF ( EPT( I,J,4)) GO TO 12
IF ( EPT( I,J,1)) GO TO 300
IF ( EPT( I,J,2)) GO TO 300
C
C ----- Fluid point ---
C
WP1( I,J) = WP1( I,J) + AVT* ( WM1( I,J+1) + WM1( I,J-1)
1 + ( WM1( I+1,J) + WM1( I-1,J) - CF* WM1( I,J))/ RS2 )
GO TO 200
C
2 IF ( WPT( I,J,4)) GO TO 8
C
C ----- Periodic Right ---
C
WP1( I,J) = WP1( I,J) + AVT* ( WM1( I,J+1) + WM1( I,J-1)
1 + ( WM1( I+1-( IX-1),J) + WM1( I-1,J) - CF* WM1( I,J)))/ RS2
GO TO 200
C
C ----- Continuative Right ---
C
8 IF ( NPB.NE.7.OR.NPB.NE.6) WP1( I,J) = WP1( I-1,J)
C
IF ( NPB.EQ.7.OR.NPB.EQ.6)
1 WP1( I,J) = WP1( I,J) + AVT* ( WM1( I,J+1) + WM1( I,J-1)
2 - 2.D0* WM1( 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) + AVT* ( WM1( I,J+1) + WM1( I,J-1)
1 + ( WM1( I+1,J) + WM1( I-1 + ( IX-1),J) - CF* WM1( I,J))
2 / RS2 )
GO TO 200
C
C ----- Continuative Left ---
C
18 WP1( I,J) = WP1( I,J) + AVT* ( WM1( I,J+1) + WM1( I,J-1)
1 + ( WM1( I,J) + WM1( I,J) - CF* WM1( I,J))/ RS2 )
GO TO 200
C
C ----- Free-SLP, No-diffusion ---
C
200 CONTINUE
C ------------------------------------------------
DO 220 J = 1, IY
DO 220 I = 1, IX
C
WM1( I,J) = PS ( I,J)
PS ( I,J) = PS1( I,J)
C
220 CONTINUE
C ------------------------------------------------
RETURN
C
300 CONTINUE
C
RETURN
END
C **********************************************************************
C
SUBROUTINE WFVR ( WP1, PS1, WPT, IX, IY, IS, JS, IE, WM1,
1 PS2, PS, WP2, EPT, NNP )
C
IMPLICIT REAL*8 (A-H,O-Z )
C
DIMENSION PS ( IX,IY), PS1( IX,IY), PS2( NNP),
1 WP1( IX,IY), WP2( NNP), WM1( IX,IY)
C
LOGICAL*1 WPT( IX,IY,4), EPT( IX,IY,4)
C
COMMON /HST/ WMX( 5000), PMA( 5000)
COMMON /CNS/ NPB, WEI, BV, RV, KCR, OMW
COMMON /IDT/ DT, LMT, UL, VU, RD, VIS, REN, LOP
COMMON /UDT/ DX, DY, DX2, DY2, RS, RS2, CF, DTX, DTY
C
C ----- Cavity problems ---
C
IF ( NPB.EQ.7.OR.NPB.EQ.6)
1 CALL VORB ( WP1, PS1, WPT, IX, IY, NNP, NEL )
C
IF ( NPB.EQ.7.OR.NPB.EQ.6) GO TO 900
C
C ----------------------------------------------------------
DO 100 J = 1, IY
DO 100 I = 1, IX
C ----------------------------------------------------------
C
IF ( WPT( I,J,1)) GO TO 100 ! Fluid poinT( B.)
IF ( WPT( I,J,2)) GO TO 50 ! SLP ( No-SLP)
IF ( WPT( I,J,3)) GO TO 150 ! Moving ( Not m.)
C
C ----- WPT( I,J,4) --- Continuative or ( Periodic ) ---
C
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 ------------------------------------------------
IF ( NPB.EQ.5.OR.NPB.EQ.4) THEN
C
IF ( J.EQ.JS .AND. I.GE.IS .AND. I.LE.IE )
C
1 WP1( I,J) = 0.5D0* ( 7.D0* PS1( I,J)
2 - 8.D0* PS1( I,J+1) + PS1( I,J+2))/ DY2
C
IF ( I.EQ.IE.AND.J.GT.1 .AND.J.LT.JS )
C
1 WP1( I,J) = 0.5D0* ( 7.D0* PS1( I,J)
2 - 8.D0* PS1( I+1,J) + PS1( I+2,J))/ DX2
C
IF ( J.EQ.1.AND.I.EQ.IE ) WP1( I,J) = 0.D0
GO TO 100
C
END IF
C ------------------------------------------------
C
C ----- No-SLP surfaces ---
C
IF ( .NOT. WPT( I+1,J,1)) GO TO 1
WA =( PS1( I,J) - PS1( I+1,J))/ DX2
C
13 IF ( .NOT. WPT( I-1,J,1)) GO TO 3
WB =( PS1( I,J) - PS1( I-1,J))/ DX2
C
12 IF ( .NOT. WPT( I,J+1, 1)) GO TO 2
WC = ( PS1( I,J) - PS1( I,J+1))/ DY2
C
14 IF (.NOT. WPT( I,J-1, 1)) GO TO 4
WD =( PS1( I,J) - PS1( I,J-1))/ DY2
110 WP1( I,J) = WA + WB + WC + WD
GO TO 100
C
1 WA = ( PS1( I,J) - PS1( I-1,J))/ DX2
GO TO 13
C
3 WB = ( PS1( I,J) - PS1( I+1,J))/ DX2
GO TO 12
C
2 WC = ( PS1( I,J) - PS1( I,J-1))/ DY2
GO TO 14
C
4 WD = ( PS1( I,J) - PS1( I,J+1))/ DY2
GO TO 110
C
C ----- No-SLP surfaces at edges ---
C
10 WP1( I,J) = 2.D0* ( PS1( I,J) - PS1( I,J-1))/ DY2
GO TO 100
C
20 WP1( I,J) = 2.D0* ( PS1( I,J) - PS1( I,J+1))/ DY2
GO TO 100
C
50 WP1( I,J) = 0.D0
GO TO 100
C
30 WA = ( PS1( I,J) - PS1( I-1,J))/ DX2
WB = WA
GO TO 12
C
40 WA = ( PS1( I,J) - PS1( I+1,J))/ DX2
WB = WA
GO TO 12
C
150 CONTINUE ! Moving
C
IF ( J.EQ.IY) GO TO 10
IF ( J.EQ. 1) GO TO 20
C
100 CONTINUE
C ----------------------------------------------------------
C
900 CONTINUE
C
IJ = 0
C ------------------------------------------------
DO 200 J = 1, IY
DO 200 I = 1, IX
C ------------------------------------------------
C
IJ = IJ + 1
PS2( IJ) = PS1( I,J)
WP2( IJ) = WP1( I,J)
WM1( I,J) = WP1( I,J)
C
200 CONTINUE
C ------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE WUC2 ( U, WP1, WM1, PS1, WPT, IX, IY, NNP, NEL )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION WP1( IX,IY), WM1( IX,IY), PS1( IX,IY), U( IX,IY)
C
LOGICAL*1 WPT( IX,IY,4)
C
COMMON /CNS/ NPB, WEI, BV, RV, KCR, OMW
COMMON /IDT/ DT, LMT, UL, VU, RD, VIS, REN, LOP
COMMON /UDT/ DX, DY, DX2, DY2, RS, RS2, CF, DTX, DTY
C
C ----------------------------------------------------------
DO 100 J = 1, IY
DO 100 I = 1, IX
C ----------------------------------------------------------
IF ( NPB.EQ.7.OR.NPB.EQ.6) THEN
C
IF ( I.EQ.IX) GO TO 100
C ----------------------------
IF ( J.EQ.IY) THEN
UF = VU
GO TO 151
END IF
C -----------------------------
ELSE
C -----------------------------
IF ( I.EQ.IX) THEN
C
IF ( J.EQ.1.OR.J.EQ.IY) GO TO 100
C
U( I,J) = ( PS1( I,J+1) - PS1( I,J-1))/( 2.D0* DY)
C =======================================================
GO TO 100
END IF
C ---------------------------
C
C ---------------------------
IF (.NOT.WPT( I,J,1)) THEN
UF = 0.D0 ! Boundary
GO TO 151
END IF
C ---------------------------
C
END IF
C ----------------------------------------------------------
C
U( I,J) = ( PS1( I,J+1) - PS1( I,J-1))/( 2.D0* DY)
C =======================================================
C
UF = ( PS1( I,J+1) - PS1( I,J-1) + PS1( I+1,J+1) - PS1( I+1,J-1))
1 /( 4.D0* DY)
C
151 CONTINUE
C
ALF = UF* DTX
C ==================
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 = ( WM1( I-1,J) + WM1( I,J))* OMW
TRM = ( WM1( I-1,J) - WM1( I,J))* OMW
C
TPL = ( WM1( I,J) + WM1( I+1,J))* WEI
TML = ( WM1( I,J) - WM1( 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
130 IF ( WPT( I,J,4)) GO TO 170
C
K = I + 1
A1 = 0.5D0* ALF
C ====================
C
A2 = A1* ALF
B2 = A2 + 2.D0* A1
C
TRP = ( WM1( K,J) + WM1( K+1,J))* OMW
TRM = ( WM1( K,J) - WM1( K+1,J))* OMW
C
TPL = ( WM1( K-1,J) + WM1( K,J))* WEI
TML = ( WM1( K-1,J) - WM1( K,J))* WEI
C
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
GO TO 100
C
150 TF = ALF* WM1( I,J)
GO TO 180
C
170 TF = ALF* WM1( 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 WVC2 ( V, WP1, WM1, PS1, WPT, EPT, IX, IY, NNP, NEL )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION WP1( IX,IY), WM1( IX,IY), PS1( IX,IY), V( IX,IY)
C
LOGICAL*1 WPT( IX,IY,4), EPT( IX,IY,4)
C
COMMON /CNS/ NPB, WEI, BV, RV, KCR, OMW
COMMON /IDT/ DT, LMT, UL, VU, RD, VIS, REN, LOP
COMMON /UDT/ DX, DY, DX2, DY2, RS, RS2, CF, DTX, DTY
C
C -----------------------------------------------------------
DO 100 J = 1, IY
DO 100 I = 1, IX
C ----------------------------------------------------------
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
V( I,J) = ( PS1( I-1,J) - PS1( I+1,J))/( 2.D0* DX )
C ========================================================
C
VF = ( PS1( I-1,J) - PS1( I+1,J) + PS1( I-1,J+1) - PS1( I+1,J+1))
1 /( 4.D0* DX )
C
C ------------------------------------------------
IF ( NPB.EQ.5.OR. NPB.EQ.4) THEN
C
IF ( .NOT.WPT( I,J,1)) THEN
V( I,J) = 0.D0 ! Boundary point
VF = 0.D0 ! Boundary point
END IF
C
END IF
C ------------------------------------------------
C
BET = VF* DTY
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 CONTINUE
A1 = 0.5D0* BET
A2 = A1* BET
B2 = A2 - 2.D0* A1
C
TRP = ( WM1( I,J-1) + WM1( I,J)) * OMW
TRM = ( WM1( I,J-1) - WM1( I,J)) * OMW
TPL = ( WM1( I,J ) + WM1( I,J+1))* WEI
TML = ( WM1( I,J ) - WM1( I,J+1))* WEI
C
TF = A1* TPL + A1* TRP + A2* TML + B2* TRM
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 332
K = J + 1
C
A1 = 0.5D0* BET
C =====================
A2 = A1* BET
B2 = A2 + 2.D0* A1
C
TRP = ( WM1( I,K ) + WM1( I,K+1))* OMW
TRM = ( WM1( I,K ) - WM1( I,K+1))* OMW
TPL = ( WM1( I,K-1) + WM1( I,K )) * WEI
TML = ( WM1( I,K-1) - WM1( 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* WM1( I,J)
GO TO 333
C
332 TF = BET* WM1( I,J+1)
C
333 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 WUC4 ( U, WP1, WM1, PS1, WPT, IX, IY )
C
IMPLICIT REAL*8 (A-H,O-Z )
C
DIMENSION WP1( IX,IY), PS1( IX,IY), WM1( IX,IY), U( IX,IY)
C
LOGICAL*1 WPT( IX,IY,4)
C
COMMON /CNS/ NPB, WEI, BV, RV, KCR, OMW
COMMON /IDT/ DT, LMT, UL, VU, RD, VIS, REN, LOP
COMMON /UDT/ DX, DY, DX2, DY2, RS, RS2, CF, DTX, DTY
C
COMMON /ABCD/ A1, A2, A3, A4, B1, B2, B3, B4, C1, C2, C3, C4,
1 D1, D2, D3, D4
C
C ----------------------------------------------------------
DO 100 J = 1, IY
DO 100 I = 1, IX
C ----------------------------------------------------------
C
C ------------------------------------------------
IF ( NPB.EQ.7.OR.NPB.EQ.6) THEN
C
IF ( I.EQ.IX) GO TO 100
C
C --------------------------------------
IF ( J.EQ.IY) THEN
UF = VU
GO TO 150
END IF
C --------------------------------------
IF ( I.EQ.1) THEN
UF = 0.D0
GO TO 150
END IF
C --------------------------------------
END IF
C --------------------------------------
IF ( I.EQ.IX) THEN
IF ( J.EQ.1.OR.J.EQ.IY) GO TO 100
C
U( I,J) = ( PS1( I,J+1) - PS1( I,J-1))/( 2.D0* DY)
GO TO 100
END IF
C --------------------------------------
IF ( .NOT.WPT( I,J,1)) THEN
UF = 0.D0 ! Boundary
GOTO 150
END IF
C --------------------------------------
C
U( I,J) = ( PS1( I,J+1) - PS1( I,J-1))/( 2.D0* DY)
C =======================================================
C
UF = ( PS1( I,J+1) - PS1( I,J-1) + PS1( I+1,J+1) - PS1( I+1,J-1))
1 /( 4.D0* DY)
C ===== Attention ========
C
150 CONTINUE
C
ALF = UF* DT/ DX
C =====================
C
IF ( ALF.LT. 0.D0) GO TO 311
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
C
C --------------------------------------
C
CALL CFSB ( ALF, 0)
C
C --------------------------------------
C
IF ( I .EQ. 1) GO TO 304
IF ( I-1.EQ. 1) GO TO 305
IF ( I+1.EQ.IX) GO TO 309
C
FDP = ( WM1( I-2,J) + WM1( I+1,J))* OMW
FDM = ( WM1( I-2,J) - WM1( I+1,J))* OMW
308 FBP = ( WM1( I-1,J) + WM1( I+2,J))* WEI
FBM = ( WM1( I-1,J) - WM1( I+2,J))* WEI
C
307 FCP = ( WM1( I-1,J) + WM1( I,J)) * OMW
FCM = ( WM1( I-1,J) - WM1( I,J)) * OMW
C
306 FAP = ( WM1( I,J) + WM1( I+1,J))* WEI
FAM = ( WM1( I,J) - WM1( 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 302
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
302 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
304 FDP = ( WM1( I-2+( IX-1),J) + WM1( I+1,J))* OMW
FDM = ( WM1( I-2+( IX-1),J) - WM1( I+1,J))* OMW
FBP = ( WM1( I-1+( IX-1),J) + WM1( I+2,J))* WEI
C
FBM = ( WM1( I-1+( IX-1),J) - WM1( I+2,J))* WEI
FCP = ( WM1( I-1+( IX-1),J) + WM1( I,J)) * OMW
FCM = ( WM1( I-1+( IX-1),J) - WM1( I,J)) * OMW
GO TO 306
C
305 FDP = ( WM1( I-2+( IX-1),J) + WM1( I+1,J))* OMW
FDM = ( WM1( I-2+( IX-1),J) - WM1( I+1,J))* OMW
GO TO 308
C
309 FBP = ( WM1( I-1,J) + WM1( I+2-( IX-1),J))* WEI
FBM = ( WM1( I-1,J) - WM1( I+2-( IX-1),J))* WEI
FDP = ( WM1( I-2,J) + WM1( I+1,J))* OMW
FDM = ( WM1( I-2,J) - WM1( I+1,J))* OMW
GO TO 307
C
311 IF ( .NOT. WPT( I+1,J,1)) GO TO 332
IF ( WPT( I+1,J,4)) GO TO 332
IF ( WPT( I, J,4)) GO TO 332
C
K = I + 1
C
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
C
IF ( .NOT. WPT( K-1,J,1)) GO TO 201
IF ( WPT( K-1,J,4)) GO TO 201
C
401 CONTINUE
C
C --------------------------------------
C
CALL CFSB ( ALF, 1 )
C
C --------------------------------------
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 = ( WM1( K-1,J) + WM1( K+2,J))* OMW
FDM = ( WM1( K-1,J) - WM1( K+2,J))* OMW
318 FBP = ( WM1( K-2,J) + WM1( K+1,J))* WEI
FBM = ( WM1( K-2,J) - WM1( K+1,J))* WEI
C
317 FCP = ( WM1( K,J) + WM1( K+1,J))* OMW
FCM = ( WM1( K,J) - WM1( K+1,J))* OMW
316 FAP = ( WM1( K-1,J) + WM1( K,J)) * WEI
FAM = ( WM1( K-1,J) - WM1( 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
C
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 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 = ( WM1( K-1,J) + WM1( K+2 -( IX-1),J))* OMW
FDM = ( WM1( K-1,J) - WM1( K+2 -( IX-1),J))* OMW
FBP = ( WM1( K-2,J) + WM1( K+1 -( IX-1),J))* WEI
C
FBM = ( WM1( K-2,J) - WM1( K+1 -( IX-1),J))* WEI
FCP = ( WM1( K,J) + WM1( K+1 -( IX-1),J))* OMW
FCM = ( WM1( K,J) - WM1( K+1 -( IX-1),J))* OMW
GO TO 316
C
315 FDP = ( WM1( K-1,J) + WM1( K+2 -( IX-1),J))* OMW
FDM = ( WM1( K-1,J) - WM1( K+2 -( IX-1),J))* OMW
GO TO 318
C
319 FBP = ( WM1( K-2+( IX-1),J) + WM1( K+1,J))* WEI
FBM = ( WM1( K-2+( IX-1),J) - WM1( K+1,J))* WEI
FDP = ( WM1( K-1,J) + WM1( K+2,J))* OMW
FDM = ( WM1( K-1,J) - WM1( K+2,J))* OMW
GO TO 317
C
600 TF = ALF* WM1( I,J)
GO TO 333
C
332 TF = ALF* WM1( I+1,J)
C
333 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) = WM1( I+1,J) + ALF* ( WM1( I,J) - WM1( I+1,J))
GO TO 100
C
200 A1 = 0.5D0* ALF
A2 = A1* ALF
B2 = A2 - 2.D0* A1
C
TRP = ( WM1( I-1,J) + WM1( I,J)) * OMW
TRM = ( WM1( I-1,J) - WM1( I,J)) * OMW
TPL = ( WM1( I,J) + WM1( I+1,J))* WEI
TML = ( WM1( I,J) - WM1( 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 = ( WM1( K,J) + WM1( K+1,J))* OMW
TRM = ( WM1( K,J) - WM1( K+1,J))* OMW
TPL = ( WM1( K-1,J) + WM1( K,J)) * WEI
TML = ( WM1( K-1,J) - WM1( 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
C
100 CONTINUE
C ----------------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE WVC4 ( V, WP1, WM1, PS1, WPT, EPT, IX, IY )
C
IMPLICIT REAL*8 (A-H,O-Z )
C
DIMENSION WP1( IX,IY), PS1( IX,IY), WM1( IX,IY), V( IX,IY)
C
LOGICAL*1 WPT( IX,IY,4), EPT( IX,IY,4)
C
COMMON /CNS/ NPB, WEI, BV, RV, KCR, OMW
COMMON /IDT/ DT, LMT, UL, VU, RD, VIS, REN, LOP
COMMON /UDT/ DX, DY, DX2, DY2, RS, RS2, CF, DTX, DTY
C
COMMON /ABCD/ A1, A2, A3, A4, B1, B2, B3, B4, C1, C2, C3, C4,
1 D1, D2, D3, D4
C
C ----------------------------------------------------------
DO 100 J = 1, IY
DO 100 I = 1, IX
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
V( I,J) = ( PS1( I-1,J) - PS1( I+1,J))/( 2.D0* DX)
C
VF = ( PS1( I-1,J) - PS1( I+1,J) + PS1( I-1,J+1) - PS1( I+1,J+1))
1 /( 4.D0* DX)
C ****** Attention *******
C
C --------------------------------------
IF ( NPB.EQ.5.OR.NPB.EQ.4) THEN
C
IF ( .NOT.WPT( I,J,1)) THEN
C
V( I,J) = 0.D0 ! Boundary point
VF = 0.D0 ! Boundary point
C
END IF
C
END IF
C --------------------------------------
C
BET = VF* DT/ DY
C =====================
C
IF ( BET .LT. 0.D0) GO TO 311
IF ( .NOT. WPT( I,J,1)) GO TO 600
IF ( EPT( I,J,2)) GO TO 400
C
IF ( .NOT. WPT( I,J-1,1)) GO TO 200
IF ( .NOT. WPT( I,J+1,1)) GO TO 200
C
400 CONTINUE
C
C --------------------------------------
C
CALL CFSB ( BET, 0 )
C
C --------------------------------------
C
FDP = ( WM1( I,J-2) + WM1( I,J+1))* OMW
FDM = ( WM1( I,J-2) - WM1( I,J+1))* OMW
FBP = ( WM1( I,J-1) + WM1( I,J+2))* WEI
FBM = ( WM1( I,J-1) - WM1( I,J+2))* WEI
C
FCP = ( WM1( I,J-1) + WM1( I,J)) * OMW
FCM = ( WM1( I,J-1) - WM1( I,J)) * OMW
FAP = ( WM1( I,J) + WM1( I,J+1))* WEI
FAM = ( WM1( I,J) - WM1( 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
C
GO TO 100
C
502 CONTINUE
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 = ( PS1( I-1,J) - PS1( I+1-( IX-1),J)
1 + PS1( I-1,J+1) - PS1( I+1-( IX-1),J+1))/( 4.D0* DX)
C
GO TO 505
C
512 CONTINUE
C
IF ( WPT( I,J+1,4).OR.WPT( I,J,4)) GO TO 100 ! Continuative
C
IF (.NOT.WPT( I,J+1,1).AND..NOT.WPT( I,J,1)) GO TO 100
C
C ----- Periodic Left ---
C
VF = ( PS1( I-1+( IX-1),J) - PS1( I+1,J) + PS1( I-1+( IX-1),J+1)
1 - PS1( I+1,J+1))/( 4.D0* DX)
C
GO TO 505
C
202 VF = ( PS1( I-1,J) - PS1( I,J) + PS1( I-1,J+1) - PS1( I,J+1))
1 /( 2.D0* DX)
C
V( I,J) = ( PS1( I-1,J) - PS1( I,J))/ DX
C =============================================
C
505 CONTINUE
C
BET = VF* DT/ DY
C =====================
C
IF ( BET.LT.0.D0) GO TO 332
GO TO 600
C
311 IF ( .NOT.WPT( I,J+1,1)) GO TO 332
C
K = J + 1
C
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 CONTINUE
C
C --------------------------------------
C
CALL CFSB ( BET,1 )
C
C --------------------------------------
C
FDP = ( WM1( I,K-1) + WM1( I,K+2))* OMW
FDM = ( WM1( I,K-1) - WM1( I,K+2))* OMW
FBP = ( WM1( I,K-2) + WM1( I,K+1))* WEI
FBM = ( WM1( I,K-2) - WM1( I,K+1))* WEI
C
FCP = ( WM1( I,K) + WM1( I,K+1))* OMW
FCM = ( WM1( I,K) - WM1( I,K+1))* OMW
FAP = ( WM1( I,K-1) + WM1( I,K)) * WEI
FAM = ( WM1( I,K-1) - WM1( 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* WM1( I,J)
GO TO 333
C
332 TF = BET* WM1( I,J+1)
C
333 CONTINUE
WP1( I,J) = WP1( I,J) - TF
WP1( I,J+1) = WP1( I,J+1) + TF
GO TO 100
C
200 A1 = 0.5D0* BET
A2 = A1* BET
B2 = A2 - 2.D0* A1
C
TRP = ( WM1( I,J-1) + WM1( I,J)) * OMW
TRM = ( WM1( I,J-1) - WM1( I,J)) * OMW
TPL = ( WM1( I,J) + WM1( I,J+1))* WEI
TML = ( WM1( I,J) - WM1( 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 = ( WM1( I,K) + WM1( I,K+1))* OMW
TRM = ( WM1( I,K) - WM1( I,K+1))* OMW
TPL = ( WM1( I,K-1) + WM1( I,K)) * WEI
TML = ( WM1( I,K-1) - WM1( 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 ----------------------------------------------------------
C
RETURN
END
C **********************************************************************