戻る

C
PROGRAM MAIN
C ******************************************************************** C
C C
C F 2 D - B S G C
C C
C ( Fourth ordered Finite Difference Method ) C
C C
C ( Version 3.8 ) C
C C
C Data: Equal length C
C C
C ( SW: DX.ne.DY / Variable length meshes: ) C
C C
C Copyright : Yasuhiro MATSUDA C
C C
C ******************************************************************** C
C
IMPLICIT REAL*8 (A-H,O-Z)
C
INTEGER*2 NDE
C
C ***** DATA ( 1 ) ***************************************************
C
*** PARAMETER ( IX= 561, IY= 101,
C
PARAMETER ( IX= 181, IY= 61,
C
*** PARAMETER ( IX = 45, IY = 15,
C
1 IX1 = IX-1, IY1 = IY-1, IXY = IX* IY, MSH = IX1* IY1 )
C
C ----- IXY: NNP & MSH ---
C
LOGICAL*1 WPT( IX,IY,4), PPT( IX,IY,4), TPT( IX,IY,4),
1 EPT( IX,IY,4), CPT( IX,IY,4), PET( IXY,4),
2 EET( IXY,4)
C
DIMENSION XE ( IXY ), YE ( IXY), PS1( IX,IY), PS2( IXY),
1 PS3( IY), PS ( IX,IY), W ( IX,IY), WP1( IX,IY),
2 WP2( IXY), FLW( IY), XS ( IXY), YS ( IXY),
3 FF ( IX,IY), WB ( IX,IY), WD ( IX,IY), CP1( IX,IY),
4 U ( IX,IY), V ( IX,IY), UB ( IX,IY), VB ( IX,IY),
5 HH ( IX,IY)
C
DIMENSION BX ( IX,IY), BY ( IX,IY), RX( IX,IY), RY ( IX,IY),
1 XZ ( IX), YZ ( IY), RSS( IX1), RAS( IY1),
2 UE ( IXY), VE ( IXY), UN ( IX,IY), VN ( IX,IY),
3 XX ( IXY), YY ( IXY), NDE( MSH,4)
C
DIMENSION PA1( IX,IY), PA2( IX,IY), PA3( IX,IY), PB1( IX,IY),
1 PB2( IX,IY), PB3( IX,IY)
C
DIMENSION TX1( IX), TX2( IX), TX3( IX), TY1( IY),
1 TY2( IY), TY3( IY)
C
DIMENSION UNA( IX,4), UNB( IX,4), UNC( IX,4), UND( IX,4),
1 UNE( IX,4), UAP( IX,4), UBP( IX,4), UCP( IX,4),
2 UDP( IX,4), UEP( IX,4), UAM( IX,4), UBM( IX,4),
3 UCM( IX,4), UDM( IX,4), UEM( IX,4), UF1( IY),
4 UF2( IY), UF3( IY)
C
DIMENSION VNA( IY,4), VNB( IY,4), VNC( IY,4), VND( IY,4),
1 VNE( IY,4), VAP( IY,4), VBP( IY,4), VCP( IY,4),
2 VDP( IY,4), VEP( IY,4), VAM( IY,4), VBM( IY,4),
3 VCM( IY,4), VDM( IY,4), VEM( IY,4), VF1( IX),
4 VF2( IX), VF3( IX)
C
C ----- Each mesh size ---
C
DIMENSION HX( IX1), HY( IY1)
C
C ----- Attention --- No. of Var. of Vorticity Data = < 100,000 ---
C
DIMENSION X1R( 100000)
C
C ----- Attention --- No. of Var. of Vorticity Data = < 100,000 ---
C
COMMON /WEH/ WEI, IHS, LTC
COMMON /CTL/ LOP, NPR, ILP, NNN
LOGICAL PRS
COMMON /IDT/ DT, VIS, DX, DY, DX2, DY2, ITT, CON, PRS, RLL
C
333 CONTINUE
C
CALL TTLE
C
C ***** DATA ( 2 ) ***************************************************
C
C --- For 560 X 100 meshes ---
C
** OPEN ( 1, FILE='560100_BSG.DATA' ) ! L = 35, RE = 1000
C
** OPEN ( 7, FILE='BSG_XRI_R1000.DATA' )
** OPEN ( 8, FILE='BSG_OUT_R1000.DATA' )
** OPEN (17, FILE='BSG_VXR_R1000.DATA' )
C
C --- For 180 X 60 meshes ---
C
OPEN ( 1, FILE='18060_BSG.DATA' ) ! L = 35, RE = 800, 100
C
OPEN ( 7, FILE='BSG_XRI_R800.DATA' )
OPEN ( 8, FILE='BSG_OUT_R800.DATA' )
OPEN (17, FILE='BSG_VXR_R800.DATA' )
C
** OPEN ( 7, FILE='BSG_XRI_R100.DATA' )
** OPEN ( 8, FILE='BSG_OUT_R100.DATA' )
** OPEN (17, FILE='BSG_VXR_R100.DATA' )
C
C --- For 44 X 14 meshes ---
C
** OPEN ( 1, FILE='4414_BSG.DATA' ) ! L = 35, RE = 100
C
C -------------------------------------------
C
** OPEN ( 7, FILE='BSG_XRI_R100.DATA' )
** OPEN ( 8, FILE='BSG_OUT_R100.DATA' )
** OPEN (17, FILE='BSG_VXR_R100.DATA' )
C
C --------------------------------------------------------------------
C
CALL INIT ( PPT, WPT, TPT, CPT, EPT, PS1, PS3, WP1, CP1,
1 PS, FF, W, U, V, UB, VB, UN, VN,
2 HH, BX, BY, RX, RY, XZ, YZ, NDE, HX,
3 HY, IX, IY )
C
CALL SFCF ( PA1, PA2, PA3, PB1, PB2, PB3, HX, HY, IX, IY )
C
CALL DFCF ( TX1, TX2, TX3, TY1, TY2, TY3, HX, HY, IX, IY )
C
CALL WUCF ( UNA, UNB, UNC, UND, UNE, UAP, UBP, UCP, UDP,
1 UEP, UAM, UBM, UCM, UDM, UEM, UF1, UF2, UF3,
2 HX, HY, IX, IY )
C
CALL WVCF ( VNA, VNB, VNC, VND, VNE, VAP, VBP, VCP, VDP,
1 VEP, VAM, VBM, VCM, VDM, VEM, VF1, VF2, VF3,
2 HX, HY, IX, IY )
C
CALL DTIN ( HX, HY, XZ, RSS, RAS, KZ, IX, IY, IX1,
1 IY1 )
C
CALL GEIN ( XZ, RSS, RAS, KZ, WP1, PS1, CP1, WPT, TPT,
1 PPT, CPT, EPT, IX, IY, IX1, IY1 )
C
CALL PRGS ( PS1, PS3, PPT, FLW, XX, YY, NDE, XZ, YZ,
1 XE, YE, XS, YS, IMP, JMP, IS, JS, IE,
2 IX, IY, IXY, MSH )
C
NNN = 0
LOP = 0
ICT = 0
C --------------------------------------------------------------------
999 LOP = LOP + 1
C
C ----- Stream Function ---
C
CALL SFCL ( PS1, PS, W, WP1, PPT, EPT, 1.5D0, ICT, YZ,
1 PA1, PA2, PA3, PB1, PB2, PB3, HX, HY, IX,
2 IY, IXY )
C
C ----- Vorticity at Boundary ---
C
CALL WFCA ( WP1, PS1, WPT, RAS, RSS, IX, IY )
C
CALL PRPS ( PS1, PS2, XE, YE, IE, JS, WP1, WP2, UB,
1 VB, U, V, FF, IX, IY, IXY )
C
CALL WUCL ( WP1, FF, PS1, WPT, EPT, U, YZ, IX, IY,
1 UNA, UNB, UNC, UND, UNE, UAP, UBP, UCP, UDP,
2 UEP, UAM, UBM, UCM, UDM, UEM, UF1, UF2, UF3,
3 HX, HY )
C
CALL WVCL ( WP1, FF, PS1, WPT, EPT, V, IX, IY, VNA,
1 VNB, VNC, VND, VNE, VAP, VBP, VCP, VDP, VEP,
2 VAM, VBM, VCM, VDM, VEM, VF1, VF2, VF3, HX,
3 HY )
C
CALL DCLC ( WP1, FF, WPT, EPT, U, V, HH, BX, BY,
1 RX, RY, YZ, TX1, TX2, TX3, TY1, TY2, TY3,
2 HX, HY, IX, IY )
C
CALL VRVR ( WP1, WB, WD, II, JJ, W, PS, X1R, PS1,
1 IX, IY, *800 )
C
GO TO 999
C
800 CONTINUE
C
CALL OUTD ( U, V, UE, VE, PPT, EPT, WP2, IS, JS,
1 IE, XS, YS, PET, EET, XE, YE, PS2, NDE,
2 IX, IY, IXY, IX1, IY1, MSH, KEY )
C
IF ( KEY.EQ.1) GO TO 333
C
CLOSE ( 1)
CLOSE ( 7)
CLOSE ( 8)
CLOSE (17)
C
STOP
END
C *********************************************************************
C
SUBROUTINE DCLC ( WP1, FF, TPT, EPT, U, V, HH, BX, BY,
1 RX, RY, YZ, TX1, TX2, TX3, TY1, TY2, TY3,
2 HX, HY, IX, IY )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
COMMON /CTL/ LOP, NPR, ILP, NNN
C
DIMENSION WP1( IX,IY), FF( IX,IY), U ( IX,IY), V ( IX,IY),
1 HH ( IX,IY), BX( IX,IY), BY( IX,IY), RX( IX,IY),
2 RY ( IX,IY), YZ( IY)
C
LOGICAL*1 TPT( IX,IY,4), EPT( IX,IY,4)
C
COMMON /WEH/ WEI, IHS, LTC
LOGICAL PRS
COMMON /IDT/ DT, VIS, DX, DY, DX2, DY2, ITT, CON, PRS, RLL
C
DIMENSION TX1( IX), TX2( IX), TX3( IX), TY1( IY), TY2( IY),
1 TY3( IY), HX ( IX-1), HY( IY-1)
C
C ----- For Diffusion term ---
C
GO TO ( 301, 302, 304), LTC
C
301 RH = CON
C -------------
GO TO 303
C
C ----- Attention ---
302 RH = 1.D0
GO TO 303
C
304 RH = 4.D0* VIS
303 CONTINUE
C
IF ( IHS.EQ.1) GO TO 81
C ----------------------------
DO 100 J = 1, IY
DO 100 I = 1, IX
C ----------------------------
C
HH( I,J) = RH* DT
100 CONTINUE
C ----------------------------
GO TO 85
C
81 CONTINUE
C
C ----- Correction coefficient ---
C
C --------------------------------------
DO 200 J = 1, IY-1
DO 200 I = 1, IX-1
C --------------------------------------
C
BX( I,J) = U( I,J)* DT/ HX( I)
BY( I,J) = V( I,J)* DT/ HY( J)
C
RX( I,J) = RH* DT/( HX( I)** 2)
RY( I,J) = RH* DT/( HY( J)** 2)
C
HH( I,J) = ( 1.D0 + ( BX( I,J)* BY( I,J))
1 /( RX( I,J) + RY( I,J)))* RH * DT
C
200 CONTINUE
C --------------------------------------
85 CONTINUE
C
C --------------------------------------
DO 300 J = 1, IY
DO 300 I = 1, IX
C
IF ( TPT( I,J,1)) GO TO 300
WP1( I,J) = FF( I,J)
C
300 CONTINUE
C --------------------------------------
DO 400 J = 1, IY
DO 400 I = 1, IX
C
IF ( TPT( I,J,1)) GO TO 1
IF ( TPT( I,J,2)) GO TO 400
C
C ----- Not Fluid, Not Conducting, Not Insulated ---
C
GO TO 400
C
1 CONTINUE
IF ( EPT( I,J,3)) GO TO 2 ! Right
IF ( EPT( I,J,4)) GO TO 12 ! Left
C
WP1( I,J) = WP1( I,J) + 2.D0 * HH( I,J)* ( TX1( I)* FF( I-1,J)
1 - TX2( I)* FF( I,J) + TX3( I)* FF( I+1,J)
C
2 + TY1( J)* FF( I,J-1) - TY2( J)* FF( I,J) + TY3( J)* FF( I,J+1))
GO TO 400
C
2 CONTINUE
WP1( I,J) = WP1( I-1,J) ! ( Outflow )
GO TO 400
C
C ---- Periodic Right ---
C
12 CONTINUE
YTR = ( YZ( J) - 1.D0)/ 2.D0
WP1( I,J) = 12.D0* ( 4.D0* YTR - 1.D0 )
C
400 CONTINUE
C --------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE DFCF ( TX1, TX2, TX3, TY1, TY2, TY3, HX, HY, IX, IY )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION TX1( IX), TX2( IX), TX3( IX), TY1(IY), TY2(IY),
1 TY3( IY), HX ( IX-1), HY( IY-1)
C
C -----------------------
DO 100 I = 1, IX
C
TX1( I) = 0.D0
TX2( I) = 0.D0
TX3( I) = 0.D0
C
100 CONTINUE
C -----------------------
DO 200 J = 1, IY
C
TY1( J) = 0.D0
TY2( J) = 0.D0
TY3( J) = 0.D0
C
200 CONTINUE
C ----------------------------------------------------------
DO 300 I = 2, IX-1
C
TX1( I) = 1.D0/( HX( I-1)* ( HX( I-1) + HX( I)))
TX2( I) = 1.D0/( HX( I-1)* HX( I))
TX3( I) = 1.D0/( HX( I) * ( HX( I-1) + HX( I)))
C
300 CONTINUE
C ----------------------------------------------------------
DO 400 J = 2, IY-1
C
TY1( J) = 1.D0/( HY( J-1)* ( HY( J-1) + HY( J)))
TY2( J) = 1.D0/( HY( J-1)* HY( J))
TY3( J) = 1.D0/( HY( J)* ( HY( J-1) + HY( J)))
C
400 CONTINUE
C ----------------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE DIST ( PS1, X, Y, IE, JS, IX, IY )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION PS1( IX,IY), X( IX,IY), Y( IX,IY)
C
COMMON /CTL/ LOP, NPR, ILP, NNN
COMMON /DDS/ DS1, DS2, DS3, DSV( 10), DIS( 100)
LOGICAL PRS
COMMON /IDT/ DT, VIS, DX, DY, DX2, DY2, ITT, CON, PRS, RLL
C
C ----- X1, X4 , X5 ---
C
C ------------------------------------
DO 100 I = IE, IX-1
C
IF ( PS1(I,2).GE.-1.D0 .AND. PS1(I+1,2).LE.-1.D0 ) THEN
DIS( 2) = ( X(I,2) + X(I+1,2))/ 2.D0
C
IF ( DABS( PS1(I,2) - PS1(I+1,2)).GT.0.0000001D0 )
1DIS( 2) = X(I,2) + ( -1.D0 - PS1(I,2))/
2 ( PS1(I+1,2) - PS1(I,2))* ( X(I+1,2) - X(I,2))
DIS( 2) = ( DIS( 2) - X(IE,2))/( Y(1,JS) - Y(1,1))
C
END IF
C ----------------------
C
IF ( PS1(I,2).LE.-1.D0 .AND. PS1(I+1,2).GE.-1.D0 ) THEN
C
DIS( 3) = DIS( 1)
DIS( 1) = ( X(I,2) + X(I+1,2))/ 2.D0
IF ( DABS( PS1( I,2) - PS1( I+1,2)).GT.0.0000001D0 )
1DIS( 1) = X(I,2) + ( -1.D0 - PS1(I,2))/
2 ( PS1(I+1,2) - PS1(I,2))* ( X(I+1,2) - X(I,2))
DIS( 1) = ( DIS(1) - X(IE,2))/( Y(1,JS ) - Y(1,1))
C
END IF
C
100 CONTINUE
C ---------------------------------
C
C ***** X 2 , X 3 ***
C
DO 150 I = IE, IX-1
C
C ---------------------------
C
IF ( PS1( I,IY-1).LE.1.D0.AND.PS1( I+1,IY-1).GE.1.D0 ) THEN
C
DSV( 1) = ( X(I,IY-1) + X(I+1,IY-1))/ 2.D0
IF ( DABS( PS1(I,IY-1) - PS1(I+1,IY-1)).GT. 1.D-7 )
1 DSV( 1)= X(I,IY-1) + ( 1.D0 - PS1(I,IY-1))/
2 ( PS1(I+1,IY-1) - PS1(I,IY-1))* ( X(I+1,IY-1) - X(I,IY-1))
DSV( 1) = ( DSV(1) - X(IE,IY-1))/( Y(1,JS ) - Y(1,1))
C
END IF
C
C ---------------------------
C
IF ( PS1( I,IY-1).GE.1.D0.AND.PS1( I+1,IY-1).LE.1.D0 ) THEN
C
DSV( 2) = ( X( I,IY-1) + X(I+1,IY-1))/ 2.D0
IF ( DABS( PS1(I,IY-1) - PS1(I+1,IY-1)).GT.1.D-7 )
1 DSV( 2) = X(I, IY-1) + ( 1.D0 - PS1(I, IY-1))
2 /( PS1(I+1,IY-1) - PS1(I, IY-1))* ( X(I+1,IY-1) - X(I,IY-1))
DSV( 2) = ( DSV( 2) - X(IE,IY-1))/( Y ( 1,JS ) - Y( 1,1))
C
END IF
C
C ----------------------------
C
150 CONTINUE
C
C -------------------------------
C
IF ( LOP/ NPR* NPR .EQ. LOP ) THEN
C
WRITE(7,2100) LOP, DT
WRITE(7,2200) 1.D0/ VIS
2100 FORMAT(1X,'Loop =', I6,2X,'DT =',F4.3)
2200 FORMAT(2X,'Re =',F7.2)
C
DS1 = DIS( 1)
DS2 = DIS( 2)
DS3 = DIS( 3)
IF ( DABS( DS1 - DS3).LE.0.01D0) DS1 = 0.D0
C
WRITE(7,2300) DS3
WRITE(7,2400) DSV( 1)
WRITE(7,2500) DSV( 2)
WRITE(7,2600) DS2
WRITE(7,2700) DS1
2300 FORMAT(' X1R= ',F8.4)
2400 FORMAT(' X2R= ',F8.4)
2500 FORMAT(' X3R= ',F8.4)
2600 FORMAT(' X4R= ',F8.4)
2700 FORMAT(' X5R= ',F8.4)
C
END IF
C --------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE DTIN ( HX, HY, XZ, RSS, RAS, KZ, IX, IY, IX1, IY1 )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION XZ( IX), RSS( IX1), RAS( IY1)
C
COMMON /WEH/ WEI, IHS, LTC
COMMON /CTL/ LOP, NPR, ILP, NNN
LOGICAL PRS
COMMON /IDT/ DT, VIS, DX, DY, DX2, DY2, ITT, CON, PRS, RLL
C
DIMENSION HX( IX1), HY( IY1)
C
C ----- Data --- ( Ex. L = 35, Re = 800 ) ---
C
IHS = 1
C
WEI = 0.D0
C ===============
DT = 0.D0
C
CON = 0.01D0
PRS = .FALSE.
C
WRITE(6,2000)
2000 FORMAT(' ',' Re = ? ')
C
*** READ (5,*) RE
RE = 800.D0
** RE = 100.D0
C
VIS = 1.D0/ RE
C ===================
C
WRITE(6,2100)
2100 FORMAT(' ',' Time step = ? ')
** READ (5,*) DT
DT = 0.025D0
C
C DT = 0.01D0
C
WRITE(6,2200)
2200 FORMAT(' Loop = ? ')
C
*** READ(5,*) ILP
ILP = 100
** ILP = 1600
ILP = 2000
ILP = 50000
ILP = 100
C
WRITE(6,2300)
2300 FORMAT(' ',' What kind of Graphics ? ')
WRITE(6,*) '1:δ-ψ-ω 2:δ & Velocity 3:δ & ω'
*** READ (5,*) KZ
KZ = 1
C
WRITE(6,2400)
2400 FORMAT(' ','Intervals of Graphics output ? ')
*** READ (5,*) NPR
NPR = 20
NPR = 100
NPR = 500
NPR = 20
C
IXY = IX * IY
MSH = IX1* IY1
C
C ----- Attention --- Comp. size: RLL X 2.0 ---
C =======================
C
DX = RLL / DFLOAT( IX-1)
DY = 2.D0/ DFLOAT( IY-1)
C
DX2 = DX* DX
DY2 = DY* DY
C
C ---------------------------------------------
WRITE(6,2500) IX1, IY1, MSH
2500 FORMAT(/,' ** IX1 =', I4,' IY1 =',I4,' Meshes =',I8)
WRITE(6,2600) DX, DY
2600 FORMAT(' * DX =',F8.4,' * DY =',F8.4)
WRITE(6,2700) RE, DT, ILP, KZ, NPR
2700 FORMAT(/,' *** Re =', F8.2,' DT =',F6.3,' ILP=',I7,' KZ =',I2,
1 ' NPR =',I4,/)
WRITE(8,2800) RE, DT, ILP, KZ, NPR
2800 FORMAT(F9.2,F6.3, I6, I2, I3)
C
C -------------------------------
DO 100 I = 1, IX-1
C
RSS( I) = HX( I)* HX( I)
100 CONTINUE
C -------------------------------
DO 200 J = 1, IY-1
C
RAS( J) = HY( J)* HY( J)
200 CONTINUE
C --------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE GEIN ( XZ, RSS, RAS, KZ, WP1, PS1, PR1, WPT,
1 TPT, PPT, PRT, EPT, IX, IY, IX1, IY1 )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION INDX( 5), INDY( 5), INDE( 5)
C
LOGICAL*1 OMG, PSI, CNT, SLP, MOV, PRS, WPT, TPT, PPT, PRT,
1 EPT
C
DIMENSION XZ( IX), RSS( IX1), RAS( IY1)
C
DIMENSION WP1( IX,IY), PS1( IX,IY), PR1( IX,IY),
1 WPT( IX,IY,4), TPT( IX,IY,4), PPT( IX,IY,4),
2 PRT( IX,IY,4), EPT( IX,IY,4)
C
NAMELIST /DATS/ ICD, OMG, PSI, SLP, MOV, CNT, PRS
C
INER = 0
C
C -----------------------------
DO 100 J = 1, IY
DO 100 I = 1, IX
C -----------------------------
C
WPT( I,J,1) = .TRUE.
PPT( I,J,1) = .TRUE.
PRT( I,J,1) = .TRUE.
TPT( I,J,1) = .TRUE.
C
100 CONTINUE
C -----------------------------
DO 200 I = 1, IX
C
EPT( I, 1,2) = .TRUE.
EPT( I,IY,1) = .TRUE.
C
200 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
333 ICD = 0
OMG = .FALSE.
PSI = .FALSE.
SLP = .FALSE.
C
MOV = .FALSE.
CNT = .FALSE.
PRS = .FALSE.
C
READ(1,DATS,END=310)
IF ( OMG ) GO TO 205
IF ( PSI ) GO TO 205
IF ( PRS ) GO TO 205
IF ( CNT ) GO TO 205
C
202 INER = 1
GO TO 333
C
205 IF ( ICD.LE.0 ) GO TO 202
C
C ---------------------------------------------------------
DO 400 IICT = 1, ICD
C
ISLP = - 10
C
C ------------------------------------------------
DO 420 L = 1, 5
C
INDX( L) = 0
INDY( L) = 0
INDE( L) = 0
C
420 CONTINUE
C ------------------------------------------------
C
WVL = 0.D0
PVL = 0.D0
TVL = 0.D0
PRL = 0.D0
C
C -------- Important ---
C
READ(1,1000, END=310, ERR=225 ) ISLP, WVL, PVL,
1 ( INDX( L), INDY( L), INDE( L), L = 1, 5 )
1000 FORMAT( I2,1X,F6.3,1X,F6.3, 5(3I3,1X))
C
IPT = 0
C
C IF ISLP = -1 Boundary SLOPES DOWNW.
C IF ISLP = 0 Boundary IS HOR.
C IF ISLP = 1 Boundary SLOPES UPW.
C IF 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
C
225 INER = 1
GO TO 400
C
230 IPT = IPT + 1
231 IPT = IPT + 1
232 IPT = IPT + 1
233 IPT = IPT + 1
C
C IPT = 1, Downward / = 2, Hor./ = 3, Upward / = 4, Vert.
C
C ------------------------------------------------
DO 440 L = 1, 5
C
IUPR = INDE( L)
GO TO ( 240, 240, 240, 235 ) IPT
235 ILOW = INDY( L)
GO TO 245
240 ILOW = INDX( L)
245 I = INDX( L)
J = INDY( L)
C
C --------------------------------------
DO 460 INDEX = ILOW, IUPR
C
IF ( I) 255, 440, 250
250 IF ( J) 255, 440, 251
251 IF ( I-IX) 252, 252, 255
252 IF ( J-IY) 260, 260, 255
C
C ----- If I or J = 0 : No Boundary specified ---
C
255 INER = 1
WRITE(16,2000) L, IX, IY
2000 FORMAT(' * Error on DATA ',I2,2X,' IX =',I3,' IY =',I3)
GO TO 440
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) = WVL
C
265 IF (.NOT. PSI ) GO TO 270
C
PPT( I,J,1) = .FALSE.
PPT( I,J,2) = SLP
PPT( I,J,3) = MOV
PS1( I,J) = PVL
C
270 CONTINUE
IF (.NOT. CNT) GO TO 278
WPT( I,J,4) = .TRUE.
PPT( I,J,4) = .TRUE.
TPT( I,J,4) = .TRUE.
PRT( I,J,4) = .TRUE.
C
278 IF (.NOT. PRS) GO TO 280
PRT( I,J,1) = .FALSE.
PRT( I,J,2) = SLP
PRT( I,J,3) = MOV
PR1( I,J) = PRL
280 GO TO ( 286, 288, 287, 285 ), IPT
C
285 J = J + 1
GO TO 460
C
286 J = J - 1
GO TO 288
C
287 J = J + 1
288 I = I + 1
460 CONTINUE
C --------------------------------------
C
C ----- From ILOW to INDE ( L) ---
C
440 CONTINUE
C ------------------------------------------------
C
400 CONTINUE
C ----------------------------------------------------------
GO TO 333
C
310 CONTINUE
C
IF ( INER.EQ.0) RETURN
C
RETURN
END
C **********************************************************************
C
SUBROUTINE INIT ( PPT, WPT, TPT, CPT, EPT, PS1, PS3, WP1,
1 CP1, PS, FF, W, U, V, UB, VB,
2 UN, VN, HH, BX, BY, RX, RY, XZ,
3 YZ, NDE, HX, HY, IX, IY )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
INTEGER*2 NDE
C
LOGICAL*1 WPT( IX,IY,4), PPT( IX,IY,4), TPT( IX,IY,4),
1 EPT( IX,IY,4), CPT( IX,IY,4)
C
DIMENSION PS1( IX,IY), PS( IX,IY), W( IX,IY), WP1( IX,IY),
1 FF ( IX,IY)
C
DIMENSION CP1( IX,IY), U ( IX,IY), UB ( IX,IY), V ( IX,IY),
1 VB ( IX,IY), HH( IX,IY), BX ( IX,IY), BY( IX,IY),
2 RX ( IX,IY), RY( IX,IY), PS3( IY), XZ( IX),
3 YZ ( IY), UN( IX,IY), VN ( IX,IY)
C
LOGICAL PRS
COMMON /IDT/ DT, VIS, DX, DY, DX2, DY2, ITT, CON, PRS, RLL
COMMON /DDS/ DS1, DS2, DS3, DSV( 10), DIS( 100)
C
DIMENSION HX( IX-1), HY( IY-1)
C
IXY = IX* IY ! NNP
MSH = ( IX-1)*( IY-1) ! MSH
C
WRITE(6,*) ' '
WRITE(6,2000)
2000 FORMAT(' ','*** RL = ? ')
*** READ (5,*) RLL
C
RLL = 30.D0
RLL = 35.D0
C
C -----------------
DS1 = 0.D0
DS2 = 0.D0
DS3 = 0.D0
C
C ----------------------------
DO 50 I = 1, 10
C
DSV( I)= 0.D0
50 CONTINUE
C ----------------------------
DO 70 I = 1, 100
C
DIS( I)= 0.D0
70 CONTINUE
C ------------------------------------------------
DO 100 I = 1, IX
C
C --------------------------------------
DO 120 J = 1, IY
C
C ----------------------------
DO 140 K = 1, 4
C
PPT( I,J,K) = .FALSE.
WPT( I,J,K) = .FALSE.
TPT( I,J,K) = .FALSE.
C
CPT( I,J,K) = .FALSE.
EPT( I,J,K) = .FALSE.
C
140 CONTINUE
C ----------------------------
C
PS1( I,J) = 0.D0
WP1( I,J) = 0.D0
CP1( I,J) = 0.D0
PS ( I,J) = 0.D0
FF ( I,J) = 0.D0
W ( 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
UN ( I,J) = 0.D0
VN ( I,J) = 0.D0
C
HH ( I,J) = 0.D0
BX ( I,J) = 0.D0
BY ( I,J) = 0.D0
RX ( I,J) = 0.D0
RY ( I,J) = 0.D0
C
XZ ( I) = 0.D0
YZ ( J) = 0.D0
PS3( J) = 0.D0
C
120 CONTINUE
C --------------------------------------
C
100 CONTINUE
C ------------------------------------------------
C
C ***** Attention *** DATA: Mesh with equal length ***
C
DO 200 I = 1, IX
C
XZ ( I) = RLL* DFLOAT(I-1)/ DFLOAT( IX-1)
200 CONTINUE
C ----------------------------------------------------
DO 300 J = 1, IY
C
YZ ( J) = 2.D0* DFLOAT( J-1)/ DFLOAT( IY-1)
300 CONTINUE
C ----------------------------------------------------
DO 400 I = 1, IX-1
C
HX( I) = DABS( XZ( I+1) - XZ( I))
400 CONTINUE
C -----------------------------------------
DO 500 J = 1, IY-1
C
HY( J) = DABS( YZ( J+1) - YZ( J))
500 CONTINUE
C -----------------------------------------
C
CALL MIN ( HX, IX, HXMN )
C
CALL MIN ( HY, IY, HYMN )
C
C -----------------------------------------
C
IF ( HXMN.GE.HYMN) GO TO 900
HXYI = ( HXMN** 3)* 0.1D0
GO TO 450
C
900 HXYI = ( HYMN** 3)* 0.1D0
450 CONTINUE
C
RETURN
END
C **********************************************************************
C
SUBROUTINE MIN ( H, IWK, HMIN)
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION H ( IWK)
C
HMIN = 1.D0
C
DO 100 I = 1, IWK
C
RMIN = HMIN - H( I)
IF ( RMIN.LE.0.D0) GO TO 100
HMIN = H( I)
C
100 CONTINUE
C
RETURN
END
C **********************************************************************
C
SUBROUTINE OUTD ( U, V, UE, VE, PPT, EPT, WP2, IS, JS,
1 IE, XS, YS, PET, EET, XE, YE, PS2, NDE,
2 IX, IY, IXY, IX1, IY1, MSH, KEY )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
INTEGER*2 NDE
C
DIMENSION XS ( IXY), YS ( IXY), XE ( IXY), YE ( IXY),
1 PS2( IXY), WP2( IXY), NDE( MSH,4), U( IX,IY),
2 V ( IX,IY), UE ( IXY), VE ( IXY)
C
LOGICAL*1 PET( IXY,4), EET( IXY,4), PPT( IX,IY,4),
1 EPT( IX,IY,4)
C
COMMON /CTL/ LOP, NPR, ILP, NNN
COMMON /DDS/ DS1, DS2, DS3, DSV( 10), DIS( 100)
LOGICAL PRS
COMMON /IDT/ DT, VIS, DX, DY, DX2, DY2, ITT, CON, PRS, RLL
C
CALL STMP ( U, UE, IX, IY, IXY )
CALL STMP ( V, VE, IX, IY, IXY )
C
CALL STMS ( PPT, PET, IX, IY, IXY )
CALL STMS ( EPT, EET, IX, IY, IXY )
C
WRITE(7,2000) DT* LOP
WRITE(7,2100) LOP, DT
WRITE(7,2200) 1.D0/ VIS
2000 FORMAT(2X,' TIME =', 1X,F10.4)
2100 FORMAT(1X,' LOP =', I6,2X,'DT=',F4.3)
2200 FORMAT(2X,' Re =',F7.2)
WRITE(7,1900) DS3
WRITE(7,1910) DSV( 1)
WRITE(7,1920) DSV( 2)
WRITE(7,1930) DS2
WRITE(7,1940) DS1
1900 FORMAT(' X1R = ',F8.4)
1910 FORMAT(' X2R = ',F8.4)
1920 FORMAT(' X3R = ',F8.4)
1930 FORMAT(' X4R = ',F8.4)
1940 FORMAT(' X5R = ',F8.4)
C
WRITE(8,2300) NNN, LOP
2300 FORMAT(2I6)
WRITE(6,2400)
2400 FORMAT(/,' *** To be continued ? ***')
WRITE(6,*) ' 1: Yes 2: No'
C
** READ(5,*) KEY
KEY = 2
C
IS = 1
JS = IDNINT( DFLOAT( 1+IY)/ 2.D0 )
IE = 1
C
MSH = ( IX-1)* ( IY-1)
C
WRITE(8,2500) IX1, IY1, RLL, IX, IS, JS, IE, IXY, MSH
2500 FORMAT(2I4,F7.2,4I4,2I6)
C
WRITE(8,2600) ( XS( I), I = 1, IXY)
WRITE(8,2600) ( YS( I), I = 1, IXY)
C
WRITE(8,2700) ( ( PET( I,J), J = 1,4), I = 1, IXY)
WRITE(8,2700) ( ( EET( I,J), J = 1,4), I = 1, IXY)
2700 FORMAT(80L1)
WRITE(8,2600) ( XE( I), I = 1, IXY)
WRITE(8,2600) ( YE( I), I = 1, IXY)
C
WRITE(8,2800) ( ( NDE( I,J), J = 1,4), I = 1, MSH )
2800 FORMAT(20I6)
C
WRITE(8,2600) ( PS2( I), I = 1, IXY)
WRITE(8,2600) ( WP2( I), I = 1, IXY)
2600 FORMAT(10F8.4)
C
RETURN
END
C **********************************************************************
C
SUBROUTINE PREP ( XE, YE, XS, YS, XZ, YZ, IMP, JMP, IS,
1 JS, IE, IX, IY, IXY )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION XE( IXY), YE( IXY), XS( IXY), YS( IXY), XZ( IX),
1 YZ( IY)
C
COMMON /CTL/ LOP, NPR, ILP, NNN
COMMON /DDS/ DS1, DS2, DS3, DSV( 10), DIS( 100)
LOGICAL PRS
COMMON /IDT/ DT, VIS, DX, DY, DX2, DY2, ITT, CON, PRS, RLL
C
IMP = DNINT( DFLOAT( IX+1)/ 2.D0 )
JMP = DNINT( DFLOAT( IY+1)/ 2.D0 )
C
IS = 1
JS = DNINT( DFLOAT( 1+IY)/ 2.D0 )
IE = 1
C
C --------------------------
NN = 0
C
DO 100 J = 1, IY
DO 100 I = 1, IX
C
NN = NN + 1
XE( NN) = XZ( I)
YE( NN) = YZ( J)
C
100 CONTINUE
C --------------------------
NN = 0
DO 200 J = 1, IY
DO 200 I = 1, IX
C
NN = NN + 1
IF( I.EQ.IX ) THEN
GO TO 250
END IF
250 CONTINUE
IF ( J.EQ.IY) THEN
GO TO 200
END IF
C
200 CONTINUE
C -------------------------------
NN = 0
DO 300 J = 1, IY
DO 300 I = 1, IX
C
NN = NN + 1
XS( NN) = DX* DFLOAT( I-1)
YS( NN) = DY* DFLOAT( J-1)
C
300 CONTINUE
C -------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE PRGS ( PS1, PS3, PPT, FLW, XX, YY, NDE, XZ,
1 YZ, XE, YE, XS, YS, IMP, JMP, IS,
2 JS, IE, IX, IY, IXY, MSH )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION PS1( IX,IY), PS3( IY), FLW( IY), NDE( MSH,4),
1 XZ ( IX), YZ ( IY), XE ( IXY), YE ( IXY),
2 XS ( IXY), YS ( IXY)
C
LOGICAL*1 PPT( IX,IY,4)
INTEGER*2 NDE
LOGICAL PRS
COMMON /IDT/ DT, VIS, DX, DY, DX2, DY2, ITT, CON, PRS, RLL
C
C ----------------------------------------
PS3( 1) = - 1.D0
N = IY - 1
DO 100 J = 2, N
C
FLW( J) = YZ ( J) - YZ( J-1)
PS3( J) = PS3( J-1) + FLW( J)
C
100 CONTINUE
C ----------------------------------------
C
C ----- Stream function variance point ---
C
DO 150 J = 2, N
C
PS1( IX, J) = PS3( J)
C
150 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)
C
200 CONTINUE
C ----------------------------------------
DO 300 J = 1, IY-1
DO 300 I = 1, IX-1
C
NM = ( J-1)* ( IX-1) + I
C
NDE( NM,1) = ( J-1)* IX + I
NDE( NM,2) = ( J-1)* IX + I + 1
NDE( NM,3) = ( J-1)* IX + I + IX + 1
NDE( NM,4) = ( J-1)* IX + I + IX
C
300 CONTINUE
C ----------------------------------------
C
CALL PREP ( XE, YE, XS, YS, XZ, YZ, IMP, JMP, IS, JS, IE, IX, IY,
1 IXY )
C
C ----------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE PRPS ( PS1, PS2, XE, YE, IE, JS, WP1, WP2, UB,
1 VB, U, V, FF, IX, IY, IXY )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
COMMON /WEH/ WEI, IHS, LTC
COMMON /CTL/ LOP, NPR, ILP, NNN
COMMON /DDS/ DS1, DS2, DS3, DSV( 10), DIS( 100)
C
DIMENSION PS1( IX,IY), PS2( IXY), XE( IXY), YE( IXY),
1 WP1( IX,IY), WP2( IXY), UB( IX,IY), VB( IX,IY),
2 U ( IX,IY), V ( IX,IY), FF( IX,IY)
C
CALL DIST ( PS1, XE, YE, IE, JS, IX, IY )
C
IF ( LOP/NPR* NPR .EQ. LOP .OR. LOP.EQ.1) THEN
C ====================================================
C
CALL STMP ( PS1, PS2, IX, IY, IXY )
C
CALL STMP ( WP1, WP2, IX, IY, IXY )
C
END IF
C
C --------------------------
DO 100 I = 1, IX
DO 100 J = 1, IY
C
UB( I,J) = U( I,J)
VB( I,J) = V( I,J)
C
100 CONTINUE
C --------------------------
C
LTC = 3
C --------------------------
DO 200 J = 1, IY
DO 200 I = 1, IX
C
FF( I,J) = WP1( I,J)
U ( I,J) = 0.D0
V ( I,J) = 0.D0
C
200 CONTINUE
C --------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE SFCF ( PA1, PA2, PA3, PB1, PB2, PB3, HX, HY, IX, IY )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION PA1( IX,IY), PA2( IX,IY), PA3( IX,IY), PB1( IX,IY),
1 PB2( IX,IY), PB3( IX,IY)
C
DIMENSION HX( IX-1), HY( IY-1)
C
C -------------------------
DO 100 J = 1, IY
DO 100 I = 1, IX
C -------------------------
C
PA1( I,J) = 0.D0
PA2( I,J) = 0.D0
PA3( I,J) = 0.D0
C
PB1( I,J) = 0.D0
PB2( I,J) = 0.D0
PB3( I,J) = 0.D0
C
100 CONTINUE
C -------------------------
DO 200 J = 2, IY-1
DO 200 I = 2, IX-1
C -------------------------
C
PA1( I,J) = 2.D0/( HX(I-1)* ( HX(I-1) + HX( I)))
PA2( I,J) = 2.D0/( HX(I-1)* HX( I))
PA3( I,J) = 2.D0/( HX( I) * ( HX(I-1) + HX( I)))
C
PB1( I,J) = 2.D0/( HY( J-1)* ( HY( J-1) + HY( J)))
PB2( I,J) = 2.D0/( HY( J-1)* HY( J))
PB3( I,J) = 2.D0/( HY( J)* ( HY( J-1) + HY( J)))
C
200 CONTINUE
C -------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE SFCL ( PS1, PS, W, WP1, PPT, EPT, ACP, ICT,
1 YZ, PA1, PA2, PA3, PB1, PB2, PB3, HX,
2 HY, IX, IY, IXY )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
COMMON /CTL/ LOP, NPR, ILP, NNN
LOGICAL PRS
COMMON /IDT/ DT, VIS, DX, DY, DX2, DY2, ITT, CON, PRS, RLL
C
DIMENSION PA1( IX,IY), PA2( IX,IY), PA3( IX,IY), PB1( IX,IY),
1 PB2( IX,IY), PB3( IX,IY)
C
DIMENSION PS1( IX,IY), PS( IX,IY), W ( IX,IY), WP1( IX,IY),
1 YZ ( IY), HX( IX-1), HY( IY-1)
C
LOGICAL*1 PPT( IX,IY,4), EPT( IX,IY,4), DRCT
C
RCF = DY2* DY2/( 2.D0* DX2* ( DX2 + DY2 ))
RS2 = DX2/ DY2
C
IUU = 0
IUD = 0
C
ITT = 0
EPS = 0.00005D0
C
IF ( ICT.GT.1) GO TO 30
ICT = ICT + 1
GO TO 40
30 CONTINUE
C
C ----- Attention --------------------------
DO 50 J = 1, IY
DO 50 I = 1, IX
C
PS1( I,J) = 2.D0* PS( I,J) - W( I,J)
C
50 CONTINUE
C ------------------------------------------
C
40 CONTINUE
DRCT = .FALSE.
JDM = IY + 1
105 IDFC = 0
ITT = ITT + 1
C
IF ( ITT.GT.IXY/10) GO TO 777
C
C ----------------------------------------------------------
DO 900 JJ = 1, IY
C
99 DRCT = .NOT. DRCT
C
IF ( DRCT ) J = JJ
IF ( .NOT.DRCT ) J = JDM - JJ
C
C ------------------------------------------------
DO 100 I = 1, IX
C
IF ( I.LE.1) GO TO 120
IF ( J.LE.1) GO TO 120
C
DXC = HX(I-1)* HX( I)* ( HX(I-1) + HX( I))
DYC = HY( J-1)* HY( J)* ( HY( J-1) + HY( J))
C
DXY = ( HY( J-1) + HY( J))* DXC + ( HX( I-1) + HX( I))* DYC
C
120 IF ( PPT( I,J,1)) GO TO 1
IF ( PPT( I,J,2)) GO TO 100
GO TO 100
C
1 CONTINUE
C
IF ( I.EQ.IX ) GO TO 130
IF ( I.EQ.1 ) GO TO 140
C
IF ( EPT( I,J,3)) GO TO 130
IF ( EPT( I,J,4)) GO TO 140
IF ( EPT( I,J,1)) GO TO 150
IF ( EPT( I,J,2)) GO TO 160
C
CONTINUE
C
RDPS = ( PA1( I,J)* PS1( I-1,J)
1 + PA3( I,J)* PS1( I+1,J) + PB1( I,J)* PS1( I,J-1)
2 + PB3( I,J)* PS1( I,J+1) + WP1( I,J))
3 /( PA2( I,J) + PB2( I,J))
C
RAB2 = DABS( RDPS - PS1( I,J))
PS1( I,J) = PS1( I,J) + ACP* ( RDPS - PS1( I,J))
C
IF ( RAB2.GT.EPS) GO TO 170
IF ( ITT .LT.10 ) GO TO 170
GO TO 100
C
170 IDFC = 1
IKX = I
IKY = J
GO TO 100
C
130 IF ( PPT( I,J,4)) GO TO 180
C
C ----- Periodic right ---
C
PS1( I,J) = PS1( I-(IX-1),J)
GO TO 100
C
C ----- Continuative right ---
C
180 PS1( I,J) = ( HY( J)* PS1( I,J-1)+ HY( J-1)* PS1( I,J+1)
1 + 0.5D0* DYC* WP1( I,J))/( HY( J-1) + HY( J))
GO TO 100
C
140 IF ( PPT( I,J,4)) GO TO 190
C
C ----- Periodic left ---
C
DELP = ( PS1( I+1,J) + PS1( I-1 + ( IX-1), J)
C
1 + ( DX2 / 2.D0* ( DX2 + DY2 ))
2 * ( PS1( I,J+1) + PS1( I,J-1) + DY2* WP1( I,J)))
C
RDLP =( ( HX(I-1)*PS1(I+1,J) + HX( I)* PS1(I-1+(IX-1),J))* DYC
C
1 + ( HY( J-1)* PS1( I,J+1) + HY( J)* PS1( I,J-1))* DXC
2 + 0.5D0* DXC* DYC* WP1( I,J))/ DXY
C
PS1( I,J) = PS1( I,J) + ACP* ( RDLP - PS1( I,J))
GO TO 100
C
C ----- Continuative Left ---
C
190 CONTINUE
C
YTR = ( YZ( J) - 1.D0 )/ 2.D0
PS1( I,J) = - 1.D0 + 8.D0* YTR* YTR* ( 3.D0 - 4.D0* YTR )
GO TO 100
C
150 IF ( PPT( I,J,4)) GO TO 192
C
PS1( I,J) = PS1( I,J-(IY-1))
GO TO 100
C
192 PS1( I,J) = ( HX( I)* PS1( I-1,J) + HX(I-1)* PS1( I+1,J)
1 + 0.5D0* DXC* WP1( I,J))/( HX(I-1) + HX( I))
GO TO 100
C
160 IF ( PPT( I,J,4)) GO TO 194
DELP = ( PS1( I,J+1) + PS1( I,J-1+( IY-1)) + ( PS1( I+1,J)
1 + PS1( I-1,J) + DX2* WP1( I,J))/ RS2 )* RCF
PS1( I,J) = PS1( I,J) + ACP* ( DELP - PS1( I,J))
GO TO 100
C
194 PS1( I,J) = ( HX( I)* PS1( I-1,J) + HX(I-1)* PS1( I+1,J)
1 + 0.5D0* DXC* WP1( I,J))/( HX(I-1) + HX( I))
GO TO 100
C
100 CONTINUE
C ------------------------------------------------
C
IF ( DRCT) GO TO 99
900 CONTINUE
C ----------------------------------------------------------
C
IF ( IUD.EQ.0) GO TO 101
PSS = 0.D0
C
C ------------------------------------------------
DO 200 I = 1, IX
C
PSS = PSS + PS1( I,2) + 0.5D0* DY2* WP1( I,1)
C
200 CONTINUE
C ------------------------------------------------
PSS = ( PSS/ IX)
C
C ------------------------------------------------
DO 300 I = 1, IX
C
PS1( I,1) = PSS
C
300 CONTINUE
C ------------------------------------------------
C
101 IF ( IUU.EQ.0) GO TO 102
PSS = 0.D0
C
C ------------------------------------------------
DO 400 I = 1, IX
C
PSS = PSS + PS1( I, IY-1) + 0.5D0* DY2* WP1( I,IY)
C
400 CONTINUE
C ------------------------------------------------
DO 500 I = 1, IX
C
PS1( I,IY) = PSS
C
500 CONTINUE
C ------------------------------------------------
102 IF ( IDFC.GT.0 ) GO TO 105
C
RETURN
C
777 CONTINUE
WRITE(6,2000) IKX, IKY, ITT, LOP
2000 FORMAT(' ** ITT.GT.IXY/10 ',I3,1X,I3,' * Psi: ???',' *ITT=',I5,
1 ' at LOP =', I6,' * STOP *')
STOP
END
C **********************************************************************
C
SUBROUTINE STMP ( DIM1, DIM2, IX, IY, IXY )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION DIM1( IX,IY), DIM2 ( IXY)
C
NN = 0
C
C ---------------------------------
DO 100 J = 1, IY
DO 100 I = 1, IX
C
NN = NN + 1
DIM2( NN) = DIM1( I,J)
C
100 CONTINUE
C ---------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE STMS ( DIM1, DIM2, IX, IY, IXY )
C
LOGICAL*1 DIM1( IX,IY,4), DIM2( IXY,4)
C
NN = 0
C
C --------------------------------------
DO 100 J = 1, IY
DO 100 I = 1, IX
C --------------------------------------
C
NN = NN + 1
C ----------------------------
DO 150 K = 1, 4
C
DIM2( NN,K) = DIM1( I,J,K)
C
150 CONTINUE
C ----------------------------
C
100 CONTINUE
C --------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE TTLE
C
WRITE(6,*) ' '
WRITE(6,2100) '************************************************* '
WRITE(6,2100) ' '
WRITE(6,2100) ' F 2 D - B S G '
WRITE(6,2100) ' '
WRITE(6,2100)' Two Dimensional Backward-Facing Step '
WRITE(6,2100)' '
WRITE(6,2100)' Fluid Flow Analysis '
WRITE(6,2100)' '
WRITE(6,2100)' by the Fourth-ordered Finite Difference Method '
WRITE(6,2100)' '
WRITE(6,2100)' ( Version : 3.8 ) '
WRITE(6,2100)' '
WRITE(6,2100)' Copyright 2011 : Yasuhiro MATSUDA '
WRITE(6,2100)' '
WRITE(6,2100)'************************************************** '
2100 FORMAT(A55)
C
RETURN
END
C **********************************************************************
C
SUBROUTINE WFCA ( WP1, PS1, WPT, RAS, RSS, IX, IY )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION WP1( IX,IY), PS1( IX,IY), RSS( IX-1), RAS( IY-1)
LOGICAL*1 WPT( IX,IY,4)
C
C ----- No SLP surfaces : Updated ---
C
C -----------------------------------------------
DO 100 J = 1, IY
DO 100 I = 1, IX
C -----------------------------------------------
C
IF ( WPT( I,J,1)) GO TO 100
IF ( WPT( I,J,2)) GO TO 50
IF ( WPT( I,J,3)) GO TO 200
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))/ RSS( I)
C
13 IF ( .NOT. WPT( I-1,J,1)) GO TO 3
WB = ( PS1( I,J) - PS1( I-1,J))/ RSS(I-1)
C
12 IF ( .NOT. WPT( I,J+1, 1)) GO TO 2
WC = ( PS1( I,J) - PS1( I,J+1))/ RAS( J)
C
14 IF ( .NOT. WPT( I,J-1, 1)) GO TO 4
WD = ( PS1( I,J) - PS1( I,J-1))/ RAS( J-1)
C
114 WP1( I,J) = WA + WB + WC + WD
GO TO 100
C
1 WA = ( PS1( I,J) - PS1( I-1,J))/ RSS(I-1)
GO TO 13
C
3 WB = ( PS1( I,J) - PS1( I+1,J))/ RSS( I)
GO TO 12
C
2 WC = ( PS1( I,J) - PS1( I,J-1))/ RAS( J-1)
GO TO 14
C
4 WD = ( PS1( I,J) - PS1( I,J+1))/ RAS( J)
GO TO 114
C
C ----- No-Slip Surfaces at Edges ---
C
10 WP1( I,J) = 3.D0* ( PS1( I,J) - PS1( I,J-1))
1 / RAS( J-1) - 0.5D0* WP1( I,J-1)
C
C ----- Unequaled Length ---
C
GO TO 100
C
20 WP1( I,J) = 3.D0* ( PS1( I,J) - PS1( I,J+1))
1 / RAS( J) - 0.5D0* WP1( I,J+1)
C
C ----- Unequaled Length ---
C
GO TO 100
C
30 WA = ( PS1( I,J) - PS1( I-1,J))/ RSS(I-1)
WB = WA
GO TO 12
C
40 WA = ( PS1( I,J) - PS1( I+1,J))/ RSS( I)
WB = WA
GO TO 12
C
50 WP1( I,J) = 0.D0
GO TO 100
C
200 IF ( J.EQ.IY) GO TO 10
IF ( J.EQ.1 ) GO TO 20
C
100 CONTINUE
C ------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE WUCF ( UNA, UNB, UNC, UND, UNE, UAP, UBP, UCP,
1 UDP, UEP, UAM, UBM, UCM, UDM, UEM, UF1,
2 UF2, UF3, HX, HY, IX, IY )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION AA( 4), BB( 4), CC( 4), DD( 4), EE( 4)
C
DIMENSION UNA( IX,4), UNB( IX,4), UNC( IX,4), UND( IX,4),
1 UNE( IX,4)
C
DIMENSION UAP( IX,4), UBP( IX,4), UCP( IX,4), UDP( IX,4),
1 UEP( IX,4)
C
DIMENSION UAM( IX,4), UBM( IX,4), UCM( IX,4), UDM( IX,4),
1 UEM( IX,4)
C
DIMENSION UF1( IY), UF2( IY), UF3( IY), HX( IX-1), HY( IY-1)
C
C ------------------------------
DO 100 K = 1, 4
C
AA( K) = 0.D0
BB( K) = 0.D0
CC( K) = 0.D0
DD( K) = 0.D0
EE( K) = 0.D0
C
100 CONTINUE
C
BB( 1) = 1.D0
CC( 2) = 2.D0
DD( 3) = 6.D0
EE( 4) = 24.D0
C
C ------------------------------
DO 200 I = 1, IX
DO 200 K = 1, 4
C ------------------------------
C
UNA( I,K) = 0.D0
UNB( I,K) = 0.D0
UNC( I,K) = 0.D0
UND( I,K) = 0.D0
UNE( I,K) = 0.D0
C
UAP( I,K) = 0.D0
UBP( I,K) = 0.D0
UCP( I,K) = 0.D0
UDP( I,K) = 0.D0
UEP( I,K) = 0.D0
C
UAM( I,K) = 0.D0
UBM( I,K) = 0.D0
UCM( I,K) = 0.D0
UDM( I,K) = 0.D0
UEM( I,K) = 0.D0
C
200 CONTINUE
C -----------------------------------
DO 300 J = 1, IY
C
UF1( J) = 0.D0
UF2( J) = 0.D0
UF3( J) = 0.D0
C
300 CONTINUE
C -----------------------------------
DO 400 J = 2, IY-1
C
UF1( J) = - HY( J)/( HY( J-1)* ( HY( J-1) + HY( J)))
UF2( J) = ( HY( J) - HY( J-1))/( HY( J-1)* HY( J))
UF3( J) = HY( J-1)/( HY( J)* ( HY( J-1) + HY( J)))
C
400 CONTINUE
C ------------------------------------------------
DO 500 I = 3, IX-2
C
HA = HX( I-2)* ( HX( I-2) + HX(I-1))
1 * ( HX( I-2) + HX(I-1) + HX( I))
2 * ( HX( I-2) + HX(I-1) + HX( I) + HX( I+1))
C
HB = HX(I-1)*( HX(I-1) + HX( I))*( HX(I-1) + HX( I) + HX( I+1))
C
HD = HX( I)* HX( I+1)
HE = HX( I) + HX( I+1)
C
C --------------------------------------
DO 600 K = 1, 4
C
UNA( I,K) = ( BB( K)* HX(I-1)* HX( I)* ( HX( I) + HX( I+1))
1 - CC( K)* ( HX(I-1)* ( 2.D0 * HX( I) + HX( I+1))
2 - HX( I)* ( HX( I) + HX( I+1)))
3 + DD( K)* ( HX(I-1) - 2.D0 * HX( I)
4 - HX( I+1)) + EE( K))/ HA
C
UNB( I,K) = ( - UNA( I,K)* ( HX( I-2) + HX(I-1))
1 * ( HX( I-2) + HX(I-1) + HX( I))
2 * ( HX( I-2) + HX(I-1) + HX( I) + HX( I+1))
3 - BB( K)* HX( I)* ( HX( I) + HX( I+1))
4 + CC( K)* ( 2.D0 * HX( I) + HX( I+1)) - DD( K))/ HB
C
UND( I,K) = ( UNA( I,K)* ( HX( I-2) + HX(I-1))
1 * ( HX( I-2) + HX(I-1) + HX( I) + HX( I+1))
2 + UNB( I,K)* HX(I-1)* ( HX(I-1) + HX( I) + HX( I+1))
3 + BB( K)* ( HX( I) + HX( I+1)) - CC( K))/ HD
C
UNE( I,K) = ( UNA( I,K)* ( HX( I-2) + HX(I-1))
1 + UNB( I,K)* HX(I-1) - UND( I,K)* HX( I)
2 + BB( K))/ HE
C
UNC( I,K) = AA( K) - ( UNA(I,K) + UNB(I,K) + UND(I,K) + UNE(I,K))
C
600 CONTINUE
C --------------------------------------
C
500 CONTINUE
C ------------------------------------------------
DO 800 I = 4, IX-1
C
HA = HX( I-3)* ( HX( I-3) + HX( I-2))
1 * ( HX( I-3) + HX( I-2) + HX(I-1))
2 * ( HX( I-3) + HX( I-2) + HX(I-1) + HX( I))
C
HB = HX(I-2)*( HX(I-2) + HX(I-1))*( HX(I-2) + HX(I-1) + HX(I))
C
HC = HX(I-1)* ( HX(I-1) + HX( I))
C
HE = HX( I)
C
C --------------------------------------
DO 900 K = 1, 4
C
UAP( I,K) = ( - BB( K)* HX(I-1)* HX( I)* ( HX( I-2) + HX(I-1))
1 + CC( K)* ( HX( I-2)* HX(I-1)
2 + ( HX(I-1)**2) - HX( I-2)* HX( I)
3 - 2.D0 * HX(I-1)* HX( I))
4 + DD( K)* ( HX( I-2) + 2.D0 * HX(I-1) - HX( I))
5 + EE( K))/ HA
C
UBP( I,K) = ( - UAP( I,K)* ( HX( I-3) + HX( I-2))* ( HX( I-3)
1 + HX( I-2) + HX(I-1))
2 * ( HX( I-3) + HX( I-2) + HX(I-1) + HX( I))
3 + BB( K)* HX(I-1)* HX( I) - CC( K)* ( HX(I-1) - HX( I))
4 - DD( K))/ HB
C
UCP( I,K) = ( - UAP( I,K)* ( HX( I-3) + HX( I-2) + HX(I-1))
1 * ( HX( I-3) + HX( I-2) + HX(I-1) + HX( I))
2 - UBP( I,K)* ( HX( I-2) + HX(I-1))
3 * ( HX( I-2) + HX(I-1) + HX( I))
4 - BB( K)* HX( I) + CC( K))/ HC
C
UEP( I,K) = ( UAP( I,K)* ( HX( I-3) + HX( I-2) + HX(I-1))
1 + UBP( I,K)* ( HX( I-2) + HX(I-1))
2 + UCP( I,K)* HX(I-1) + BB( K))/ HE
C
UDP( I,K) = AA(K) - ( UAP(I,K) + UBP(I,K) + UCP(I,K) + UEP(I,K))
C
900 CONTINUE
C --------------------------------------
C
800 CONTINUE
C ------------------------------------------------
DO 950 I = 2, IX-3
C
HE = HX( I+2)* ( HX( I+1) + HX( I+2))
1 * ( HX( I) + HX( I+1) + HX( I+2))
2 * ( HX(I-1) + HX( I) + HX( I+1) + HX( I+2))
C
HD = HX( I+1)* ( HX(I) + HX(I+1))* ( HX(I-1) + HX(I) + HX(I+1))
C
HC = HX( I)* ( HX(I-1) + HX( I))
C
HA = HX(I-1)
C
C --------------------------------------
DO 970 K = 1, 4
C
UEM( I,K) = ( BB( K)* HX(I-1)* HX( I)
1 * ( HX( I) + HX( I+1)) + CC( K)* ( ( HX( I)**2)
2 + HX( I)* HX( I+1) - HX(I-1)* HX( I+1)
3 - 2.D0 * HX(I-1)* HX( I)) + DD( K)* ( HX(I-1)
4 - 2.D0 * HX( I) - HX( I+1)) + EE( K))/ HE
C
UDM( I,K) = ( - UEM( I,K)* ( HX( I+1) + HX( I+2))
1 * ( HX( I) + HX( I+1) + HX( I+2))
2 * ( HX(I-1) + HX( I) + HX( I+1)
3 + HX( I+2)) - BB( K)* HX(I-1)* HX( I)
4 - CC( K)* ( HX( I) - HX(I-1)) + DD( K))/ HD
C
UCM( I,K) = ( - UEM( I,K)* ( HX( I)+ HX( I+1) + HX( I+2))
1 * ( HX(I-1) + HX( I) + HX( I+1) + HX( I+2))
2 - UDM( I,K)* ( HX( I) + HX( I+1))
3 * ( HX(I-1) + HX( I) + HX( I+1))
4 + BB( K)* HX(I-1) + CC( K))/ HC
C
UAM( I,K) = ( UEM( I,K)* ( HX( I) + HX( I+1) + HX( I+2))
1 + UDM( I,K)* ( HX( I) + HX( I+1))
2 + UCM( I,K)* HX( I) - BB( K))/ HA
C
UBM( I,K) = AA( K) - ( UAM(I,K) + UCM(I,K) + UDM(I,K) + UEM(I,K))
C
970 CONTINUE
C --------------------------------------
C
950 CONTINUE
C ------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE WUCL ( WP1, WM1, PS1, WPT, EPT, U, YZ, IX, IY,
1 UNA, UNB, UNC, UND, UNE, UAP, UBP, UCP, UDP,
2 UEP, UAM, UBM, UCM, UDM, UEM, UF1, UF2, UF3,
3 HX, HY )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION WP1( IX,IY), WM1( IX,IY), PS1( IX,IY), U( IX,IY),
1 YZ ( IY), AA ( 4), BB ( 4), CC ( 4),
2 DD ( 4), EE ( 4), UPA( 4), UPB( 4),
3 UPC( 4), UPD( 4)
C
DIMENSION UNA( IX,4), UNB( IX,4), UNC( IX,4), UND( IX,4),
1 UNE( IX,4)
C
DIMENSION UAP( IX,4), UBP( IX,4), UCP( IX,4), UDP( IX,4),
1 UEP( IX,4)
C
DIMENSION UAM( IX,4), UBM( IX,4), UCM( IX,4), UDM( IX,4),
1 UEM( IX,4)
C
DIMENSION UF1( IY), UF2( IY), UF3( IY), HX( IX-1), HY( IY-1)
C
LOGICAL*1 WPT( IX,IY,4), EPT( IX,IY,4)
C
LOGICAL PRS
COMMON /IDT/ DT, VIS, DX, DY, DX2, DY2, ITT, CON, PRS, RLL
COMMON /WEH/ WEI, IHS, LTC
C
C ---------------------------------------
DO 100 K = 1, 4
C
AA( K) = 0.D0
BB( K) = 0.D0
CC( K) = 0.D0
DD( K) = 0.D0
EE( K) = 0.D0
C
100 CONTINUE
C
BB( 1) = 1.D0
CC( 2) = 2.D0
DD( 3) = 6.D0
EE( 4) = 24.D0
C
C ----------------------------------------------------------
DO 200 J = 1, IY
DO 200 I = 1, IX
C ----------------------------------------------------------
C
IF ( WPT( I,J,1)) THEN
IF (.NOT.EPT( I,J,3)) GO TO 10
U( I,J) = UF1( J)* PS1( I,J-1) + UF2( J)* PS1( I,J)
1 + UF3( J)* PS1( I,J+1)
GO TO 200
10 IF ( .NOT.EPT( I,J,4)) GO TO 30
YTR = ( YZ( J) - 1.D0 )/ 2.D0
U( I,J) = 24.D0* YTR* ( 1.D0 - 2.D0* YTR )
GO TO 200
C
30 CONTINUE
U( I,J) = UF1( J)* PS1( I,J-1) + UF2( J)* PS1( I,J)
1 + UF3( J)* PS1( I,J+1)
GO TO 115
C
ELSE
IF ( .NOT.EPT( I,J,1)) GO TO 40 ! UE
U( I,J) = 0.D0
GO TO 200
40 IF ( .NOT.EPT( I,J,2)) GO TO 50 ! SITA
U( I,J) = 0.D0
GO TO 200
C
50 IF ( .NOT.EPT( I,J,4)) GO TO 61 ! STEP
U( I,J) = 0.D0
GO TO 200
C
61 CONTINUE
END IF
C
115 CONTINUE
UT = U( I,J)* DT
C
IF ( EPT( I,J,4)) GO TO 80
IF ( EPT( I,J,3)) GO TO 80
IF ( EPT( I-1,J,4)) GO TO 29
IF ( EPT( I+1,J,3)) GO TO 29
C
CENT = - UT* ( UNA( I,1)* WM1( I-2,J) + UNB( I,1)* WM1( I-1,J)
1 + UNC( I,1)* WM1( I,J) + UND( I,1)* WM1( I+1,J)
2 + UNE( I,1)* WM1( I+2,J))
C
3 + ( UT** 2)* ( UNA( I,2)* WM1( I-2,J) + UNB( I, 2)* WM1( I-1,J)
4 + UNC( I,2)* WM1( I,J) + UND( I, 2)* WM1( I+1,J)
5 + UNE( I,2)* WM1( I+2,J))/ 2.D0
C
6 - ( UT** 3)* ( UNA( I,3)* WM1( I-2,J) + UNB( I,3)* WM1( I-1,J)
7 + UNC( I,3)* WM1( I,J) + UND( I,3)* WM1( I+1,J)
8 + UNE( I,3)* WM1( I+2,J))/ 6.D0
C
9 + ( UT** 4)* ( UNA( I,4)* WM1( I-2,J) + UNB( I,4)* WM1( I-1,J)
A + UNC( I,4)* WM1( I,J) + UND( I,4)* WM1( I+1,J)
B + UNE( I,4)* WM1( I+2,J))/ 24.D0
C
GO TO 80
C
29 CONTINUE
A1 = HX( I)/( HX(I-1)* ( HX(I-1) + HX( I)))
B1 = ( HX(I-1) - HX( I))/( HX(I-1)* HX( I))
C1 = HX(I-1)/( HX( I)* ( HX(I-1) + HX( I)))
C
A2 = 1.D0/( HX(I-1)* ( HX(I-1) + HX( I)))
B2 = 1.D0/( HX(I-1)* HX( I))
C2 = 1.D0/( HX( I)* ( HX(I-1) + HX( I)))
C
CENT = - UT* ( - A1* WM1(I-1,J) - B1* WM1(I,J) + C1* WM1(I+1,J))
1 + ( UT** 2)*( A2* WM1(I-1,J) - B2* WM1(I,J) + C2* WM1(I+1,J))
C
80 CONTINUE
IF ( U( I,J).LT.0) GO TO 90
IF ( EPT( I,J,4)) GO TO 900
IF ( EPT( I,J,3)) GO TO 130
C
IF ( EPT( I-1,J,4)) GO TO 110
IF ( EPT( I-2,J,4)) GO TO 120
C
UPW = - UT* ( UAP( I,1)* WM1( I-3,J) + UBP( I,1)* WM1( I-2,J)
1 + UCP( I,1)* WM1( I-1,J) + UDP( I,1)* WM1( I,J)
2 + UEP( I,1)* WM1( I+1,J))
C
3 + ( UT** 2 )* ( UAP( I,2)* WM1( I-3,J) + UBP( I,2)* WM1( I-2,J)
4 + UCP( I,2)* WM1( I-1,J) + UDP( I,2)* WM1( I,J)
5 + UEP( I,2)* WM1( I+1,J))/ 2.D0
C
6 - ( UT** 3 )* ( UAP( I,3)* WM1( I-3,J) + UBP( I,3)* WM1( I-2,J)
7 + UCP( I,3)* WM1( I-1,J) + UDP( I,3)* WM1( I,J)
8 + UEP( I,3)* WM1( I+1,J))/ 6.D0
C
9 + ( UT** 4 )* ( UAP( I,4)* WM1( I-3,J) + UBP( I,4)* WM1( I-2,J)
A + UCP( I,4)* WM1( I-1,J) + UDP( I,4)* WM1( I,J)
B + UEP( I,4)* WM1( I+1,J))/ 24.D0
GO TO 150
C
900 WP1( I,J) = WP1( I,J) - UT* WM1( I,J)/ HX( I)
GO TO 200
C
110 UPW = - UT* ( - WM1( I-1,J) + WM1( I,J))/ HX(I-1)
GO TO 150
C
120 CONTINUE
C
C ------------------------------------------------
DO 160 K = 1, 3
C
HA = HX( I-2)* ( HX(I-2) + HX(I-1))*( HX(I-2) + HX(I-1) + HX(I))
C
HB = HX(I-1)* ( HX(I-1) + HX( I))
HD = HX( I)
C
UPA( K) = ( BB( K)* HX(I-1)* HX( I)
1 + CC( K)* ( HX( I) - HX(I-1)) - DD( K))/ HA
C
UPB( K) = ( - UPA( K)* ( HX( I-2) + HX(I-1))
1 * ( HX( I-2) + HX(I-1) + HX( I))
2 - BB( K)* HX( I) + CC( K))/ HB
C
UPD( K) = ( UPA( K)* ( HX( I-2) + HX(I-1))
1 + UPB( K)* HX(I-1) + BB( K))/ HD
C
UPC( K) = AA( K) - ( UPA( K) + UPB( K) + UPD( K))
C
160 CONTINUE
C ------------------------------------------------
C
UPW = - UT* ( UPA( 1)* WM1( I-2,J) + UPB( 1)* WM1( I-1,J)
1 + UPC( 1)* WM1( I,J) + UPD( 1)* WM1( I+1,J))
C
2 + ( UT**2 )* ( UPA( 2)* WM1( I-2,J) + UPB( 2)* WM1( I-1,J)
3 + UPC( 2)* WM1( I,J) + UPD( 2)* WM1( I+1,J))/ 2.D0
C
5 - ( UT**3 )* ( UPA( 3)* WM1( I-2,J) + UPB( 3)* WM1( I-1,J)
6 + UPC( 3)* WM1( I,J) + UPD( 3)* WM1( I+1,J))/ 6.D0
C
GO TO 150
C
130 A1 = HX(I-1)/( HX( I-2)* ( HX( I-2) + HX(I-1)))
B1 = ( HX( I-2) + HX(I-1))/( HX( I-2)* HX(I-1))
C1 = ( HX( I-2) + 2.D0* HX(I-1))/( HX(I-1)
1 * ( HX( I-2) + HX(I-1)))
C
A2 = 1.D0 /( HX( I-2)* ( HX( I-2) + HX(I-1)))
B2 = 1.D0 /( HX( I-2)* HX(I-1))
C2 = 1.D0 /( HX( I-1)* ( HX( I-2) + HX(I-1)))
C
UPW = - UT* ( A1* WM1( I-2,J) - B1* WM1( I-1,J) + C1* WM1( I,J))
1 + ( UT** 2 )* ( A2* WM1( I-2,J) - B2* WM1( I-1,J)
2 + C2* WM1( I,J))
C
WP1( I,J) = WP1( I,J) + UPW
GO TO 200
C
90 IF ( EPT( I,J,4)) GO TO 170
IF ( EPT( I,J,3)) GO TO 333
IF ( EPT( I+1,J,3)) GO TO 190
IF ( EPT( I+2,J,3)) GO TO 180
C
UPW = - UT* ( UAM( I,1)* WM1( I-1,J) + UBM( I,1)* WM1( I,J)
1 + UCM( I,1)* WM1( I+1,J) + UDM( I,1)* WM1( I+2,J)
2 + UEM( I,1)* WM1( I+3,J))
C
3 + ( UT** 2)* ( UAM( I,2)* WM1( I-1,J) + UBM( I,2)* WM1( I,J)
4 + UCM( I,2)* WM1( I+1,J) + UDM( I,2)* WM1( I+2,J)
5 + UEM( I,2)* WM1( I+3,J))/ 2.D0
C
6 - ( UT** 3)* ( UAM( I,3)* WM1( I-1,J) + UBM( I,3)* WM1( I,J)
7 + UCM( I,3)* WM1( I+1,J) + UDM( I,3)* WM1( I+2,J)
8 + UEM( I,3)* WM1( I+3,J))/ 6.D0
C
9 + ( UT** 4)* ( UAM( I,4)* WM1( I-1,J) + UBM( I,4)* WM1( I,J)
A + UCM( I,4)* WM1( I+1,J) + UDM( I,4)* WM1( I+2,J)
B + UEM( I,4)* WM1( I+3,J))/ 24.D0
GO TO 220
C
170 A1 = ( 2.D0* HX( I) + HX( I+1))/( HX( I)* ( HX( I) + HX( I+1)))
B1 = ( HX( I) + HX( I+1))/( HX( I)* HX( I+1))
C1 = HX( I)/( HX( I+1)* ( HX( I) + HX( I+1)))
C
A2 = 1.D0 /( HX( I)* ( HX( I) + HX( I+1)))
B2 = 1.D0 /( HX( I)* HX( I+1))
C2 = 1.D0 /( HX( I+1)* ( HX( I) + HX( I+1)))
C
UPW = - UT* ( - A1* WM1(I,J) + B1* WM1(I+1,J) - C1* WM1(I+2,J))
1 + ( UT** 2)* ( A2* WM1(I,J) - B2* WM1(I+1,J) + C2* WM1(I+2,J))
C
WP1( I,J) = WP1( I,J) + UPW
GO TO 200
C
180 CONTINUE
C
C ------------------------------------------------
DO 230 K = 1, 3
C
HD = HX( I+1)* ( HX( I) + HX( I+1))* ( HX(I-1)
1 + HX( I) + HX( I+1))
HC = HX( I)* ( HX(I-1) + HX( I))
HA = HX(I-1)
C
UPD( K) = ( - BB( K)* HX(I-1)* HX( I)
1 - CC( K)* ( HX( I) - HX(I-1)) + DD( K))/ HD
C
UPC( K) = ( - UPD( K)* ( HX( I) + HX( I+1))
1 * ( HX(I-1) + HX( I) + HX( I+1))
2 + BB( K)* HX(I-1) + CC( K))/ HC
C
UPA( K) = ( UPD( K)* ( HX( I) + HX( I+1))
1 + UPC( K)* HX( I) - BB( K))/ HA
C
UPB( K) = AA( K) - ( UPA( K) + UPC( K) + UPD( K))
C
230 CONTINUE
C ------------------------------------------------
C
UPW = - UT* ( UPA( 1)* WM1( I-1,J)
1 + UPB( 1)* WM1( I,J)+ UPC( 1)* WM1( I+1,J)
2 + UPD( 1)* WM1( I+2,J))
C
3 + ( UT** 2 )* ( UPA( 2)* WM1( I-1,J)
4 + UPB( 2)* WM1( I,J) + UPC( 2)* WM1( I+1,J)
5 + UPD( 2)* WM1( I+2,J))/ 2.D0
C
6 - ( UT** 3 )* ( UPA( 3)* WM1( I-1,J)
7 + UPB( 3)* WM1( I,J) + UPC( 3)* WM1( I+1,J)
8 + UPD( 3)* WM1( I+2,J))/ 6.D0
GO TO 220
C
190 UPW = - UT* ( - WM1( I,J) + WM1( I+1,J))/ HX( I)
GO TO 220
C
333 WP1( I,J) = WP1( I,J) + UT* WM1( I,J)/ HX(I-1)
GO TO 200
C
220 CONTINUE
C
150 CONTINUE
C
WP1( I,J) = WP1( I,J) + WEI* CENT + ( 1.D0 - WEI )* UPW
C
200 CONTINUE
C ----------------------------------------------------------
IF ( LTC.NE.1) GO TO 950
C
C -----------------------------------------------
DO 800 J = 1, IY
DO 800 I = 1, IX
C
IF ( .NOT.WPT( I,J,4)) GO TO 800
C
IF ( EPT( I,J,1)) WP1( I,J) = WP1( I,J-1)
IF ( EPT( I,J,2)) WP1( I,J) = WP1( I,J+1)
IF ( EPT( I,J,3)) WP1( I,J) = WP1( I-1,J)
IF ( EPT( I,J,4)) WP1( I,J) = WP1( I+1,J)
C
800 CONTINUE
C ------------------------------------------------
C
950 CONTINUE
C
RETURN
END
C **********************************************************************
C
SUBROUTINE WVCL ( WP1, WM1, PS1, WPT, EPT, V, IX, IY,
1 VNA, VNB, VNC, VND, VNE, VAP, VBP, VCP,
2 VDP, VEP, VAM, VBM, VCM, VDM, VEM,
3 VF1, VF2, VF3, HX, HY )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION WP1( IX,IY), WM1( IX,IY), PS1( IX,IY), V( IX,IY)
C
DIMENSION AA ( 4), BB ( 4), CC ( 4), DD ( 4), EE ( 4),
1 UPA( 4), UPB( 4), UPC( 4), UPD( 4)
C
LOGICAL*1 WPT( IX,IY,4) , EPT( IX,IY,4)
C
LOGICAL PRS
COMMON /IDT/ DT, VIS, DX, DY, DX2, DY2, ITT, CON, PRS, RLL
C
COMMON /WEH/ WEI, IHS, LTC
C
DIMENSION VNA( IY,4), VNB( IY,4), VNC( IY,4), VND( IY,4),
1 VNE( IY,4)

