戻る

PROGRAM MAIN
C
C ******************************************************************** C
C C
C C V 2 D - F D M C
C C
C ( Convection Equation Solver ) C
C C
C ( The 4th & 2nd orderd Finite Difference Method ) C
C C
C ( Version 5.0 ) C
C C
C Copyright : Yasuhiro Matsuda C
C C
C ******************************************************************** C
C
IMPLICIT REAL*8(A-H,O-Z)
C
C ----- DATA ( 1) -------------------------------------
C
** PARAMETER ( IX = 21, IY = 21,
C
** PARAMETER ( IX = 51, IY = 51,
C
PARAMETER ( IX = 101, IY = 101,
C
1 NNP = IX* IY, NEL = ( IX-1)* ( IY-1))
C -----------------------------------------------------
C
DIMENSION WVL( IX,IY), WS( IX,IY), U ( IX,IY), V ( IX,IY),
1 CE ( IX,IY), XE( IX,IY), YE( IX,IY), DIF( IX,IY)
C
DIMENSION XX ( 10201), YY ( 10201), WW ( 10201),
1 NDE( 10000,4), PSX( 50000), PSY( 50000)
C
LOGICAL*1 WPT( IX,IY,4)
C
C ----- For Analytical Solution ---
C
DIMENSION AWW( 10201) ! Max.: 100* 100
C -----------------------------------------------
C
COMMON /CNST/ WMX, ERV, ERL, CMX, DDT, RLA, CMA, CMN, MSC
COMMON /DAT2/ DX2, FS, CF, DY, DY2, TM, R, AL, BT, CFI
COMMON /DAT1/ DX, F, DT, NLT, NPR, UL, UU, D, VIS, WEI,
1 ICR, ISC
C
CHARACTER*6 MSC
C
OPEN ( 7, FILE='CASE.DATA')
C
C ----- DATA ( 2) -------------------------------------
C
** OPEN ( 8, FILE='CV2D_20.DATA' ) ! 20* 20 DT = 0.05
** OPEN ( 10, FILE='2ND_CV2D_20_Fig.DATA' )
** OPEN ( 10, FILE='4TH_CV2D_20_Fig.DATA' )
C
C ---------------------------
C
** OPEN ( 8, FILE='CV2D_50.DATA' )
C
** OPEN ( 10, FILE='2ND_CV2D_50_Fig.DATA' ) ! 50* 50 DT = 0.005
** OPEN ( 10, FILE='4TH_CV2D_50_Fig.DATA' ) ! 50* 50 DT = 0.005
C
C ---------------------------
C
OPEN ( 8, FILE='CV2D_100.DATA' ) ! 100* 100 DT = 0.0025 0.0005
C
** OPEN ( 10, FILE='2ND_CV2D_100_Fig.DATA' )
OPEN ( 10, FILE='4TH_CV2D_100_Fig.DATA' )
C
C -----------------------------------------------------
C
KCP = 0
C
9999 LOP = 0
C
CALL ISB1 ( WVL, WS, WPT, IX, IY, XE, YE, CE, U, V,
1 NNP, XX, YY, WW )
C
CALL ISB2 ( U, V, XE, YE, CE, WVL, XX, YY, WW, PSX,
1 PSY, NDE, IX, IY, LOP, NNP, NEL, II, JJ )
C
C ------------------------------------------------
999 CONTINUE
C
CALL ANSB ( XE, YE, CE, IX, IY, LOP )
C
LOP = LOP + 1
DDT = DDT + DT
C
DO 100 I = 1, IX
DO 100 J = 1, IY
C
WS( I,J) = WVL( I,J)
100 CONTINUE
C
IF ( ISC) 350, 350, 370
350 CONTINUE
C
C ----- 4th Ordered Scheme ---
C
C ----------------------------------------------------------
C
CALL WU4C ( WVL, WS, WPT, IX, IY, U )
C
C ----------------------------------------------------------
C
CALL WV4C ( WVL, WS, WPT, IX, IY, V )
C
C ----------------------------------------------------------
GO TO 390
C
370 CONTINUE
C
C ----- 2nd Ordered Scheme ---
C
C ----------------------------------------------------------
C
CALL WU2C ( WVL, WS, WPT, IX, IY, U )
C
C ----------------------------------------------------------
C
CALL WV2C ( WVL, WS, WPT, IX, IY, V )
C
C ----------------------------------------------------------
C
390 CONTINUE
C
IF ( ICR.EQ.1) CALL WDCP ( U, V, WVL, WS, WPT, IX, IY, DIF )
C ==============================
C
C ----- Error computation ---
C
CALL ERCP ( WVL, CE, IX, IY, LOP )
C
IF ( LOP.GE.NLT) GO TO 1000
C
IF ( MOD( LOP,NPR).EQ.0)
C
1 CALL PLTS ( WVL, XE, YE, CE, XX, YY, WW, AWW, PSX,
2 PSY, IX, IY, LOP, II, JJ, NNP )
C
GO TO 999
C ------------------------------------------------
C
1000 CONTINUE
C
CALL PSMV ( WVL, XE, YE, IX, IY, PSX, PSY, LOP, II, JJ )
C
CALL SBLT ( WVL, CE, XX, YY, WW, AWW, PSX, PSY, DIF, IX,
1 IY, LOP, NNP, KCP, *9999 )
C
CLOSE (10)
CLOSE ( 8)
CLOSE ( 7)
C
STOP
END
C **********************************************************************
C
SUBROUTINE ANSB ( X, Y, CB, IX, IY, LOP )
C
IMPLICIT REAL*8(A-H,O-Z)
C
COMMON /DAT1/ DX, F, DT, NLT, NPR, UL, UU, D, VIS, WEI, ICR, ISC
C
DIMENSION X( IX,IY), Y( IX,IY), CB( IX,IY)
C
PAI = 3.1415926535 8979323846 D0
DXX = 2.D0/ DFLOAT( IX-1)
C
J = 0
999 CONTINUE
K = ( IX-1)/ 2 + 1
J = J + 1
C ------------------------------------------------
DO 100 I = 1, K
C
X( I,J) = -1.D0 + DX* DFLOAT( I-1)
100 CONTINUE
C ------------------------------------------------
DO 200 I = K, IX
C
X( I,J) = DX* DFLOAT( I-K)
200 CONTINUE
C ------------------------------------------------
C
IF ( J.LT.IY) GO TO 999
I = 0
9999 CONTINUE
K = ( IY - 1)/ 2.D0 + 1
I = I + 1
C ------------------------------------------------
DO 300 J = 1, K
C
Y( I,J) = - 1.D0 + DX* DFLOAT( J-1)
300 CONTINUE
C ------------------------------------------------
DO 400 J = K, IY
C
Y( I,J) = DX* DFLOAT( J-K)
400 CONTINUE
C ------------------------------------------------
C
IF ( I.LT.IX) GO TO 9999
T = DFLOAT( LOP)* DT
PX = 0.5D0* DCOS( T - PAI/ 2.D0)
PY = 0.5D0* DSIN( T - PAI/ 2.D0)
C
C --------------------------------------------------------------
DO 500 I = 1, IX
DO 500 J = 1, IY
C
R = DSQRT(( X( I,J) - PX)** 2 + ( Y( I,J) - PY)** 2 )
C ===========================================================
C
CB ( I,J) = 0.5D0 + 0.5D0* DCOS(( 10.D0/ 3.D0)* PAI* R )
IF ( R.GT.0.3D0) CB( I,J) = 0.D0
C
500 CONTINUE
C --------------------------------------------------------------
CBMX = 0.D0
DO 600 I = 1, IX
DO 600 J = 1, IY
C
IF ( CB( I,J).LE.CBMX) GO TO 600
CBMX = CB( I,J)
II = I
JJ = J
C
600 CONTINUE
C ---------------------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE ERCP ( WVL, CE, IX, IY, LOP )
C
C ----- Numerical Error for each LOP ---
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION WVL( IX,IY), CE( IX,IY)
C
COMMON /CNST/ WMX, ERV, ERL, CMX, DDT, RLA, CMA, CMN, MSC
CHARACTER*6 MSC
C
C ----- Attention Loop < = 50000 ---
C
COMMON /ECHR/ ERR( 50000), SCMA
COMMON /DAT1/ DX, F, DT, NLT, NPR, UL, UU, D, VIS, WEI, ICR, ISC
C
NEL2 = ( IX-2)* ( IY-2)
C
N = IDNINT( DFLOAT( IX+1)/ 2.D0)
C
RMS = 0.D0
C ------------------------------------------------
DO 100 I = 1, IX
DO 100 J = 1, IY
C
RMS = RMS + DABS( WVL( I,J) - CE( I,J))
C
100 CONTINUE
C ------------------------------------------------
C
RMS = RMS/ DFLOAT( NEL2)
RLA = RLA + RMS/ CMX
C
ERR( LOP) = RMS* 100.D0
ERV = ERR( LOP)
TRLA = RLA/ DFLOAT( LOP)
C
CMA = - 100.D0
CMN = 100.D0
C
C ------------------------------------------------
DO 200 I = 1, IX
DO 200 J = 1, IY
C
IF ( CMA.LT.WVL( I,J)) CMA = WVL( I,J)
IF ( CMN.GT.WVL( I,J)) CMN = WVL( I,J)
C
200 CONTINUE
C ------------------------------------------------
C
IF ( CMA.GT.1.1D0) THEN
WRITE(6,*) ' '
WRITE(6,*) ' ** CMA. > 1.1 * STOP '
STOP
END IF
C
C ----------------------------------------------------------
IF ( LOP.EQ.NLT) THEN
C
CMA = - 100.D0
CMN = 100.D0
C ------------------------------------------------
DO 300 I = 1, IX
DO 300 J = 1, IY
C
IF ( CMA.LT.WVL( I,J)) CMA = WVL( I,J)
IF ( CMN.GT.WVL( I,J)) CMN = WVL( I,J)
C
300 CONTINUE
C ------------------------------------------------
PRMS = RMS * 100.D0
ERL = TRLA* 100.D0
C
END IF
C ----------------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE INPB ( ISC, WEI, ICR, DT, NLT, NPR, MSC, IX, IY )
C
IMPLICIT REAL*8(A-H,O-Z)
C
CHARACTER*6 MSC
C
COMMON /CNTL/ IFC
C
1 CONTINUE
C
C ----- Attention ---
C
WRITE(6,2000) IX, IY
2000 FORMAT(' *** Comp. Area : 2 x 2 : IX =',I3,' IY =',I3,//,
1 ' * Which FDM ? 4 : (4th) FWA 2 : (2nd) SWA)')
READ(5,*) MTH
WRITE(6,*) ' MTH =', MTH
C
C ------------------------------
IF ( MTH.EQ.4) IFC = 1
IF ( MTH.EQ.2) IFC = 3
C ------------------------------
C
IF ( IFC.EQ.1) GO TO 11 ! FWA-Method
IF ( IFC.EQ.2) GO TO 22
IF ( IFC.EQ.3) GO TO 33 ! SWA-Method
IF ( IFC.EQ.4) GO TO 44
C
GO TO 1
C
11 ISC = 0
ICR = 0
MSC = ' FWA '
GO TO 2
C
22 ISC = 0
ICR = 1
MSC = 'FWA(C)'
GO TO 2
C
33 ISC = 1
ICR = 0
MSC = ' SWA '
GO TO 2
C
44 ISC = 1
ICR = 1
MSC = 'SWA(C)'
GO TO 2
C
2 CONTINUE
WRITE(6,2100)
2100 FORMAT(/,' ** Weighting parameter ? (0〜1)')
C
C -----------------------
C
*** READ(5,*) WEI
C
WEI = 0.5D0
C
C -----------------------
C
WRITE(6,2200) MSC, WEI
2200 FORMAT(' ** Method = ',A6,' * WEI =',F5.2 )
C
IF ( WEI.GT.1.D0) GO TO 2
IF ( WEI.LT.0.D0) GO TO 2
WRITE(6,2300)
2300 FORMAT(/,' ** Time step ?')
*** READ(5,*) DT
DT = 0.005D0
DT = 0.0025D0
C =================
C
NLT = IDNINT( 2.D0* 3.14159 26535 897932 D0/ DT )
C ========
C
NTPR = IDNINT( DFLOAT( NLT)/ 6.D0)
WRITE(6,2400) DT, NLT, NTPR
2400 FORMAT(/,' ** DT =',F8.4,' * Loop Limit =',I5,
1 ' Temp. Est. NTPR=', I5)
C
C ------ Attention Loop < = 50000 ---
C
IF ( NLT.GE.50000) WRITE(6,*) ' *** Error, Loop < = 50000 ***'
IF ( NLT.GE.50000) STOP
C
4 CONTINUE
WRITE(6,2600)
2600 FORMAT(' * Plot Intervals ? (* ZERO for the Last Only)')
** READ(5,*) NPR
C
NPR = NTPR
WRITE(6,*) ' * NPR =', NPR
IF ( NPR.EQ.0) NPR = NLT
IF ( NPR.LT.0) GO TO 4
IF ( NPR.GT.NLT) GO TO 4
C
RETURN
END
C **********************************************************************
C
SUBROUTINE INPD ( WVL, WPT, IX, IY )
C
IMPLICIT REAL*8(A-H,O-Z)
C
LOGICAL*1 OMG, PSI, CNT, SLP, MVE, WPT
C
DIMENSION WVL ( IX,IY), WPT( IX,IY,4),
2 INDX( 5), INDY( 5), INDE( 5)
C
COMMON /CNST/ WMX, ERV, ERL, CMX, DDT, RLA, CMA, CMN, MSC
CHARACTER*6 MSC
COMMON /DAT2/ DX2, FS, CF, DY, DY2, TM, R, AL, BT, CFI
COMMON /DAT1/ DX, F, DT, NLT, NPR, UL, UU, D, VIS, WEI, ICR,
1 ISC
C
NAMELIST /DAB/ R, AL, BT, TM
NAMELIST /DAT/ IDT, OMG, PSI, SLP, MVE, CNT
NAMELIST /DAA/ DX, F, DT, NLT, NPR, UL, UU, D, VIS, WEI, ICR,
1 ISC
C
DX = 2.D0/ DFLOAT( IX-1 ) ! IX * IY meshes
DY = 2.D0/ DFLOAT( IY-1)
F = DY/ DX
C
UL = 1.D0
UU = -1.D0
D = 1.D0
C -------------------------------------------------------------
C
CALL INPB ( ISC, WEI, ICR, DT, NLT, NPR, MSC, IX, IY )
C
C -------------------------------------------------------------
DX2 = DX* DX
DY2 = DY* DY
C
FS = F* F
CF = 2.D00* ( FS + 1.D0)
CFI = 1.D0/ CF
TM = 0.D0
C
AL = DT/ DX
BT = DT/ DY
C
C ----- Attention ---
C
VIS = 0.D0 ! Original Diffusivity for Convection Eq. ---
C ==============================================================
C
FR = VIS* DT/ DX2
C
IF ( VIS) 950, 950, 220
C
220 R = UL* D/ VIS
950 CONTINUE
WRITE(6,2000) DX, F, UL, UU, D
2000 FORMAT(/,' * DX =', F6.3,' F =',F6.3,' UL =',F6.3,' UU =',
1 F6.3,' D =',F6.3 )
WRITE(6,2100) R, AL, BT, FR, TM
2100 FORMAT(' * R =',D10.3,' AL =',F6.3,' BT =',F6.3,' Fr =',F6.3,
1 ' TM =',F6.3 )
C
LED = 4
LOP = 0
IDER = 0
C
DO 100 I = 1, IX
DO 100 J = 1, IY
C
WPT( I,J,1) =.TRUE.
100 CONTINUE
C
C ----- Edge points ---
C
900 CONTINUE
IDT = 0
OMG =.FALSE.
PSI =.FALSE.
SLP =.FALSE.
MVE =.FALSE.
CNT =.FALSE.
C
READ(8,DAT,ERR= 202)
C
LOP = LOP + 1
C
IF ( OMG) GO TO 205
IF ( PSI) GO TO 205
IF ( CNT) GO TO 205
C
202 IDER = 1
WRITE(6,2200)
2200 FORMAT(/,' * Error in DATA',/)
WRITE (6,2300) IDT, OMG, PSI, SLP, MVE, CNT
2300 FORMAT(/,'*** IDT, OMG, PSI, SLP, MVE, CNT=',I2,7L2 )
GO TO 900
C
205 IF ( IDT.LE.0) GO TO 202
C ----------------------------------------------------------
DO 300 ICT = 1, IDT
C
ISLP = - 10
C ------------------------------------------------
DO 320 L = 1, 5
C
INDX( L) = 0
INDY( L) = 0
INDE( L) = 0
C
320 CONTINUE
C ------------------------------------------------
WVR = 0.D0
PVL = 0.D0
C
READ(8,215,ERR=225) ISLP, WVR, PVL,
1 ( INDX( L), INDY( L), INDE( L), L = 1, 5 )
C
215 FORMAT(I2,1X,2(F6.3,1X),5(3I3,1X))
C
IPT = 0
C
C IF ISLP = -1 Boundary : Downward.
C IF ISLP = 0 Boundary : Horizontal
C
C IF ISLP = 1 Boundary : Upward
C IF ISLP = 10 Boundary : Vertical
C
IF ( ISLP + 1 ) 225, 233, 221
221 IF ( ISLP ) 225, 232, 222
222 IF ( ISLP - 1 ) 225, 231, 223
223 IF ( ISLP - 10) 225, 230, 225
C
225 IDER = 1
WRITE (6,2400)
2400 FORMAT(/,' * Read Error or Error in slope spec.',/ )
WRITE (6,2500) ISLP, WVR, PVL, ( INDX( L),INDY( L),INDE( L),
1 L = 1, 5 )
2500 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 : LINE slopes Downward. IPT = 2 : Horiz.
C
C IPT = 3 : LINE slopes Upward. IPT = 4 : Vertical
C
C ------------------------------------------------
DO 340 L = 1, 5
C
IUR = INDE( L)
GO TO ( 240, 240, 240, 235), IPT
C
235 ILW = INDY( L)
GO TO 245
C
240 ILW = INDX( L)
245 I = INDX( L)
J = INDY( L)
C
C ----- X & Y coordinate : Assigned ---
C
C --------------------------------------
DO 290 IDX = ILW, IUR
C
IF ( I) 255, 340, 250
250 IF ( J) 255, 340, 251
251 IF ( I-IX) 252, 252, 255
252 IF ( J-IY) 260, 260, 255
C
C ----- I or J = 0 : No boundary specified ---
C
255 IDER = 1
WRITE(6,2600) L, IX, IY
2600 FORMAT(/,' Error in coordinate spec. number', I1,'. IX =',I3,
1 ' IY =',I3)
WRITE(6,2500) ISLP, WVR, PVL,
1 ( INDX( M), INDY( M), INDE( M), M = 1, 5 )
GO TO 340
C
260 IF ( .NOT.OMG ) GO TO 265
WPT( I,J,1) =.FALSE.
WPT( I,J,2) = SLP
WPT( I,J,3) = MVE
WVL( I,J) = WVR
C
265 IF ( .NOT.PSI) GO TO 270
C
270 IF ( .NOT.CNT ) GO TO 280
WPT( I,J,4) =.TRUE.
C
280 GO TO ( 286, 288, 287, 285 ), IPT
C
285 J = J + 1
GO TO 290
C
286 J = J - 1
GO TO 288
C
287 J = J + 1
288 I = I + 1
290 CONTINUE
C --------------------------------------
C
C ----- From ILW to INDE( L) ---
C
340 CONTINUE
C ------------------------------------------------
C
300 CONTINUE
C ----------------------------------------------------------
C
IF ( LOP.GE.LED) GO TO 310
GO TO 900
C
310 CONTINUE
C
IF ( IDER.EQ.0) RETURN
C
RETURN
END
C **********************************************************************
C
SUBROUTINE ISB1 ( WVL, S, WPT, IX, IY, XE, YE, CE, U, V,
1 NNP, XX, YY, WW )
C
IMPLICIT REAL*8(A-H,O-Z)
C
LOGICAL*1 WPT( IX,IY,4)
C
DIMENSION WVL( IX,IY), S ( IX,IY), U ( IX,IY), V ( IX,IY),
1 CE ( IX,IY), XE( IX,IY), YE( IX,IY), XX( NNP),
2 YY ( NNP), WW( NNP)
C
COMMON /CNST/ WMX, ERV, ERL, CMX, DDT, RLA, CMA, CMN, MSC
CHARACTER*6 MSC
COMMON /DAT2/ DX2, FS, CF, DY, DY2, TM, R, AL, BT, CFI
COMMON /DAT1/ DX, F, DT, NLT, NPR, UL, UU, D, VIS, WEI, ICR,
1 ISC
C
C ----------------------------
C
CALL TTLE
C
C ----------------------------
C
DDT = 0.D0
RLA = 0.D0
C --------------------------------------
DO 100 J = 1, IY
DO 100 I = 1, IX
C --------------------------------------
C
C ----------------------------
DO 150 K = 1, 4
C
WPT( I,J, K) = .FALSE.
150 CONTINUE
C ----------------------------
C
CE ( I,J) = 0.D0
XE ( I,J) = 0.D0
YE ( I,J) = 0.D0
c
U ( I,J) = 0.D0
V ( I,J) = 0.D0
C
WVL( I,J) = 0.D0
S ( I,J) = 0.D0
C
100 CONTINUE
C --------------------------------------
C
DO 200 I = 1, NNP
C
XX( I) = 0.D0
YY( I) = 0.D0
WW( I) = 0.D0
C
200 CONTINUE
C -----------------------------------------------------------
C
CALL INPD ( WVL, WPT, IX, IY )
C
C -----------------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE ISB2 ( U, V, XE, YE, CE, WVL, XX, YY,
1 WW, PSX, PSY, NDE, IX, IY, LOP, NNP,
2 NEL, II, JJ )
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION U ( IX,IY), V ( IX,IY), XE ( IX,IY),
1 YE ( IX,IY), CE ( IX,IY), WVL( IX,IY),
2 XX ( NNP), YY ( NNP), WW ( NNP),
3 PSX( 50000), PSY( 50000), NDE( NEL,4)
C
COMMON /CNST/ WMX, ERV, ERL, CMX, DDT, RLA, CMA, CMN, MSC
CHARACTER*6 MSC
C
C ----------------------------------------------------------------------
C
CALL SDAT ( U, V, XE, YE, WVL, XX, YY, WW, NDE, IX, IY, NNP, NEL )
C
C ----------------------------------------------------------------------
C
CALL ANSB ( XE, YE, CE, IX, IY, LOP )
C
C --------------------------------------------
DO 100 I = 1, IX
DO 100 J = 1, IY
C
WVL( I,J) = CE( I,J)
100 CONTINUE
C ---------------------------------------------------------------
C
CALL PSMV ( WVL, XE, YE, IX, IY, PSX, PSY, LOP, II, JJ )
C
C ---------------------------------------------------------------
C
C ----- Attention ---
C
CMX = 1.D0
C ===============
C
RETURN
END
C **********************************************************************
C
SUBROUTINE PLTS ( WVL, XE, YE, CE, XX, YY, WW, AWW, PSX,
1 PSY, IX, IY, LOP, II, JJ, NNP )
C
IMPLICIT REAL*8(A-H,O-Z)
C
COMMON /CNTL/ IFC
COMMON /CNST/ WMX, ERV, ERL, CMX, DDT, RLA, CMA, CMN, MSC
CHARACTER*6 MSC
COMMON /DAT1/ DX, F, DT, NLT, NPR, UL, UU, D, VIS, WEI, ICR,ISC
C
DIMENSION WVL( IX,IY), XE ( IX,IY), YE( IX,IY), CE ( IX,IY),
1 XX ( NNP), YY ( NNP), WW( NNP), AWW( NNP),
2 PSX( 50000), PSY( 50000)
C
C ----------------------------------------------------------------
C
CALL PSMV ( WVL, XE, YE, IX, IY, PSX, PSY, LOP, II, JJ )
C
C ----------------------------------------------------------------
C
CALL STMP ( WVL, WW, IX, IY, NNP)
C
C ----------------------------------------------------------------
C
CALL STMP ( CE, AWW, IX, IY, NNP)
C
C ----------------------------------------------------------------
WRITE(10,2000) IFC
2000 FORMAT(I2)
WRITE(10,2100) MSC, WEI, LOP, DDT, WMX
2100 FORMAT(A6,F8.4,I5,2F8.4)
WRITE(10,2200) ( XX( I), YY( I), WW( I), AWW( I), I = 1, NNP)
2200 FORMAT(15F8.4)
C
RETURN
END
C **********************************************************************
C
SUBROUTINE PSCP ( WVL, IX, IY, PSX, PSY, LOP )
C
C ----- At the Last Loop ---
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION WVL( IX,IY), SX ( 101), SY ( 101), SWX( 101),
1 SWY( 101), PSX( 50000), PSY( 50000)
C
COMMON /DAT1/ DX, F, DT, NLT, NPR, UL, UU, D, VIS, WEI, ICR, ISC
C
C --------------------------------------
WMX = 0.D0
C
DO 100 I = 1, IX
DO 100 J = 1, IY
C
IF ( WVL( I,J).GT.WMX) THEN
WMX = WVL( I,J)
II = I
JJ = J
END IF
C
100 CONTINUE
C ---------------------------------------
C
N = IDNINT( DFLOAT( IX+1)/ 2.D0)
C
C ---------------------------------------
DO 200 I = 1, IX
DO 200 J = 1, IY
C
C ----- A: DX & DY ---
C
SX ( I) = DX* DFLOAT( I-N)
SY ( J) = DX* DFLOAT( J-N)
SWX( I) = WVL( I, JJ)
SWY( J) = WVL( II, J )
C
200 CONTINUE
C -----------------------------------------------
C
CALL SPOS ( PSX, CMAX, SX, SWX, IX, LOP )
C
C -----------------------------------------------
C
CALL SPOS ( PSY, CMAY, SY, SWY, IY, LOP )
C
C -----------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE PSMV ( WVL, XE, YE, IX, IY, PSX, PSY, LOP, II, JJ )
C
C ----- Position of max. value ---
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION WVL( IX,IY), XE ( IX,IY), YE ( IX,IY),
1 PSX( 50000), PSY( 50000), SX ( 101),
2 SY ( 101), SWX( 101), SWY( 101)
C
COMMON /ECHR/ ERR( 50000), SCMA
COMMON /CNST/ WMX, ERV, ERL, CMX, DDT, RLA, CMA, CMN, MSC
CHARACTER*6 MSC
COMMON /DAT1/ DX, F, DT, NLT, NPR, UL, UU, D, VIS, WEI,
1 ICR, ISC

C
WMX = 0.D0
C ------------------------------------------------
DO 100 I = 1, IX
DO 100 J = 1, IY
C
IF ( DABS( WVL(I,J)).LE.WMX) GO TO 100
C -------------------------------------------
C
WMX = WVL( I,J)
II = I
JJ = J
C
100 CONTINUE
C ------------------------------------------------
C
IF ( LOP.EQ.0) GO TO 900
N = DFLOAT( IX+1)/ 2.D0
C
C ----- A : Mesh size in X - Direction ---
C
C ------------------------------------------------
DO 200 I = 1, IX
DO 200 J = 1, IY
C
SX ( I) = DX* DFLOAT( I-N)
SY ( J) = DX* DFLOAT( J-N)
SWX( I) = WVL( I,JJ )
SWY( J) = WVL( II, J )
C
200 CONTINUE
C ------------------------------------------------
C
CALL SPOS ( PSX, CMAX, SX, SWX, IX, LOP )
C
C ------------------------------------------------
C
CALL SPOS ( PSY, CMAY, SY, SWY, IY, LOP )
C
C ------------------------------------------------
900 CONTINUE
C
IF ( LOP.EQ.0)
1WRITE(6,2000) LOP, WMX, II, JJ, XE( II,JJ), YE( II,JJ)
2000 FORMAT(/,'* Loop=',I5,' Cmax=',F8.4,' (', I3,',', I3,')',
1 ' X =', F8.4,' Y =', F8.4,/)
IF ( LOP.EQ.0) SCMA = WMX
C
IF ( LOP.NE.0)
1WRITE(6,2100) LOP, ERV, WMX, II, JJ, PSX( LOP), PSY( LOP)
2100 FORMAT('* Loop=',I5,' Av_E=',F7.4,'(%)',' Cmax=',F8.4,
1 ' (', I3,',', I3,')',' X =',F8.4,' Y =',F8.4)
C
RETURN
END
C **********************************************************************
C
SUBROUTINE SBLT ( WVL, CE, XX, YY, WW, AWW, PSX, PSY,
1 DIF, IX, IY, LOP, NNP, KCP, * )
C
IMPLICIT REAL*8(A-H,O-Z)
C
COMMON /CNTL/ IFC
COMMON /ECHR/ ERR( 50000), SCMA
COMMON /CNST/ WMX, ERV, ERL, CMX, DDT, RLA, CMA, CMN, MSC
CHARACTER*6 MSC
COMMON /DAT1/ DX, F, DT, NLT, NPR, UL, UU, D, VIS, WEI,
1 ICR, ISC
C
DIMENSION WVL( IX,IY), CE ( IX,IY), XX ( NNP), YY ( NNP),
1 WW ( NNP), AWW( NNP), PSX( 50000), PSY( 50000),
2 DIF( IX,IY)
C
C -------------------------------------------
C
CALL STMP ( WVL, WW, IX, IY, NNP )
C
C -------------------------------------------
C
C ----- Anal. Sol. ----
C
C ------------------------------------------------
C
CALL STMP ( CE, AWW, IX, IY, NNP )
C
C ------------------------------------------------
C
CALL PSCP ( WVL, IX, IY, PSX, PSY, LOP )
C
C ------------------------------------------------
C
C ----- Error Computation ---
C
DSM = SCMA - CMA
WRITE(6,2000) LOP, ERL, CMA, CMN, DSM
2000 FORMAT(/,'* Last=',I5,' Err/Loop=',F8.4,'(%)',
1 ' Cmax=',F8.4,' Cmin=',F8.4,/,11X,' Diff_Cmax.=',F9.5,/)
C
IF ( CMA.GE.1.1D0) WRITE(6,*) ' ** Diverged & STOP * CMA >=1.1'
IF ( CMA.GE.1.1D0) STOP
C
IF ( ICR.EQ.0) GO TO 3333 ! NO Correction
C
RNMX = 0.D0
C ------------------------------------------------
DO 100 I = 1, IX
DO 100 J = 1, IY
C
IF ( DIF( I,J).LE.RNMX) GO TO 100
RNMX = DIF( I,J)
II = I
JJ = J
C
100 CONTINUE
C ------------------------------------------------
** AC_DIF = RNMX/ DT
AC_DIF = RNMX
C
WRITE(6,*) ' *** Correction ** Max. DIF( I,J)=', AC_DIF, II, JJ
C
3333 CONTINUE ! Attention : At the Last Loop
WRITE(10,2100) IFC
2100 FORMAT(I2)
WRITE(10,2200) MSC, WEI, LOP, DDT, WMX
2200 FORMAT(A6,F8.4,I5,2F8.4)
WRITE(10,2300) ( XX( I), YY( I), WW( I), AWW( I), I = 1, NNP)
2300 FORMAT(15F8.4)
C
WRITE(6,2400)
2400 FORMAT('*** To be Continued ? ** YES: "1" : ** NO: "0" **')
C
*** READ(5,*) KCP
KCP = 0
IF ( KCP.EQ.0) THEN
WRITE(7,2500) MSC, WEI, DT, LOP, DSM
2500 FORMAT(A6,2F8.4,I5,F9.5)
RETURN
END IF
C
REWIND ( 8)
RETURN 1
C
END
C **********************************************************************
C
SUBROUTINE SDAT ( U, V, XE, YE, WVL, XX, YY, WW,
1 NDE, IX, IY, NNP, NEL )
C
IMPLICIT REAL*8(A-H,O-Z)
C
COMMON /DAT1/ DX, F, DT, NLT, NPR, UL, UU, D, VIS, WEI, ICR,ISC
C
DIMENSION U ( IX,IY), V ( IX,IY), XE( IX,IY), YE( IX,IY),
1 WVL( IX,IY), XX( NNP), YY( NNP), WW( NNP),
2 NDE( NEL,4)
C
N = DFLOAT( IX+1)/ 2.D0
C
C ----- Velocity Field --- "DX " : DX or DY ---
C
DO 100 I = 1, IX
DO 100 J = 1, IY
C
XE( I,J) = DFLOAT( I-N )* DX
YE( I,J) = DFLOAT( J-N )* DX
U ( I,J) = - YE( I,J)
V ( I,J) = XE( I,J)
C
100 CONTINUE
C ---------------------------------------------
UMX = 0.D0
VMX = 0.D0
C
DO 200 I = 1, IX
DO 200 J = 1, IY
C
IF ( UMX.LE.U( I,J)) UMX = U( I,J)
IF ( VMX.LE.V( I,J)) VMX = V( I,J)
C
200 CONTINUE
C ----------------------------------------------
C
C ----- Check for Velocity Field ---
C
WRITE(6,2000) UMX, VMX
2000 FORMAT(' * Umax.=',F7.3,' Vmax.=',F7.3)
C
C -------------------------------------------------
C
CALL STMP ( XE, XX, IX, IY, NNP)
C
CALL STMP ( YE, YY, IX, IY, NNP)
C
CALL STMP ( WVL, WW, IX, IY, NNP)
C
CALL STND ( IX, IY, XX, YY, NDE, NNP, NEL )
C
C -------------------------------------------------
WRITE(10,2100) IX, IY, NNP, NEL
2100 FORMAT(4I5)
WRITE(10,2200) (( NDE( I,J),I = 1, NEL),J = 1, 4 )
2200 FORMAT(20I5)
C
RETURN
END
C **********************************************************************
C
SUBROUTINE SPOS ( PSX, CMA, X,C, NNP, LOP )
C
IMPLICIT REAL * 8(A-H,O-Z)
C
DIMENSION X( NNP), C( NNP), PSX( 50000)
C
PSX( LOP) = 0.D0
CMA = 0.D0
C
C ------------------------------------------------
DO 100 I = 1, NNP
C
IF ( CMA.LT.C( I)) PSX( LOP) = X( I)
IF ( CMA.LT.C( I)) CMA = C( I)
C
100 CONTINUE
C ------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE STMP ( DM1, DM2, IX, IY, NNP )
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION DM1( IX,IY), DM2( NNP)
C
NN = 0
C
C --------------------------------------
DO 100 I = 1, IX
DO 100 J = 1, IY
C
NN = NN + 1
DM2( NN) = DM1( I,J)
C
100 CONTINUE
C --------------------------------------
RETURN
END
C **********************************************************************
C
SUBROUTINE STND ( IX, IY, X, Y, NDE, NNP, NEL )
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION X( NNP), Y( NNP), NDE( NEL,4)
C
IF ( IX.EQ.IY) GO TO 90
C
STOP
C ==========
C
90 CONTINUE
C
DX = 1.D0/ DFLOAT( IX-1)
DY = 1.D0/ DFLOAT( IY-1)
C
NN = 0
C --------------------------------------------------
DO 100 I = 1, IX
DO 100 J = 1, IY
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 TTLE
C
IMPLICIT REAL*8(A-H,O-Z)
C
WRITE(6,*) ' '
WRITE(6,*) '*****************************************************'
WRITE(6,*) '* *'
WRITE(6,*) '* C V 2 D - F D M *'
WRITE(6,*) '* *'
WRITE(6,*) '* ( V.5.0 ) *'
WRITE(6,*) '* *'
WRITE(6,*) '* ( 4th & 2nd ordered FDM ) *'
WRITE(6,*) '* *'
WRITE(6,*) '* ( 2D-Convection Equation ) *'
WRITE(6,*) '* *'
WRITE(6,*) '* Copyright 2011 : Yasuhiro MATSUDA *'
WRITE(6,*) '* *'
WRITE(6,*) '*****************************************************'
WRITE(6,*) ' '
C
RETURN
END
C **********************************************************************
C
SUBROUTINE WDCP ( U, V, WVL, W, WPT, IX, IY, DIF )
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION WVL( IX,IY), W( IX,IY), U( IX,IY), V( IX,IY),
1 DIF( IX,IY)
C
LOGICAL*1 WPT( IX,IY,4)
C
COMMON /DAT2/ DX2, FS, CF, DY, DY2, TM, R, AL, BT, CFI
COMMON /DAT1/ DX, F, DT, NLT, NPR, UL, UU, D, VIS, WEI,ICR,
1 ISC
C
C ----- With Correction ----- Attention in case of DX = DY ---
C
C ----- Cotrrected Diffusivity -----------------------------------------
C
DO 100 I = 1, IX
DO 100 J = 1, IY
C
DIF( I,J) = VIS + 0.5D0* U( I,J)* V( I,J)* DT
C
100 CONTINUE
C ----------------------------------------------------------------------
DO 200 J = 1, IY
DO 200 I = 1, IX
C
IF ( WPT( I,J,1)) GO TO 200
WVL( I,J) = W( I,J)
C
200 CONTINUE
C ---------------------------------------------------------------
DO 300 J = 1, IY
DO 300 I = 1, IX
C
IF ( WPT( I,J,1)) GO TO 1 ! Fluid Point
IF ( WPT( I,J,2)) GO TO 333 ! Slip BC
IF ( WPT( I,J,3)) GO TO 300 ! Moving BC
C
C ----- Rigid, No-SLP, Not moving ---
C
GO TO 333
C
1 CONTINUE
C
C ----- Normal Interior to fluid ---
C
WVL( I,J) = WVL( I,J)
C
C ----- DIF(i,J) : Diffusion Coefficient with Correction ---
C
1 + DIF( I,J)* DT* (( W( I,J+1) + W( I,J-1) - 2.D0* W( I,J))/ FS
2 + ( W( I+1,J) + W( I-1,J) - 2.D0* W( I,J))/ DX2)
C
C ----- FS = 1 for DX = DY ---
C
GO TO 300
C
333 CONTINUE
C
300 CONTINUE
C ---------------------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE WU4C ( WVL, WS, WPT, IX, IY, UF )
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION WVL( IX,IY), WS( IX,IY), UF( IX,IY)
C
LOGICAL *1 WPT( IX,IY,4)
C
COMMON /DAT2/ DX2, FS, CF, DY, DY2, TM, R, AL, BT, CFI
COMMON /DAT1/ DX, F, DT, NLT, NPR, UL, UU, D, VIS, WEI, ICR,
1 ISC
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
DATA A11 / 0.5833333333 333333 D0/, ! 7/12
1 A12 / 0.625 D0/, ! 5/8
2 A13 / - 0.0833333333 333333 D0/, ! -1/12
3 A14 / - 0.125 D0/, ! -1/8
C
4 A21 / - 0.0833333333 333333 D0/, ! -1/12
5 A22 / - 0.0416666666 666667 D0/, ! -1/24
6 A23 / 0.0833333333 333333 D0/, ! 1/12
7 A24 / 0.0416666666 666667 D0/ ! 1/24
C
C --------------------------------------------------------------------
DO 100 I = 1, IX
DO 100 J = 1, IY
C --------------------------------------------------------------------
C
IF ( I.EQ.IX) GO TO 100
C
ALF = UF( I,J)* AL ! UF* DT/DX
C =======================
C
IF ( ALF.LT.0.D0) GO TO 311 ! U < 0
C ==========================================
C
IF ( .NOT. WPT( I,J,1)) GO TO 600 ! Boundary
C
IF ( WPT( I,J,4)) GO TO 600 ! Cont.
IF ( WPT( I+1,J,4)) GO TO 600 ! Cont.
IF ( WPT( I-1,J,4)) GO TO 600 ! Cont.
IF ( WPT( I+1,J,4)) GO TO 600 ! Cont.
C
IF ( I.EQ.1) GO TO 400
C
IF ( .NOT.WPT( I-1,J,1)) GO TO 200 ! Boundary
IF ( WPT( I-1,J,4)) GO TO 200 ! Cont.
IF ( .NOT.WPT( I+1,J,1)) GO TO 200 ! Boundary
IF ( WPT( I+1,J,4)) GO TO 200 ! Cont.
C
C ------------------------------------------------
400 CONTINUE
C
AA = ALF* ALF ! B** 2
AAA = ALF* AA ! B** 3
AAAA= ALF* AAA ! B** 4
C
A1 = A11* ALF ! A11* B
A2 = A12* AA ! A12* B** 2
A3 = A13* AAA ! A13* B** 3
A4 = A14* AAAA ! A14* B** 4
C
B1 = A21* ALF ! A21* B
B2 = A22* AA ! A22* B** 2
B3 = A23* AAA ! A23* B** 3
B4 = A24* AAAA ! A24* B** 4
C
C --------------------------------------
C
B = ALF - 1.D00
C ====================
C
BB = B* B
BBB = B* BB
BBBB= B* BBB
C
* C1 = A11* B
* C2 = A12* BB
* C3 = A13* BBB
* C4 = A14* BBBB
C
C1 = A11* B + A11
C2 = A12* BB - A12
C3 = A13* BBB + A13
C4 = A14* BBBB - A14
C
D1 = A21* B
D2 = A22* BB
D3 = A23* BBB
D4 = A24* BBBB
C
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 = ( WS( I-2,J) + WS( I+1,J))* ( 1.D0 - WEI)
FDM = ( WS( I-2,J) - WS( I+1,J))* ( 1.D0 - WEI)
C
308 FBP = ( WS( I-1,J) + WS( I+2,J))* WEI
FBM = ( WS( I-1,J) - WS( I+2,J))* WEI
C
307 FCP = ( WS( I-1,J) + WS( I,J))* ( 1.D0 - WEI)
FCM = ( WS( I-1,J) - WS( I,J))* ( 1.D0 - WEI)
C
306 FAP = ( WS( I,J) + WS( I+1,J))* WEI
FAM = ( WS( I,J) - WS( I+1,J))* WEI
C
TF = A1* FAP + A2* FAM + A3* FAP + A4* FAM
1 + B1* FBP + B2* FBM + B3* FBP + B4* FBM
C
2 + C1* FCP + C2* FCM + C3* FCP + C4* FCM
3 + D1* FDP + D2* FDM + D3* FDP + D4* FDM
C
WVL( I,J) = WVL( I,J) - TF
WVL( I+1,J) = WVL( I+1,J) + TF
C
IF ( I+1.LT.IX) GO TO 302
IF ( .NOT. WPT( I+1,J,1)) GO TO 100 ! Boundary
C
C ----- Periodic Right ---
C
WVL( I+1-(IX-1),J) = WVL( I+1-(IX-1),J) + TF
GO TO 100
C
302 CONTINUE
IF ( I.GT.1) GO TO 100
IF ( .NOT. WPT( I,J,1)) GO TO 100 ! Boundary
C
C ---- Periodic Left ---
C
WVL( I+(IX-1),J) = WVL( I+(IX-1),J) - TF
GO TO 100
C
304 CONTINUE
FDP = ( WS( I-2 + ( IX-1),J) + WS( I+1,J))* ( 1.D0 - WEI)
FDM = ( WS( I-2 + ( IX-1),J) - WS( I+1,J))* ( 1.D0 - WEI)
FBP = ( WS( I-1 + ( IX-1),J) + WS( I+2,J))* WEI
C
FBM = ( WS( I-1 + ( IX-1),J) - WS( I+2,J))* WEI
FCP = ( WS( I-1 + ( IX-1),J) + WS( I,J)) * ( 1.D0 - WEI)
FCM = ( WS( I-1 + ( IX-1),J) - WS( I,J)) * ( 1.D0 - WEI)
GO TO 306
C
305 CONTINUE
FDP = ( WS( I-2 + (IX-1),J) + WS(I+1,J))* ( 1.D0 - WEI)
FDM = ( WS( I-2 + (IX-1),J) - WS(I+1,J))* ( 1.D0 - WEI)
GO TO 308
C
309 CONTINUE
FBP = ( WS( I-1,J) + WS( I+2-( IX-1),J))* WEI
FBM = ( WS( I-1,J) - WS( I+2-( IX-1),J))* WEI
FDP = ( WS( I-2,J) + WS( I+1,J))* ( 1.D0 - WEI)
FDM = ( WS( I-2,J) - WS( I+1,J))* ( 1.D0 - WEI)
GO TO 307
C
C ----------------------------------------------------------
C
311 CONTINUE
IF ( .NOT. WPT( I+1,J,1)) GO TO 601 ! Boundary
IF ( WPT( I+1,J,4)) GO TO 601 ! Cont.
IF ( WPT( I,J,4)) GO TO 601 ! Cont.
C
K = I + 1
IF ( WPT( K+1,J,4)) GO TO 201
IF ( WPT( K-1,J,4)) GO TO 201
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
B = ALF + 1.D0
BB = B* B
BBB = B* BB
BBBB= B* BBB
C
C1 = A11* B - A11
C2 = A12* BB - A12
C3 = A13* BBB - A13
C4 = A14* BBBB - A14
C
D1 = A21* B
D2 = A22* BB
D3 = A23* BBB
D4 = A24* BBBB
C
IF ( K.EQ.IX) GO TO 314
IF ( K+1 .EQ.IX) GO TO 315
IF ( K-1 .EQ.1 ) GO TO 319
FDP = ( WS( K-1,J) + WS( K+2,J))* ( 1.D0 - WEI)
FDM = ( WS( K-1,J) - WS( K+2,J))* ( 1.D0 - WEI)
C
318 CONTINUE
FBP = ( WS( K-2,J) + WS( K+1,J))* WEI
FBM = ( WS( K-2,J) - WS( K+1,J))* WEI
C
317 CONTINUE
FCP = ( WS( K,J) + WS( K+1,J))* ( 1.D0 - WEI)
FCM = ( WS( K,J) - WS( K+1,J))* ( 1.D0 - WEI)
C
316 CONTINUE
FAP = ( WS( K-1,J) + WS( K,J))* WEI
FAM = ( WS( K-1,J) - WS( 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
WVL( K-1,J) = WVL( K-1,J) - TF
WVL( K,J) = WVL( 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
WVL( K -(IX-1),J) = WVL( K -(IX-1),J) + TF
GO TO 100
C
312 IF ( K-1.GT.1) GO TO 100
C
IF ( .NOT. WPT( K-1,J,1)) GO TO 100
C
C ----- Periodic Left ---
C
WVL( K-1+(IX-1),J) = WVL( K-1+(IX-1),J) - TF
C
GO TO 100
C
314 CONTINUE
C
FDP = ( WS( K-1,J) + WS( K+2-(IX-1),J))* ( 1.D0 - WEI)
FDM = ( WS( K-1,J) - WS( K+2-(IX-1),J))* ( 1.D0 - WEI)
C
FBP = ( WS( K-2,J) + WS( K+1-(IX-1),J))* WEI
FBM = ( WS( K-2,J) - WS( K+1-(IX-1),J))* WEI
C
FCP = ( WS( K,J ) + WS( K+1-(IX-1),J))* ( 1.D0 - WEI)
FCM = ( WS( K,J ) - WS( K+1-(IX-1),J))* ( 1.D0 - WEI)
C
GO TO 316
C
315 CONTINUE
C
FDP = ( WS( K-1,J) + WS( K+2-( IX-1),J))* ( 1.D0 - WEI)
FDM = ( WS( K-1,J) - WS( K+2-( IX-1),J))* ( 1.D0 - WEI)
C
GO TO 318
C
319 CONTINUE
C
FBP = ( WS( K-2 + ( IX-1),J) + WS( K+1,J))* WEI
FBM = ( WS( K-2 + ( IX-1),J) - WS( K+1,J))* WEI
C
FDP = ( WS( K-1,J) + WS( K+2,J))* ( 1.D0 - WEI)
FDM = ( WS( K-1,J) - WS( K+2,J))* ( 1.D0 - WEI)
C
GO TO 317
C
600 TF = ALF* WS( I,J)
C
GO TO 602
C
601 TF = ALF* WS( I+1,J)
C
602 WVL( I,J) = WVL( I,J) - TF
WVL( I+1,J) = WVL( I+1,J) + TF
C
GO TO 100
C
200 CONTINUE
A1 = 0.5D0* ALF
A2 = A1* ALF
B2 = A2 - 2.D0* A1
C
TRP = ( WS( I-1,J) + WS( I,J))* ( 1.D0 - WEI)
TRM = ( WS( I-1,J) - WS( I,J))* ( 1.D0 - WEI)
TPL = ( WS( I,J) + WS( I+1,J))* WEI
TML = ( WS( I,J) - WS( I+1,J))* WEI
TF = A1* TPL + A1* TRP + A2* TML + B2* TRM
C
WVL( I+1,J) = WVL( I+1,J) + TF
WVL( I,J) = WVL( I,J) - TF
C
GO TO 100
C
201 A1 = 0.5D0* ALF
A2 = A1* ALF
B2 = A2 + 2.D0* A1
C
TRP = ( WS( K,J) + WS( K+1,J))* ( 1.D0 - WEI)
TRM = ( WS( K,J) - WS( K+1,J))* ( 1.D0 - WEI)
TPL = ( WS( K-1,J) + WS( K,J)) * WEI
TML = ( WS( K-1,J) - WS( K,J)) * WEI
TF = A1* TPL + A1* TRP + A2* TML + B2* TRM
C
WVL( K,J) = WVL( K,J) + TF
WVL( K-1,J) = WVL( K-1,J) - TF
C
100 CONTINUE
C -------------------------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE WU2C ( WVL, WS, WPT, IX, IY, UF )
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION WVL( IX,IY), WS( IX,IY), UF( IX,IY)
C
LOGICAL*1 WPT( IX,IY,4)
C
COMMON /DAT2/ DX2, FS, CF, DY, DY2, TM, R, AL, BT, CFI
COMMON /DAT1/ DX, F, DT, NLT, NPR, UL, UU, D, VIS, WEI, ICR,
1 ISC
C
C ----------------------------------------------------------
DO 100 I = 1, IX
DO 100 J = 1, IY
C ----------------------------------------------------------
C
IF ( I.EQ.IX) GO TO 100
C
ALF = UF( I,J)* AL
IF ( ALF.LT.0.D0) GO TO 201
IF ( .NOT. WPT( I,J,1)) GO TO 150 ! Boundary
IF ( WPT( I,J,4)) GO TO 150 ! Cont.
IF ( WPT( I+1,J,4)) GO TO 150 ! Cont.
C
A1 = 0.5D0* ALF
A2 = A1* ALF
B2 = A2 - 2.D0* A1
C
TRP = ( WS( I-1,J) + WS( I,J))* ( 1.D0 - WEI)
TRM = ( WS( I-1,J) - WS( I,J))* ( 1.D0 - WEI)
TPL = ( WS( I,J) + WS( I+1,J))* WEI
TML = ( WS( I,J) - WS( I+1,J))* WEI
TF = A1* TPL + A1* TRP + A2* TML + B2* TRM
C
WVL( I+1,J) = WVL( I+1,J) + TF
WVL( I,J) = WVL( I,J) - TF
C
GO TO 100
C
201 CONTINUE
C
IF ( .NOT. WPT( I+1,J,1)) GO TO 170 ! Boundary
IF ( WPT( I+1,J,4)) GO TO 170 ! Cont.
IF ( WPT( I,J,4)) GO TO 170 ! Cont.
C
K = I + 1
A1 = 0.5D0* ALF
A2 = A1* ALF
B2 = A2 + 2.D0* A1
C
TRP = ( WS( K,J) + WS( K+1,J))* ( 1.D0 - WEI)
TRM = ( WS( K,J) - WS( K+1,J))* ( 1.D0 - WEI)
TPL = ( WS( K-1,J) + WS( K,J)) * WEI
TML = ( WS( K-1,J) - WS( K,J)) * WEI
TF = A1* TPL + A1* TRP + A2* TML + B2* TRM
C
WVL( K,J) = WVL( K,J) + TF
WVL( K-1,J) = WVL( K-1,J) - TF
C
GO TO 100
C
150 TF = ALF* WS( I,J)
C
GO TO 190
C
170 TF = ALF* WS( I+1,J)
C
190 WVL( I,J) = WVL( I,J) - TF
WVL( I+1,J) = WVL( I+1,J) + TF
C
100 CONTINUE
C ----------------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE WV4C ( WVL, WS, WPT, IX, IY, VF )
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION WVL( IX,IY), WS( IX,IY), VF( IX,IY)
C
LOGICAL *1 WPT( IX,IY, 4)
C
COMMON /DAT2/ DX2, FS, CF, DY, DY2, TM, R, AL, BT, CFI
COMMON /DAT1/ DX, F, DT, NLT, NPR, UL, UU, D, VIS, WEI,ICR,
1 ISC
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
DATA A11 / 0.5833333333 333333 D0/, ! 7/12
1 A12 / 0.625 D0/, ! 5/8
2 A13 / - 0.0833333333 333333 D0/, ! -1/12
3 A14 / - 0.125 D0/, ! -1/8
C
4 A21 / - 0.0833333333 333333 D0/, ! -1/12
5 A22 / - 0.0416666666 666667 D0/, ! -1/24
6 A23 / 0.0833333333 333333 D0/, ! 1/12
7 A24 / 0.0416666666 666667 D0/ ! 1/24
C
C ----------------------------------------------------------
DO 100 I = 1, IX
DO 100 J = 1, IY
C ----------------------------------------------------------
C
IF ( J.EQ.IY) GO TO 100
C
BET = VF( I,J)* BT
C
IF (BET .LT. 0.D0) GO TO 311
IF ( .NOT.WPT( I,J,1)) GO TO 600 ! Boundary
IF ( WPT( I,J,4)) GO TO 600 ! Cont.
IF ( WPT( I,J+1, 4)) GO TO 600 ! Cont.
IF ( WPT( I,J-1, 4)) GO TO 200 ! Cont.
IF ( WPT( I,J+1, 4)) GO TO 200 ! Cont.
C
AA = BET* BET
AAA = BET* AA
AAAA = BET* AAA
A1 = A11* BET
A2 = A12* AA
A3 = A13* AAA
A4 = A14* AAAA
C
B1 = A21* BET
B2 = A22* AA
B3 = A23* AAA
B4 = A24* AAAA
B = BET - 1.D0
BB = B* B
BBB = B* BB
BBBB = B* BBB
C
C1 = A11* B + A11
C2 = A12* BB - A12
C3 = A13* BBB + A13
C4 = A14* BBBB - A14
C
D1 = A21* B
D2 = A22* BB
D3 = A23* BBB
D4 = A24* BBBB
C
IF ( J .EQ. 1) GO TO 304
IF ( J-1.EQ. 1) GO TO 305
IF ( J+1.EQ.IY) GO TO 309
C
FDP = ( WS( I,J-2 ) + WS( I,J+1))* ( 1.D0 - WEI)
FDM = ( WS( I,J-2 ) - WS( I,J+1))* ( 1.D0 - WEI)
C
308 FBP = ( WS( I,J-1) + WS( I,J + 2 ))* WEI
FBM = ( WS( I,J-1) - WS( I,J + 2 ))* WEI
C
307 FCP = ( WS( I,J-1) + WS( I,J))* ( 1.D0 - WEI)
FCM = ( WS( I,J-1) - WS( I,J))* ( 1.D0 - WEI)
C
306 FAP = ( WS( I,J) + WS( I,J+1))* WEI
FAM = ( WS( I,J) - WS( I,J+1))* WEI
C
TF = A1* FAP + A2* FAM + A3* FAP + A4* FAM
1 + B1* FBP + B2* FBM + B3* FBP + B4* FBM
C
2 + C1* FCP + C2* FCM + C3* FCP + C4* FCM
3 + D1* FDP + D2* FDM + D3* FDP + D4* FDM
C
WVL( I,J) = WVL( I,J) - TF
WVL( I,J+1) = WVL( I,J+1) + TF
C
IF ( J+1.LT.IY) GO TO 302
IF ( .NOT. WPT( I,J+1,1)) GO TO 100
C
C ----- Periodic Right ---
C
WVL( I,J+1-(IY-1)) = WVL( I,J+1-(IY-1)) + TF
GO TO 100
C
302 CONTINUE
C
IF ( J .GT.1) GO TO 100
IF ( .NOT. WPT( I,J,1)) GO TO 100
C
C ----- Periodic Left ---
C
WVL( I,J+(IY-1)) = WVL( I,J +(IY-1)) - TF
GO TO 100
C
304 FDP = ( WS( I,J-2 + ( IY-1)) + WS( I,J+1))* ( 1.D0 - WEI)
FDM = ( WS( I,J-2 + ( IY-1)) - WS( I,J+1))* ( 1.D0 - WEI)
FBP = ( WS( I,J-1 + ( IY-1)) + WS( I,J + 2 ))* WEI
C
FBM = ( WS( I,J-1 + ( IY-1)) - WS( I,J + 2 ))* WEI
FCP = ( WS( I,J-1 + ( IY-1)) + WS( I,J))* ( 1.D0 - WEI)
FCM = ( WS( I,J-1 + ( IY-1)) - WS( I,J))* ( 1.D0 - WEI)
GO TO 306
C
305 FDP = ( WS( I,J-2 + ( IY-1)) + WS( I,J+1))* ( 1.D0 - WEI)
FDM = ( WS( I,J-2 + ( IY-1)) - WS( I,J+1))* ( 1.D0 - WEI)
GO TO 308
C
309 CONTINUE
FBP = ( WS( I,J-1) + WS( I,J + 2-( IY-1)))* WEI
FBM = ( WS( I,J-1) - WS( I,J + 2-( IY-1)))* WEI
FDP = ( WS( I,J-2) + WS( I,J+1)) * ( 1.D0 - WEI)
FDM = ( WS( I,J-2) - WS( I,J+1)) * ( 1.D0 - WEI)
GO TO 307
C
311 CONTINUE
IF ( .NOT. WPT( I,J+1,1)) GO TO 601
IF ( WPT( I,J+1,4)) GO TO 601
IF ( WPT( I,J,4)) GO TO 601
C
K = J + 1
IF ( WPT( I,K+1,4)) GO TO 201
IF ( WPT( I,K-1,4)) GO TO 201
C
AA = BET* BET
AAA = BET* AA
AAAA = BET* AAA
C
A1 = A11* BET
A2 = A12* AA
A3 = A13* AAA
A4 = A14* AAAA
C
B1 = A21* BET
B2 = A22* AA
B3 = A23* AAA
B4 = A24* AAAA
C
B = BET + 1.D0
BB = B* B
BBB = B* BB
BBBB = B* BBB
C
C1 = A11* B - A11
C2 = A12* BB - A12
C3 = A13* BBB - A13
C4 = A14* BBBB - A14
C
D1 = A21* B
D2 = A22* BB
D3 = A23* BBB
D4 = A24* BBBB
C
IF ( K.EQ. IY) GO TO 314
IF ( K+1 .EQ.IY) GO TO 315
IF ( K-1 .EQ.1 ) GO TO 319
C
FDP = ( WS( I,K-1) + WS( I,K+2))* ( 1.D0 - WEI)
FDM = ( WS( I,K-1) - WS( I,K+2))* ( 1.D0 - WEI)
C
318 FBP = ( WS( I,K-2) + WS( I,K+1))* WEI
FBM = ( WS( I,K-2) - WS( I,K+1))* WEI
C
317 FCP = ( WS( I,K ) + WS( I,K+1))* ( 1.D0 - WEI)
FCM = ( WS( I,K ) - WS( I,K+1))* ( 1.D0 - WEI)
C
316 FAP = ( WS( I,K-1) + WS( I,K ))* WEI
FAM = ( WS( I,K-1) - WS( I,K ))* WEI
C
TF = A1* FAP + A2* FAM + A3* FAP + A4* FAM
1 + B1* FBP + B2* FBM + B3* FBP + B4* FBM
C
2 + C1* FCP + C2* FCM + C3* FCP + C4* FCM
3 + D1* FDP + D2* FDM + D3* FDP + D4* FDM
C
WVL( I,K-1) = WVL( I,K-1) - TF
WVL( I,K ) = WVL( I,K) + TF
C
IF ( K.LT.IY) GO TO 312
IF ( .NOT.WPT( I,K,1)) GO TO 100
C
C ----- Periodic Right ---
C
WVL( I,K-(IY-1)) = WVL( I,K-(IY-1)) + TF
C
GO TO 100
C
312 CONTINUE
C
IF ( K-1.GT.1) GO TO 100
IF ( .NOT. WPT( I,K-1,1)) GO TO 100
C
C ----- Periodic Left ---
C
WVL( I,K-1+(IY-1)) = WVL( I,K-1+(IY-1)) - TF
GO TO 100
C
314 CONTINUE
FDP = ( WS( I,K-1) + WS( I,K+2-( IY-1)))* ( 1.D0 - WEI)
FDM = ( WS( I,K-1) - WS( I,K+2-( IY-1)))* ( 1.D0 - WEI)
C
FBP = ( WS( I,K-2) + WS( I,K+1-( IY-1)))* WEI
FBM = ( WS( I,K-2) - WS( I,K+1-( IY-1)))* WEI
C
FCP = ( WS( I,K) + WS( I,K+1-( IY-1)))* ( 1.D0 - WEI)
FCM = ( WS( I,K) - WS( I,K+1-( IY-1)))* ( 1.D0 - WEI)
GO TO 316
C
315 CONTINUE
C
FDP = ( WS( I,K-1) + WS( I,K+2-( IY-1)))* ( 1.D0 - WEI)
FDM = ( WS( I,K-1) - WS( I,K+2-( IY-1)))* ( 1.D0 - WEI)
GO TO 318
C
319 CONTINUE
C
FBP = ( WS( I,K-2 + ( IY-1)) + WS( I,K+1))* WEI
FBM = ( WS( I,K-2 + ( IY-1)) - WS( I,K+1))* WEI
FDP = ( WS( I,K-1) + WS( I,K+2 ))* ( 1.D0 - WEI)
FDM = ( WS( I,K-1) - WS( I,K+2 ))* ( 1.D0 - WEI)
C
GO TO 317
C
600 TF = BET* WS( I,J)
C
GO TO 602
C
601 TF = BET* WS( I,J+1)
C
602 WVL( I,J) = WVL( I,J) - TF
WVL( I,J+1) = WVL( I,J+1) + TF
C
GO TO 100
C
200 A1 = 0.5D0* BET
A2 = A1* BET
B2 = A2 - 2.D0* A1
C
TRP = ( WS( I,J-1) + WS( I,J)) * ( 1.D0 - WEI)
TRM = ( WS( I,J-1) - WS( I,J)) * ( 1.D0 - WEI)
TPL = ( WS( I,J) + WS( I,J+1))* WEI
TML = ( WS( I,J) - WS( I,J+1))* WEI
C
TF = A1* TPL + A1* TRP + A2* TML + B2* TRM
C
WVL( I,J+1) = WVL( I,J+1) + TF
WVL( I,J) = WVL( I,J) - TF
C
GO TO 100
C
201 A1 = 0.5D0* BET
A2 = A1* BET
B2 = A2 + 2.D0* A1
C
TRP = ( WS( I,K) + WS( I,K+1))* ( 1.D0 - WEI)
TRM = ( WS( I,K) - WS( I,K+1))* ( 1.D0 - WEI)
TPL = ( WS( I,K-1) + WS( I,K)) * WEI
TML = ( WS( I,K-1) - WS( I,K)) * WEI
TF = A1* TPL + A1* TRP + A2* TML + B2* TRM
C
WVL( I,K) = WVL( I,K) + TF
WVL( I,K-1) = WVL( I,K-1) - TF
C
100 CONTINUE
C ----------------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE WV2C ( WVL, WS, WPT, IX, IY, VF )
C
IMPLICIT REAL * 8(A-H,O-Z)
C
DIMENSION WVL( IX,IY), WS( IX,IY), VF( IX,IY)
C
LOGICAL*1 WPT( IX,IY,4)
C
COMMON /DAT2/ DX2, FS, CF, DY, DY2, TM, R, AL, BT, CFI
COMMON /DAT1/ DX, F, DT, NLT, NPR, UL, UU, D, VIS, WEI,ICR,
1 ISC
C
C ----------------------------------------------------------
DO 100 I = 1, IX
DO 100 J = 1, IY
C ----------------------------------------------------------
C
IF ( J.EQ.IY) GO TO 100
C
BET = VF( I,J)* BT
C
IF ( BET.LT.0.D0) GO TO 201
IF ( .NOT. WPT( I,J,1)) GO TO 150 ! Boundary
IF ( WPT( I,J,4)) GO TO 150 ! Cont.
IF ( WPT( I,J+1,4)) GO TO 150 ! Cont.
C
A1 = 0.5D0* BET
A2 = A1* BET
B2 = A2 - 2.D0* A1
C
TRP = ( WS( I,J-1) + WS( I,J))* ( 1.D0 - WEI)
TRM = ( WS( I,J-1) - WS( I,J))* ( 1.D0 - WEI)
TPL = ( WS( I,J) + WS( I,J+1))* WEI
TML = ( WS( I,J) - WS( I,J+1))* WEI
C
TF = A1* TPL + A1* TRP + A2* TML + B2* TRM
C
WVL( I,J+1) = WVL( I,J+1) + TF
WVL( I,J) = WVL( I,J) - TF
C
GO TO 100
C
201 CONTINUE
IF ( .NOT. WPT( I,J+1,1)) GO TO 170
IF ( WPT( I,J+1,4)) GO TO 170
IF ( WPT( I,J, 4)) GO TO 170
C
K = J + 1
A1 = 0.5D0* BET
A2 = A1* BET
B2 = A2 + 2.D0* A1
C
TRP = ( WS( I,K) + WS( I,K+1))* ( 1.D0 - WEI)
TRM = ( WS( I,K) - WS( I,K+1))* ( 1.D0 - WEI)
TPL = ( WS( I,K-1) + WS( I,K)) * WEI
TML = ( WS( I,K-1) - WS( I,K)) * WEI
C
TF = A1* TPL + A1* TRP + A2* TML + B2* TRM
C
WVL( I,K) = WVL( I,K) + TF
WVL( I,K-1) = WVL( I,K-1) - TF
C
GO TO 100
C
150 TF = BET* WS( I,J)
C
GO TO 602
C
170 TF = BET* WS( I,J+1)
602 WVL( I,J) = WVL( I,J) - TF
WVL( I,J+1) = WVL( I,J+1) + TF
C
100 CONTINUE
C ----------------------------------------------------------
C
RETURN
END
C ****************************************************************