DIMENSION VAP( IY,4), VBP( IY,4), VCP( IY,4), VDP( IY,4),
1 VEP( IY,4)

DIMENSION VAM( IY,4), VBM( IY,4), VCM( IY,4), VDM( IY,4),
1 VEM( IY,4)

DIMENSION VF1( IX), VF2( IX), VF3( IX), HX( IX-1), HY( IY-1)
C
VR = 0.D0
VL = 0.D0
C
C --------------------------------------
DO 100 K = 1, 4
C
AA( K) = 0.D0
BB( K) = 0.D0
CC( K) = 0.D0
DD( K) = 0.D0
EE( K) = 0.D0
C
100 CONTINUE
C
BB( 1) = 1.D0
CC( 2) = 2.D0
DD( 3) = 6.D0
EE( 4) = 24.D0
C
C --------------------------------------
DO 200 J = 1, IY
DO 200 I = 1, IX
C --------------------------------------
C
IF ( WPT( I,J,1)) THEN
IF (.NOT.EPT( I,J,3)) GO TO 10
C
V( I,J) = V( I-1,J)
C V( I,J) = 0.D0 ! Right
GO TO 200
C
10 IF (.NOT.EPT( I,J,4)) GO TO 30
V( I,J) = 0.D0 ! Left
GO TO 200
C
30 CONTINUE
C
C ----- Interior ---
C
V( I,J) = - ( VF1( I)* PS1( I-1,J) + VF2( I)* PS1( I,J)
1 + VF3( I)* PS1( I+1,J))
GO TO 115
C
ELSE
IF ( .NOT.EPT( I,J,1)) GO TO 40 ! Upper
V( I,J) = 0.D0
GO TO 200
C
40 IF ( .NOT.EPT( I,J,2)) GO TO 50 ! Below
V( I,J) = 0.D0
GO TO 200
C
50 IF ( .NOT.EPT( I,J,4)) GO TO 61
V( I,J) = 0.D0 ! Step
GO TO 115
C
61 CONTINUE
END IF
C
115 CONTINUE
C
VT = V( I,J)* DT
IF ( EPT( I,J,2)) GO TO 80
IF ( EPT( I,J,1)) GO TO 80
IF ( EPT( I,J-1,2)) GO TO 29
IF ( EPT( I,J+1,1)) GO TO 29
C
CENT = - VT* ( VNA( J,1)* WM1( I,J-2) + VNB( J,1)* WM1( I,J-1)
1 + VNC( J,1)* WM1( I,J) + VND( J,1)* WM1( I,J+1)
2 + VNE( J,1)* WM1( I,J+2))
C
3 + ( VT** 2)* ( VNA( J,2)* WM1( I,J-2) + VNB( J,2)* WM1( I,J-1)
4 + VNC( J,2)* WM1( I,J) + VND( J,2)* WM1( I,J+1)
5 + VNE( J,2)* WM1( I,J+2))/ 2.D0
C
6 - ( VT** 3)* ( VNA( J,3)* WM1( I,J-2) + VNB( J,3)* WM1( I,J-1)
7 + VNC( J,3)* WM1( I,J) + VND( J,3)* WM1( I,J+1)
8 + VNE( J,3)* WM1( I,J+2))/ 6.D0
C
9 + ( VT** 4)* ( VNA( J,4)* WM1( I,J-2) + VNB( J,4)* WM1( I,J-1)
A + VNC( J,4)* WM1( I,J) + VND( J,4)* WM1( I,J+1)
B + VNE( J,4)* WM1( I,J+2))/ 24.D0
C
GO TO 80
C
29 CONTINUE
C
A1 = HY( J)/( HY( J-1)* ( HY( J-1) + HY( J)))
B1 = ( HY( J-1)-HY( J))/( HY( J-1)* HY( J))
C1 = HY( J-1)/( HY( J)* ( HY( J-1) + HY( J)))
C
A2 = 1.D0/( HY( J-1)* ( HY( J-1) + HY( J)))
B2 = 1.D0/( HY( J-1)* HY( J))
C2 = 1.D0/( HY( J)* ( HY( J-1) + HY( J)))
C
CENT = - VT* ( -A1* WM1(I,J-1) - B1* WM1(I,J) + C1* WM1(I,J+1))
C
1 + ( VT** 2)*( A2* WM1(I,J-1) - B2* WM1(I,J) + C2* WM1(I,J+1))
C
80 CONTINUE
C
IF ( V( I,J).LT.0) GO TO 90
IF ( EPT( I,J,1)) GO TO 130
IF ( EPT( I,J,2)) GO TO 900
IF ( EPT( I,J-1,2)) GO TO 110
IF ( EPT( I,J-2,2)) GO TO 120
C
UPW = - VT* ( VAP( J,1)* WM1( I,J-3) + VBP( J,1)* WM1( I,J-2)
1 + VCP( J,1)* WM1( I,J-1) + VDP( J,1)* WM1( I,J)
2 + VEP( J,1)* WM1( I,J+1))
C
3 + (VT** 2)* ( VAP( J,2)* WM1( I,J-3) + VBP( J,2)* WM1( I,J-2)
4 + VCP( J,2)* WM1( I,J-1) + VDP( J,2)* WM1( I,J)
5 + VEP( J,2)* WM1( I,J+1))/ 2.D0
C
6 - (VT** 3)* ( VAP( J,3)* WM1( I,J-3) + VBP( J,3)* WM1( I,J-2)
7 + VCP( J,3)* WM1( I,J-1) + VDP( J,3)* WM1( I,J)
8 + VEP( J,3)* WM1( I,J+1))/ 6.D0
C
9 + (VT** 4)* ( VAP( J,4)* WM1( I,J-3) + VBP( J,4)* WM1( I,J-2)
A + VCP( J,4)* WM1( I,J-1) + VDP( J,4)* WM1( I,J)
B + VEP( J,4)* WM1( I,J+1))/ 24.D0
GO TO 150
C
900 WP1( I,J) = WP1( I,J) - VT* WM1( I,J)/ HY( J)
GO TO 200
C
110 UPW = - VT* ( - WM1( I,J-1) + WM1( I,J))/ HY( J-1)
GO TO 150
C
120 CONTINUE
C ------------------------------------------------
DO 160 K=1, 3
C
HA = HY( J-2)*( HY( J-2) + HY( J-1))
1 *( HY( J-2) + HY( J-1) + HY( J))
C
HB = HY( J-1)*( HY( J-1) + HY( J))
C
HD = HY( J)
C
UPA( K) = ( BB( K)* HY( J-1)* HY( J)
1 + CC( K)* ( HY( J) - HY( J-1)) - DD( K))/ HA
C
UPB( K) = ( - UPA( K)* ( HY( J-2) + HY( J-1))
1 * ( HY( J-2) + HY( J-1) + HY( J))
2 - BB( K)* HY( J) + CC( K))/ HB
C
UPD( K) = ( UPA( K)* ( HY( J-2) + HY( J-1)) + UPB( K)* HY( J-1)
1 + BB( K))/ HD
C
UPC( K) = AA( K) - ( UPA( K) + UPB( K) + UPD( K))
C
160 CONTINUE
C ------------------------------------------------
C
UPW =
C
1 - VT* ( UPA( 1)* WM1( I,J-2) + UPB( 1)* WM1( I,J-1)
2 + UPC( 1)* WM1( I,J) + UPD( 1)* WM1( I,J+1))
C
3 + (VT** 2 )* ( UPA( 2)* WM1( I,J-2) + UPB( 2)* WM1( I,J-1)
4 + UPC( 2)* WM1( I,J) + UPD( 2)* WM1( I,J+1))/ 2.D0
C
5 - (VT** 3 )* ( UPA(3)* WM1( I,J-2)+ UPB(3)* WM1( I,J-1)
6 + UPC(3)* WM1( I,J) + UPD(3)* WM1( I,J+1))/ 6.D0
GO TO 150
C
130 A1 = HY( J-1)/( HY( J-2)* ( HY( J-2) + HY( J-1)))
B1 = ( HY( J-2) + HY( J-1))/( HY( J-2)* HY( J-1))
C1 = ( HY( J-2) + 2.D0 * HY( J-1))
1 /( HY( J-1)* ( HY( J-2) + HY( J-1)))
C
A2 = 1.D0 /( HY( J-2)* ( HY( J-2) + HY( J-1)))
B2 = 1.D0 /( HY( J-2)* HY( J-1))
C2 = 1.D0 /( HY( J-1)* ( HY( J-2) + HY( J-1)))
C
UPW = - VT*( A1* WM1( I,J-2) - B1* WM1( I,J-1) + C1* WM1( I,J))
1 + ( VT** 2)* ( A2* WM1( I,J-2) - B2* WM1( I,J-1)
2 + C2* WM1( I,J))
C
WP1( I,J) = WP1( I,J) + UPW
GO TO 200
C
90 IF ( EPT( I,J,2)) GO TO 170
IF ( EPT( I,J,1)) GO TO 333
IF ( EPT( I,J+1, 1)) GO TO 190
IF ( EPT( I,J+2, 1)) GO TO 180
C
UPW = - VT * ( VAM( J,1)* WM1( I,J-1) + VBM( J,1)* WM1( I,J)
1 + VCM( J,1)* WM1( I,J+1) + VDM( J,1)* WM1( I,J+2)
2 + VEM( J,1)* WM1( I,J+3))
C
3 + ( VT**2 )* ( VAM( J,2)* WM1( I,J-1) + VBM( J,2)* WM1( I,J)
4 + VCM( J,2)* WM1( I,J+1) + VDM( J,2)* WM1( I,J+2)
5 + VEM( J,2)* WM1( I,J+3))/ 2.D0
C
6 - ( VT**3)* ( VAM( J,3)* WM1( I,J-1) + VBM( J,3)* WM1( I,J)
7 + VCM( J,3)* WM1( I,J+1) + VDM( J,3)* WM1( I,J+2)
8 + VEM( J,3)* WM1( I,J+3))/ 6.D0
C
9 + ( VT**4)* ( VAM( J,4)* WM1( I,J-1) + VBM( J,4)* WM1( I,J)
A + VCM( J,4)* WM1( I,J+1) + VDM( J,4)* WM1( I,J+2)
B + VEM( J,4)* WM1( I,J+3))/ 24.D0
GO TO 220
C
170 CONTINUE
C
A1 = ( 2.D0* HY( J) + HY( J+1))/( HY( J)* ( HY( J) + HY( J+1)))
B1 = ( HY( J) + HY( J+1))/( HY( J)* HY( J+1))
C1 = HY( J)/( HY( J+1)* ( HY( J) + HY( J+1)))
C
A2 = 1.D0/( HY( J)* ( HY( J) + HY( J+1)))
B2 = 1.D0/( HY( J)* HY( J+1))
C2 = 1.D0/( HY( J+1)* ( HY( J) + HY( J+1)))
C
UPW = - VT* ( - A1* WM1(I,J) + B1* WM1(I,J+1) - C1* WM1(I,J+2))
1 + ( VT** 2)* ( A2* WM1(I,J) - B2* WM1(I,J+1) + C2* WM1(I,J+2))
C
WP1( I,J) = WP1( I,J) + UPW
GO TO 200
C
180 DO 230 K = 1, 3
C
HD = HY( J+1)* ( HY(J) + HY(J+1))* ( HY(J-1) + HY(J) + HY(J+1))
HC = HY( J) * ( HY( J-1) + HY( J))
HA = HY( J-1)
C
UPD( K) = ( - BB( K)* HY( J-1)* HY( J)
1 - CC( K)* ( HY( J) - HY( J-1)) + DD( K))/ HD
C
UPC( K) = ( - UPD( K)* ( HY( J) + HY( J+1))
1 * ( HY( J-1) + HY( J) + HY( J+1))
2 + BB( K)* HY( J-1) + CC( K))/ HC
C
UPA( K) = ( UPD( K)* ( HY( J) + HY( J+1))
1 + UPC( K)* HY( J) - BB( K))/ HA
C
UPB( K) = AA( K) - ( UPA( K) + UPC( K) + UPD( K))
C
230 CONTINUE
C
UPW = - VT* ( UPA( 1)* WM1( I,J-1) + UPB( 1)* WM1( I,J)
1 + UPC( 1)* WM1( I,J+1) + UPD( 1)* WM1( I,J+2))
C
2 + ( VT**2 )* ( UPA( 2)* WM1( I,J-1) + UPB( 2)* WM1( I,J)
3 + UPC( 2)* WM1( I,J+1) + UPD( 2)* WM1( I,J+2))/2.D0
C
4 - ( VT**3)* ( UPA( 3)* WM1( I,J-1) + UPB( 3)* WM1( I,J)
5 + UPC( 3)* WM1( I,J+1) + UPD( 3)* WM1( I,J+2))/6.D0
GO TO 220
C
190 UPW = - VT* ( - WM1( I,J) + WM1( I,J+1))/ HY( J)
GO TO 220
C
333 WP1( I,J) = WP1( I,J) + VT* WM1( I,J)/ HY( J-1)
GO TO 200
C
220 CONTINUE
C
150 CONTINUE
C
WP1( I,J) = WP1( I,J) + WEI* CENT + ( 1.D0 - WEI )* UPW
C
200 CONTINUE
C --------------------------------------
C
IF ( LTC.NE.1) GO TO 605
C
C --------------------------------------
DO 800 J = 1, IY
DO 800 I = 1, IX
C
IF ( .NOT.WPT( I,J,4 )) GO TO 800
C
IF ( EPT( I,J,1)) WP1( I,J) = WP1( I,J-1)
IF ( EPT( I,J,2)) WP1( I,J) = WP1( I,J+1)
IF ( EPT( I,J,3)) WP1( I,J) = WP1( I-1,J)
IF ( EPT( I,J,4)) WP1( I,J) = WP1( I+1,J)
C
800 CONTINUE
C --------------------------------------
C
605 CONTINUE
C
RETURN
END
C **********************************************************************
C
SUBROUTINE VARV ( WP1, WB, WD, VRV, II,JJ, IX, IY )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION WP1( IX,IY), WB( IX,IY), WD( IX,IY)
C
C ---------------------------
DO 100 J = 1, IY
DO 100 I = 1, IX
C ---------------------------
C
WD( I,J) = 0.D0
C
100 CONTINUE
C ---------------------------
C
WL1 = 0.D0
WL2 = 0.D0
WDL = 0.D0
C
ABSW = 0.D0
VRV = 0.D0
C
C ---------------------------
DO 200 J = 1, IY
DO 200 I = 1, IX
C ---------------------------
C
WZS = DABS( WP1( I,J))
IF ( WZS.GT.WL1) WL1 = WZS
C
200 CONTINUE
C ---------------------------
WL2 = WL1* 0.001D0
C
C ------------------------------------------------
DO 300 J = 1, IY
DO 300 I = 1, IX
C ------------------------------------------------
C
ABSW = DABS( WP1( I,J))
IF ( ABSW.LE.WL2) GO TO 300
WD( I,J) = DABS( ( WP1( I,J) - WB( I,J))/ WP1( I,J))
IF ( WD( I,J).GT.WDL ) THEN
WDL = WD( I,J)
II = I
JJ = J
END IF
C
300 CONTINUE
C ------------------------------------------------
VRV = WDL* 100.D0 ! (%)
C
RETURN
END
C **********************************************************************
C
SUBROUTINE VRVR ( WP1, WB, WD, II, JJ, W, PS, X1R, PS1, IX, IY,
1 * )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
COMMON /CTL/ LOP, NPR, ILP, NNN
COMMON /DDS/ DS1, DS2, DS3, DSV( 10), DIS( 100)
C
DIMENSION WP1( IX,IY), WB ( IX,IY), WD ( IX,IY), W( IX,IY),
1 PS ( IX,IY), X1R( 100000), PS1( IX,IY)
C
C ---------------------------
DO 100 J = 1, IY
DO 100 I = 1, IX
C
W ( I,J) = PS ( I,J)
PS( I,J) = PS1( I,J)
C
100 CONTINUE
C ---------------------------
C
C ----- X1R ---
C
IF ( LOP.EQ.1) DS3 = 0.D0
X1R( LOP) = DS3
C
C ----- VRVR ---
C
CALL VARV ( WP1, WB, WD, VRV, II, JJ, IX, IY )
C
C ----------------------------------------------
C
IF ( LOP/NPR* NPR .EQ. LOP) THEN
C !
C ----- NNN : No. of plot timing ---
NNN = NNN + 1
WRITE(6, 2000) LOP, VRV, II, JJ, X1R( LOP)
WRITE(17,2100) LOP, VRV, II, JJ, X1R( LOP)
C
END IF
C
2000 FORMAT(' * LOP =', I6,' Var.of(W)=',F10.3,' (%) I,J of Max.',
1 2I4,' X1=',F6.3)
2100 FORMAT( I6,F9.3,2I4,F6.3)
C
C ----------------------------
DO 200 J = 1, IY
DO 200 I = 1, IX
C
WB( I,J) = WP1( I,J)
200 CONTINUE
C ----------------------------
C
C ----- Stationarity Check ---
C
IF ( VRV.LE.0.1D0) THEN
WRITE(6,2200) LOP
2200 FORMAT(' *** Stationary Condition Reached *** at Loop =', I6)
C
NNN = NNN + 1
WRITE(6, 2000) LOP, VRV, II, JJ, X1R( LOP)
WRITE(17,2100) LOP, VRV, II, JJ, X1R( LOP)
END IF
C -----------------------------
C
IF ( VRV.LE.0.1D0 .OR. LOP.EQ.ILP) RETURN 1
C
RETURN
END
C **********************************************************************
C
SUBROUTINE WVCF ( VNA, VNB, VNC, VND, VNE, VAP, VBP, VCP,
1 VDP, VEP, VAM, VBM, VCM, VDM, VEM, VF1,
2 VF2, VF3, HX, HY, IX, IY )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION AA( 4), BB( 4), CC( 4), DD( 4), EE( 4)
C
DIMENSION VNA( IY,4), VNB( IY,4), VNC( IY,4), VND( IY,4),
1 VNE( IY,4)
C
DIMENSION VAP( IY,4), VBP( IY,4), VCP( IY,4), VDP( IY,4),
1 VEP( IY,4)
C
DIMENSION VAM( IY,4), VBM( IY,4), VCM( IY,4), VDM( IY,4),
1 VEM( IY,4)
C
DIMENSION VF1( IX), VF2( IX), VF3( IX), HX( IX-1), HY( IY-1)
C
C ------------------------
C
DO 100 K = 1, 4
C
AA( K) = 0.D0
BB( K) = 0.D0
CC( K) = 0.D0
DD( K) = 0.D0
EE( K) = 0.D0
C
100 CONTINUE
C
BB( 1) = 1.D0
CC( 2) = 2.D0
DD( 3) = 6.D0
EE( 4) = 24.D0
C
C --------------------------------------
DO 200 J = 1, IY
DO 200 K = 1, 4
C --------------------------------------
C
VNA( J,K) = 0.D0
VNB( J,K) = 0.D0
VNC( J,K) = 0.D0
VND( J,K) = 0.D0
VNE( J,K) = 0.D0
C
VAP( J,K) = 0.D0
VBP( J,K) = 0.D0
VCP( J,K) = 0.D0
VDP( J,K) = 0.D0
VEP( J,K) = 0.D0
C
VAM( J,K) = 0.D0
VBM( J,K) = 0.D0
VCM( J,K) = 0.D0
VDM( J,K) = 0.D0
VEM( J,K) = 0.D0
C
200 CONTINUE
C --------------------------------------
DO 300 I = 1, IX
C
VF1( I) = 0.D0
VF2( I) = 0.D0
VF3( I) = 0.D0
C
300 CONTINUE
C --------------------------------------
DO 400 I = 2, IX-1
C
VF1( I) = - HX( I)/( HX(I-1)* ( HX(I-1) + HX( I)))
VF2( I) = ( HX( I) - HX(I-1))/( HX(I-1)* HX( I))
VF3( I) = HX(I-1)/( HX( I)* ( HX(I-1) + HX( I)))
C
400 CONTINUE
C
C --------------------------------------
DO 500 J = 3, IY-2
C
HA = HY( J-2)* ( HY( J-2) + HY( J-1))
1 * ( HY( J-2) + HY( J-1) + HY( J))
2 * ( HY( J-2) + HY( J-1) + HY( J) + HY( J+1))
C
HB = HY( J-1)*( HY( J-1)+HY( J))*( HY( J-1)+HY( J)+HY( J+1))
C
HD = HY( J)* HY( J+1)
HE = HY( J) + HY( J+1)
C
C ----------------------------
DO 700 K = 1,4
C
VNA( J,K) = ( BB( K)* HY( J-1)* HY( J)* ( HY( J) + HY( J+1))
1 - CC( K)* ( HY( J-1)* ( 2.D0 * HY( J) + HY( J+1))
2 - HY( J)* ( HY( J) + HY( J+1)))
C
3 + DD( K)* ( HY( J-1) - 2.D0 * HY( J) - HY( J+1))
4 + EE( K))/ HA
C
VNB( J,K) = ( - VNA( J,K)* ( HY( J-2) + HY( J-1))
1 * ( HY( J-2) + HY( J-1) + HY( J))
2 * ( HY( J-2) + HY( J-1) + HY( J) + HY( J+1))
C
3 - BB( K)* HY( J)* ( HY( J) + HY( J+1))
4 + CC( K)* ( 2.D0 * HY( J) + HY( J+1)) - DD( K))/ HB
C
VND( J,K) = ( VNA( J,K)* ( HY( J-2) + HY( J-1))
1 * ( HY( J-2) + HY( J-1) + HY( J) + HY( J+1))
C
2 + VNB( J,K)* HY( J-1)* ( HY( J-1) + HY( J) + HY( J+1))
3 + BB( K)* ( HY( J) + HY( J+1)) - CC( K))/ HD
C
VNE( J,K) = ( VNA( J,K)* ( HY( J-2) + HY( J-1))
1 + VNB( J,K)* HY(J-1) - VND(J,K)* HY(J) + BB( K))/ HE
C
VNC( J,K) = AA( K)
1 - ( VNA( J,K) + VNB( J,K) + VND( J,K) + VNE( J,K))
C
700 CONTINUE
C ---------------------------
C
500 CONTINUE
C --------------------------------------
DO 800 J = 4, IY-1
C
HA = HY( J-3)* ( HY( J-3) + HY( J-2))
1 * ( HY( J-3) + HY( J-2) + HY( J-1))
2 * ( HY( J-3) + HY( J-2) + HY( J-1) + HY( J))
C
HB = HY( J-2)*( HY( J-2)+HY( J-1))*( HY( J-2)+HY( J-1)+HY( J))
C
HC = HY( J-1)* ( HY( J-1) + HY( J))
HE = HY( J)
C
C ----------------------------
DO 900 K = 1, 4
C
VAP( J,K) = ( - BB( K)*HY( J-1)*HY( J)*( HY( J-2) + HY( J-1))
1 + CC( K)* ( HY( J-2)* HY( J-1) + ( HY( J-1)** 2)
C
2 - HY( J-2)* HY( J) - 2.D0 * HY( J-1)* HY( J))
3 + DD( K)* ( HY( J-2) + 2.D0 * HY( J-1) - HY( J))
4 + EE( K))/ HA
C
VBP( J,K) = ( - VAP( J,K)*( HY( J-3) + HY( J-2))
1 * ( HY( J-3) + HY( J-2) + HY( J-1))
2 * ( HY( J-3) + HY( J-2) + HY( J-1) + HY( J))
C
3 + BB( K)*HY( J-1)*HY( J)
4 - CC( K)*( HY( J-1) - HY( J)) - DD( K))/ HB
C
VCP( J,K) = ( - VAP( J,K)
1 * ( HY( J-3) + HY( J-2) + HY( J-1))
2 * ( HY( J-3) + HY( J-2) + HY( J-1) + HY( J))
C
3 - VBP( J,K)* ( HY( J-2) + HY( J-1))
4 * ( HY( J-2) + HY( J-1) + HY( J))
5 - BB( K)* HY( J) + CC( K))/ HC
C
VEP( J,K) = ( VAP( J,K)* ( HY( J-3)+ HY( J-2) + HY( J-1))
1 + VBP( J,K)* ( HY( J-2) + HY( J-1))
2 + VCP( J,K)* HY( J-1) + BB( K))/ HE
C
VDP( J,K) = AA( K)
1 - ( VAP( J,K) + VBP( J,K) + VCP( J,K) + VEP( J,K))
C
900 CONTINUE
C ----------------------------
C
800 CONTINUE
C ------------------------------------------------
DO 950 J = 2, IY-3
C
HE = HY( J+2)*( HY( J+1) + HY( J+2))
1 *( HY( J) + HY( J+1) + HY( J+2))
2 *( HY( J-1) + HY( J) + HY( J+1) + HY( J+2))
C
HD = HY(J+1)* ( HY(J) + HY(J+1))* ( HY(J-1) + HY(J) + HY(J+1))
C
HC = HY( J)* ( HY( J-1) + HY( J))
C
HA = HY( J-1)
C
C --------------------------------------
DO 970 K = 1, 4
C
VEM( J,K) = ( BB( K)* HY( J-1)* HY( J)* ( HY( J)
1 + HY( J+1)) + CC( K)* ( ( HY( J)** 2)
2 + HY( J)* HY( J+1) - HY( J-1)* HY( J+1)
C
3 - 2.D0 * HY( J-1)* HY( J)) + DD( K)* ( HY( J-1)
4 - 2.D0 * HY( J) - HY( J+1)) + EE( K))/ HE
C
VDM( J,K) = ( - VEM( J,K)* ( HY( J+1) + HY( J+2))
1 * ( HY( J) + HY( J+1) + HY( J+2))
2 * ( HY( J-1) + HY( J) + HY( J+1) + HY( J+2))
C
3 - BB( K)* HY( J-1)* HY( J)
4 - CC( K)* ( HY( J) - HY( J-1)) + DD( K))/ HD
C
VCM( J,K) = ( - VEM( J,K)* ( HY( J)+ HY( J+1) + HY( J+2))
1 * ( HY( J-1) + HY( J) + HY( J+1) + HY( J+2))
2 - VDM( J,K)* ( HY( J) + HY( J+1))
C
3 * ( HY( J-1) + HY( J) + HY( J+1))
4 + BB( K)* HY( J-1) + CC( K))/ HC
C
VAM( J,K) = ( VEM( J,K)* ( HY( J) + HY( J+1) + HY( J+2))
1 + VDM( J,K)* ( HY( J) + HY( J+1))
2 + VCM( J,K)* HY( J) - BB( K))/ HA
C
VBM( J,K) = AA( K)
1 - ( VAM( J,K) + VCM( J,K) + VDM( J,K) + VEM( J,K))
C
970 CONTINUE
C --------------------------------------
C
950 CONTINUE
C ------------------------------------------------
C
RETURN
END
C **********************************************************************