戻る

C ********************************************************* C
C C
C F 2 D - FDM ( SV ) C
C C
C 2D Fluid Analysis by the 4th order FDM C
C C
C ( Version 3.5 ) C
C C
C Cavity Flow Problem C
C C
C Copyright : Yasuhiro MATSUDA C
C C
C ( Meshes of Equal length & Unequal length ) C
C C
C ********************************************************* C
C
PROGRAM MAIN
C
IMPLICIT REAL*8 (A-H,O-Z)
C
INTEGER X, Y
C
C ----- DATA ( 1 ) ---
C
PARAMETER ( X = 21, Y = 21, IX1 = X-1, IY1 = Y-1,
1 IXY = X* Y, IXY1= IX1* IY1)
C
*** PARAMETER ( X= 41, Y= 41, IX1=X-1, IY1=Y-1, IXY=X*Y, IXY1=IX1*IY1 )
*** PARAMETER ( X= 51, Y= 51, IX1=X-1, IY1=Y-1, IXY=X*Y, IXY1=IX1*IY1)
*** PARAMETER ( X= 61, Y= 61, IX1=X-1, IY1=Y-1, IXY=X*Y, IXY1=IX1*IY1 )
*** PARAMETER ( X= 81, Y= 81, IX1=X-1, IY1=Y-1, IXY=X*Y, IXY1=IX1*IY1 )
*** PARAMETER ( X=101, Y=101, IX1=X-1, IY1=Y-1, IXY=X*Y, IXY1=IX1*IY1 )
C
C ---------------------------------------
C
C MSH = 1 : Mesh of equal length
C MSH = 2 : Mesh of unequal length
C
C IHS = 1 : With Correction
C WEI : Weighting parameter
C
C ---------------------------------------
C
LOGICAL*1 WPT( X,Y,4), PPT( X,Y,4), EPT( X,Y,4)
C
DIMENSION PS1( X,Y), PS ( X,Y), PSM( X,Y), WP1( X,Y),
1 FF ( X,Y), WD ( X,Y), WZ ( X,Y), U ( X,Y),
2 V ( X,Y), UB( X,Y), VB ( X,Y), HH ( X,Y)
C
DIMENSION RDP( X,Y), BX( X,Y), BY( X,Y), XZ ( X),
1 YZ ( Y), HX( IX1), HY( IY1)
C
DIMENSION UF1( Y), UF2( Y), UF3( Y),
1 UNA( X,4), UNB( X,4), UNC( X,4), UND( X,4), UNE( X,4),
2 UAP( X,4), UBP( X,4), UCP( X,4), UDP( X,4), UEP( X,4),
3 UAM( X,4), UBM( X,4), UCM( X,4), UDM( X,4), UEM( X,4)
C
DIMENSION VF1( X), VF2( X), VF3( X),
1 VNA( Y,4), VNB( Y,4), VNC( Y,4), VND( Y,4), VNE( Y,4),
2 VAP( Y,4), VBP( Y,4), VCP( Y,4), VDP( Y,4), VEP( Y,4),
3 VAM( Y,4), VBM( Y,4), VCM( Y,4), VDM( Y,4), VEM( Y,4)
C
DIMENSION RAF( IX1), RAS( IY1), PE( IXY), XX( IXY), YY( IXY),
1 NDE( IXY1,4)
C
COMMON /CTRL/ MSH, IHS, WEI
COMMON /IDAT/ DT, LOPE, REN, LOP, NPL
C
C ----- Attention --- X, Y < = 51 x 51 ---
C
COMMON /CM1/ PA1( 51,51), PA2( 51,51), PA3( 51,51),
1 PB1( 51,51), PB2( 51,51), PB3( 51,51)
C
COMMON /CM2/ TM1( 51), TM2( 51), TM3( 51), TY1( 51), TY2( 51),
1 TY3( 51)
C
C ------------------------------------------------
C
*** OPEN ( 7, FILE='F2D_FDM_R100_UM20_FIG.DAT' )
*** OPEN ( 8, FILE='F2D_FDM_R100_UM20_OUT.DAT' )
*** OPEN ( 10, FILE='F2D_FDM_20.DAT' )
C
*** OPEN ( 7, FILE='F2D_FDM_R1000_UM50_FIG.DAT' )
*** OPEN ( 8, FILE='F2D_FDM_R1000_UM50_OUT.DAT' )
*** OPEN ( 10, FILE='F2D_FDM_50.DAT' )
C
OPEN ( 7, FILE='F2D_FDM_R400_UM20_FIG.DAT' ) ! DT = 0.02D0
OPEN ( 8, FILE='F2D_FDM_R400_UM20_OUT.DAT' )
OPEN ( 10, FILE='F2D_FDM_20.DAT' )
C
*** OPEN ( 7, FILE='F2D_FDM_R400_UM50_FIG.DAT' )
*** OPEN ( 8, FILE='F2D_FDM_R400_UM50_OUT.DAT' )
*** OPEN ( 10, FILE='F2D_FDM_50.DAT' )
C
C ------------------------------------------------
C
CALL PRSB ( PPT, WPT, EPT, PS1, WP1, PS, FF, PSM, U,
1 V, UB, VB, WD, WZ, HH, BX, BY, XZ,
2 YZ, RAF, RAS, HX, HY, X, Y, IX1, IY1 )
C
CALL WUCF ( HX, HY, UF1, UF2, UF3, UNA, UNB, UNC, UND,
1 UNE, UAP, UBP, UCP, UDP, UEP, UAM, UBM, UCM,
2 UDM, UEM, X, Y, IX1, IY1 )
C
CALL WVCF ( HX, HY, VF1, VF2, VF3, VNA, VNB, VNC, VND,
1 VNE, VAP, VBP, VCP, VDP, VEP, VAM, VBM, VCM,
2 VDM, VEM, X, Y, IX1, IY1 )
C
LOP = 0
C
C --------------------------------------------------------------------
999 LOP = LOP + 1
C
CALL PSCL ( PS1, PS, PSM, WP1, PPT, EPT, 1.5D0, HX,
1 HY, RDP, X, Y, IX1, IY1 )
C
CALL WFCA ( WP1, PS1, WPT, HX, HY, RAS, RAF, X, Y,
1 IX1, IY1 )
C
CALL VARC ( WP1, FF, WD, WZ, VRW, II, JJ, X, Y )
C
CALL VARC ( PS1, PSM, WD, WZ, VRP, II, JJ, X, Y )
C
CALL WUCL ( WP1, FF, PS1, WPT, EPT, U, HX, HY,
1 UF1, UF2, UF3, UNA, UNB, UNC, UND, UNE,
2 UAP, UBP, UCP, UDP, UEP, UAM, UBM, UCM,
3 UDM, UEM, X, Y, IX1, IY1 )
C
CALL WVCL ( WP1, FF, PS1, WPT, EPT, V, HX, HY,
1 VF1, VF2, VF3, VNA, VNB, VNC, VND, VNE,
2 VAP, VBP, VCP, VDP, VEP, VAM, VBM, VCM,
3 VDM, VEM, X, Y, IX1, IY1 )
C
CALL TMCL ( WP1, FF, WPT, EPT, U, V, HH, HX,
1 HY, BX, BY, X, Y, IX1, IY1 )
C
CALL CHCK ( U, V, UB, VB, VRW, VRT, VRP, X, Y,
1 *900 )
C
GO TO 999
900 CONTINUE
C
CALL POUT ( XX, YY, NDE, PE, PS, UB, VB, PS1, U,
1 V, XZ, YZ, WP1, HX, HY, JU, IV, II,
2 JJ, X, Y, IXY, IXY1, IX1, IY1 )
C
CLOSE ( 7)
CLOSE ( 8)
CLOSE (10)
C
STOP
END
C **********************************************************************
C
SUBROUTINE CHCK ( U, V, UB, VB, VRW, VRT, VRP, IX, IY, * )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
COMMON /UDAT/ AS, ITN
COMMON /CTRL/ MSH, IHS, WEI
COMMON /DEF / VAW(50000), VARP(50000)
COMMON /IDAT/ DT, LOPE, REN, LOP, NPL

C
DIMENSION U( IX,IY), V( IX,IY), UB( IX,IY), VB( IX,IY)
C
IF ( LOP.EQ.1) NPL = 0
IF ( LOP.EQ.1) GO TO 123
C
IF ( VRW.LE.0.01D0 .AND. VRP.LE.0.01D0 ) THEN
C
WRITE(6,*) ' '
WRITE(6,*) ' *** Stationary condition reached at Loop =', LOP
WRITE(8,*) ' '
WRITE(8,*) ' *** Stationary condition reached at Loop =', LOP
C
NPL = NPL + 1
VAW( NPL) = VRW
VARP( NPL) = VRP
C
WRITE(6,2000) LOP, ITN, VRW, VRP
WRITE(8,2000) LOP, ITN, VRW, VRP
C
RETURN 1
C
END IF
C ------------------------------
C
NPR = 100
C ===============
C
IF ( MOD( LOP,NPR).EQ.0 ) THEN
C ===================================
C
NPL = NPL + 1
VAW( NPL) = VRW
VARP( NPL) = VRP
C
WRITE(6,2000) LOP, ITN, VRW, VRP
WRITE(8,2000) LOP, ITN, VRW, VRP
2000 FORMAT(' LOP=', I6,' It=', I3,' VrW=',F7.2,'% ',' VrP=',F7.2,
1 '% ')
END IF
C
IF ( LOP.GE.LOPE) RETURN 1
C
123 CONTINUE
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
RETURN
END
C **********************************************************************
C
SUBROUTINE TTLE
C
WRITE(6,2000)
WRITE(6,2100)'***************************************************'
WRITE(6,2100)' '
WRITE(6,2100)' *** F 2 D - FDM ( SV ) *** '
WRITE(6,2100)' '
WRITE(6,2100)' 2D Fluid Analysis by the 4th order FDM '
WRITE(6,2100)' '
WRITE(6,2100)' ( V.3.5 ) '
WRITE(6,2100)' '
WRITE(6,2100)' Copyright 2011 : Yasuhiro MATSUDA '
WRITE(6,2100)' '
WRITE(6,2100)'***************************************************'
C
2000 FORMAT(/)
2100 FORMAT(A55)
C
RETURN
END
C **********************************************************************
C
SUBROUTINE DINP ( HX, HY, RAF, RAS, IX, IY, IX1, IY1 )
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION HX( IX1), HY( IY1), RAF( IX1), RAS( IY1)
C
COMMON /IDAT/ DT, LOPE, REN, LOP, NPL
C
C --------------------------------
DO 100 I = 1, IX-1
C
RAF( I) = HX( I)* HX( I)
C
100 CONTINUE
C --------------------------------
C
DO 200 J = 1, IY-1
C
RAS( J) = HY( J)* HY( J)
C
200 CONTINUE
C --------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE GINP ( WP1, PS1, WPT, PPT, EPT, IX, IY )
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION INDX( 5), INDY( 5), INDE( 5)
DIMENSION WP1( IX,IY), PS1( IX,IY)
C
LOGICAL*1 WPT, PPT, EPT
DIMENSION WPT( IX,IY,4), PPT( IX,IY,4), EPT( IX,IY,4)
C
NAMELIST /GDAT/ OMG, PSI, SLP, MOV, CNT, PRS
LOGICAL*1 OMG, PSI, SLP, MOV, CNT, PRS
C
C ----- SPEC Data READ ---
C
IER = 0
C ----------------------------
DO 100 I = 1, IX
DO 100 J = 1, IY
C
WPT( I,J,1) = .TRUE.
PPT( I,J,1) = .TRUE.
C
100 CONTINUE
C ----------------------------
C
C ----- All points are initially in the field ---
DO 120 I = 1, IX
C
EPT( I, 1,2) = .TRUE.
EPT( I,IY,1) = .TRUE.
C
120 CONTINUE
C -----------------------------------------------
DO 130 I = 1, IY
C
EPT( 1, I,4) = .TRUE.
EPT( IX,I,3) = .TRUE.
C
130 CONTINUE
C -----------------------------------------------
C
C ----- Edge points ---
C
200 CONTINUE
OMG = .FALSE.
PSI = .FALSE.
SLP = .FALSE.
MOV = .FALSE.
CNT = .FALSE.
PRS = .FALSE.
C
IDA = 1
C
READ(10,GDAT,END=999)
C ==========================
C
IF ( OMG) GO TO 205
IF ( PSI) GO TO 205
IF ( PRS) GO TO 205
IF ( CNT) GO TO 205
IF ( MOV) GO TO 205
C
205 CONTINUE
C ----------------------------------------------------------
DO 300 IJK = 1, IDA
C
ISLP = -10
C
C ------------------------------------------------
DO 320 L = 1, 5
C
INDX( L) = 0
INDY( L) = 0
INDE( L) = 0
320 CONTINUE
C ------------------------------------------------
C
WVAL = 0.D0
PVAL = 0.D0
TVAL = 0.D0
PRVL= 0.D0
C
READ(10,1000,END=999,ERR=225 ) ISLP, WVAL, PVAL, TVAL,
1 ( INDX( L), INDY( L), INDE( L), L =1,5 ), PRVL
1000 FORMAT( I2,1X,F6.3,1X,2(F6.3,1X),5(3I3,1X), F6.3)
C
IPT = 0
C
C IF ISLP = -1 Boundary : Downward
C IF ISLP = 0 Boundary : Horizontal
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 IER = 1
WRITE(6,2100) ISLP, WVAL, PVAL, TVAL, ( INDX( L), INDY( L),
1 INDE( L), L= 1, 5 ), PRVL
2100 FORMAT(5X, I2, 1X, 3( F6.3,1X) , 5( 3I3,1X) ,F6.3)
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 : Downward. IPT = 2, Horizontal
C IPT = 3, LINE : Upward. IPT = 4, Vertical
C
C -----------------------------------------------
DO 340 L = 1, 5
C
IUPR = INDE( L)
GO TO ( 240, 240, 240, 235 ) IPT
C
235 CONTINUE
ILOW = INDY( L)
GO TO 245
C
240 ILOW = INDX( L)
245 I = INDX( L)
J = INDY( L)
C
C --------------------------------------
DO 290 INDEX = ILOW, IUPR
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 ---- IF I or J = 0 : No Boundary specified ---
C
255 IER = 1
WRITE(6,2200) L, IX, IY
2200 FORMAT(' * Error in SPEC No.', I1, 2X,' IX =', I3,' IY =', I3)
WRITE ( 16,2100) ISLP, WVAL, PVAL, TVAL,
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) = MOV
WP1( I,J) = WVAL
C
265 IF ( .NOT.PSI) GO TO 270
PPT( I,J,1) = .FALSE.
PPT( I,J,2) = SLP
PPT( I,J,3) = MOV
PS1( I,J) = PVAL
C
270 CONTINUE
IF ( .NOT.CNT) GO TO 278
PPT( I,J,4) = .TRUE.
278 CONTINUE
IF (.NOT.PRS) GO TO 279
279 CONTINUE
IF (.NOT.MOV) GO TO 280
WPT( I,J,3) = .TRUE.
280 GO TO ( 286, 288, 287, 285 ), IPT
285 J = J + 1
GO TO 290
C
286 J = J - 1
GO TO 288
C
287 J = J + 1
288 I = I + 1
C
290 CONTINUE
C --------------------------------------
C
340 CONTINUE
C ------------------------------------------------
C
300 CONTINUE
C ----------------------------------------------------------
GO TO 200
C
999 RETURN
END
C **********************************************************************
C
SUBROUTINE IGSS ( PS1, PPT, XZ, YZ, IX, IY )
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION PS1( IX,IY), XZ ( IX), YZ ( IY)
C
LOGICAL*1 PPT( IX,IY,4)
C
COMMON /CTRL/ MSH, IHS, WEI
COMMON /IDAT/ DT, LOPE, REN, LOP, NPL
C
N = IY - 1
C
DO 100 I = 1, IX
DO 100 J = 1, IY
C
IF (.NOT.PPT( I,J,1) ) GO TO 100
C
PS1( I,J) = PS1( IX,J)
C
100 CONTINUE
C
RETURN
END
C **********************************************************************
C
SUBROUTINE MAXS ( U, IX, IY, UMX, IU, JU, IH )
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION U( IX,IY)
C
IF ( IH.NE.1) GO TO 112
C
UMX = -1000.D0
C
C -------------------------------------
DO 100 I = 1, IX
DO 100 J = 1, IY
C
IF ( U(I,J).LT.UMX) GO TO 100
IU = I
JU = J
UMX = U( I,J)
C
100 CONTINUE
C -------------------------------------
GO TO 900
C
112 CONTINUE
C
IF ( IH.NE.2) GO TO 113
UMX = 0.D0
I = IU
C -------------------------------------
DO 200 J = 1, IY
C
UV = U( I,J) - UMX
IF ( UV.LE.0.D0 ) GO TO 200
IU = I
JU = J
UMX = U( I,J)
C
200 CONTINUE
C ------------------------------------
GO TO 900
C
113 CONTINUE
UMX = 0.D0
J = JU
C ------------------------------------
DO 300 I = 1, IX
C
UV = U( I,J) - UMX
IF ( UV.LE.0.D0 ) GO TO 300
IU = I
JU = J
UMX = U( I,J)
C
300 CONTINUE
C ------------------------------------
C
900 RETURN
END
C **********************************************************************
C
SUBROUTINE POUT( XX, YY, NDE, PE, PS, UB, VB, PS1, U,
1 V, XZ, YZ, WP1, HX, HY, JU, IV, II,
2 JJ, IX, IY, IXY, IXY1, IX1, IY1 )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
COMMON /CTRL/ MSH, IHS, WEI
COMMON /IDAT/ DT, LOPE, REN, LOP, NPL
COMMON /DEF / VAW(50000), VARP(50000)
C
DIMENSION XX( IXY), YY( IXY), PE( IXY), PS( IX,IY),
1 UB( IX,IY), VB( IX,IY), PS1( IX,IY), U ( IX,IY),
2 V ( IX,IY), NDE( IXY1,4)
C
DIMENSION XZ( IX), YZ( IY), WP1( IX,IY), HX( IX1), HY( IY1)
C
NNP = IX* IY
NEL = ( IX-1)* ( IY-1)
C
CALL STND ( XX, YY, NDE, NNP, NEL, XZ, YZ, IX, IY )
C
CALL CKUV ( U, UB, DT, DUX, IPU, JPU, IX, IY )
C
CALL CKUV ( V, VB, DT, DVX, IPV, JPV, IX, IY )
C
CALL STMP ( IX, IY, PS1, PE, NNP )
C
C ----- Plotting Data ---
C
WRITE(7,2000) REN, DT, LOP, MSH, NPL
WRITE(8,2000) REN, DT, LOP, MSH, NPL
2000 FORMAT( F8.2, F7.5, I5, I2, I3)
WRITE(7,2100) ( PE( I), I = 1, IXY )
WRITE(7,2100) ( VAW( I), I = 1, NPL )
WRITE(7,2100) ( VARP( I), I = 1, NPL )
WRITE(7,2100) ( XX( I), I = 1, IXY )
WRITE(7,2100) ( YY( I), I = 1, IXY )
WRITE(7,2100) ( ( U( I,J), J = 1, IY ), I = 1, IX)
WRITE(7,2100) ( ( V( I,J), J = 1, IY ), I = 1, IX)
WRITE(7,2100) ( HX( I), I = 1, IX-1 )
WRITE(7,2100) ( HY( I), I = 1, IY-1 )
2100 FORMAT(10F8.3)
WRITE(7,2400) ( ( NDE( I,J), J = 1, 4), I = 1, IXY1 )
2400 FORMAT(20I4)
C
C ------------------------------------------------
C
CALL DTOP ( PS1, U, V, XZ, YZ, IX, IY )
C
C ------------------------------------------------
C
WRITE(6,*) ' '
WRITE(6,*) ' *** End ***'
WRITE(8,*) ' '
WRITE(8,*) ' *** End ***'
C
RETURN
END
C **********************************************************************
C
SUBROUTINE DTOP ( PS1, U, V, XZ, YZ, IX, IY )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
COMMON /IDAT/ DT, LOPE, REN, LOP, NPL
COMMON /CTRL/ MSH, IHS, WEI
C
DIMENSION PS1( IX,IY), U( IX,IY), V( IX,IY), XZ( IX), YZ( IY)
C
C ----- Detailed output ---
C
CALL MAXS ( U, IX, IY, RLU, IAU, JAU, 1 )
C
CALL MAXS ( V, IX, IY, RLV, IAV, JAV, 1 )
C
C ------ Attention : Approx. for equal length of meshes ---
C
DY = 1.D0/ DFLOAT( IY-1)
C ==============================
C
RW = DT/( DY**2* REN )
BV = RLU* DT / DY
C ==============================
C
IF ( RLV.GT.RLU) BV = RLV* DT/ DY
C
WRITE(6,*) ' '
WRITE(6,2000) RW, BV
WRITE(8,*) ' '
WRITE(8,2000) RW, BV
2000 FORMAT(' * r2, b =',2F9.4)
C
C ------------------------------------
DO 100 I = 1, IX
DO 100 J = 1, IY
C
PS1( I,J) = DABS( PS1( I,J) )
C
100 CONTINUE
C ------------------------------------
C
IMU = ( IX+1)/ 2
JMV = ( IY+1)/ 2
PMID = PS1( IMU, JMV )
C
C -----------------------------------------------------
C
CALL MAXS ( PS1, IX, IY, PMX, IP, JP, 1 )
C
C -----------------------------------------------------
C
XP = XZ( IP)
YP = YZ( JP)
C
WRITE(6,*) ' '
WRITE(6,2100) PMX, XP, YP, PMID
WRITE(8,*) ' '
WRITE(8,2100) PMX, XP, YP, PMID
2100 FORMAT(' * P_max.=',F8.4,' at Xp,Yp =',2F8.4, ' * P_mid.=',F8.4)
C
XAU = XZ( IAU)
YAU = YZ( JAU)
C
XAV = XZ( IAV)
YAV = YZ( JAV)
C
WRITE(6,2200) RLU, XAU, YAU, RLV, XAV, YAV
WRITE(8,2200) RLU, XAU, YAU, RLV, XAV, YAV
2200 FORMAT(/,' * Umax.=',F8.3,' at X, Y =',2F8.4,/,' * Vmax.=',
1 F8.3,' at X, Y =',2F8.4)
C
C -------------------------------------------------
C
CALL MAXS( U, IX, IY, UMX, IMU, JU, 2)
C
CALL MAXS( V, IX, IY, VMX, IV, JMV, 3)
C
C -------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE PSCL ( PS1, PS, PSM, WP1, PPT, EPT, Q, HX, HY,
1 RDP, IX, IY, IX1, IY1 )
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION PS1( IX,IY), PS( IX,IY), PSM( IX,IY), WP1( IX,IY),
1 RDP( IX,IY), HX( IX1), HY ( IY1)
C
LOGICAL*1 PPT( IX,IY,4), EPT( IX,IY,4), DRT
C
COMMON /UDAT/ AS, ITN
COMMON /IDAT/ DT, LOPE, REN, LOP, NPL
COMMON /CM1 / PA1( 51,51), PA2( 51,51), PA3( 51,51),
1 PB1( 51,51), PB2( 51,51), PB3( 51,51)
C
FS = 1.D0
CFI = 1.D0
C
AFS = AS* FS
C =================
C
RUU = 0.D0
RUD = 0.D0
C
ITN = 0
EPS = 0.00005D0
C
C --------------------------------------------
DO 100 I = 1, IX
DO 100 J = 1, IY
C
PS1( I,J) = 2.D0* PS( I,J) - PSM( I,J)
C
100 CONTINUE
C --------------------------------------------
C
DRT = .FALSE.
JDY = IY + 1
C
105 CONTINUE
IDC = 0
ITN = ITN + 1
C
IF ( ITN.GT.500) GO TO 98
C
C ----------------------------------------------------------
DO 200 JJ = 1, IY
C
10 CONTINUE
C
DRT = .NOT. DRT
C
IF ( DRT) J = JJ
IF ( .NOT.DRT) J = JDY - JJ
C
C ------------------------------------------------
DO 300 I = 1, IX
C
DX = HX( I-1)* HX( I)* ( HX( I-1) + HX( I) )
DY = HY( J-1)* HY( J)* ( HY( J-1) + HY( J) )
C
DIXY = ( HY( J-1) + HY( J) )*DX + ( HX(I-1) + HX(I) )*DY
C
IF ( PPT( I,J,1) ) GO TO 1
IF ( PPT( I,J,2) ) GO TO 300
GO TO 300
C
1 CONTINUE
C
IF ( EPT( I,J,3) ) GO TO 2
IF ( EPT( I,J,4) ) GO TO 12
IF ( EPT( I,J,1) ) GO TO 22
IF ( EPT( I,J,2) ) GO TO 32
C
RDP( I,J) = ( PA1( I,J)*PS1( I-1,J) + PA3( I,J)*PS1( I+1,J)
1 + PB1( I,J)*PS1( I,J-1) + PB3( I,J)*PS1( I,J+1)
2 + WP1( I,J) )/ ( PA2( I,J) + PB2( I,J) )
C
RAB2 = DABS( RDP( I,J) - PS1( I,J) )
C
PS1( I,J) = PS1( I,J) + Q *( RDP( I,J) - PS1( I,J) )
C
IF ( RAB2.GT.EPS ) GO TO 5
IF ( ITN.LT.10 ) GO TO 5
GO TO 300
C
5 IDC = 1
IXN = I
IYN = J
C
GO TO 300
C
2 IF ( PPT( I,J,4)) GO TO 8
C
C ----- Periodic right ---
C
PS1( I,J) = PS1( I-(IX-1),J)
GO TO 300
C
C ----- Continuative right ---
C
8 CONTINUE
C
PS1( I,J)
1 = ( HY( J)*PS1(I,J-1) + HY( J-1)*PS1(I,J+1)
2 + 0.5D0*DY* WP1(I,J) )/ ( HY( J-1) + HY( J) )
GO TO 300
C
12 IF ( PPT( I,J,4) ) GO TO 18
C
C ----- Periopdic left ---
C
DELP = ( PS1( I+1,J)+PS1( I-1+( IX-1), J)
1 + FS*( PS1( I,J+1)+PS1( I,J-1)+AS* WP1( I,J) ) )*CFI
C
RDL = ( ( HX( I-1)*PS1( I+1,J) + HX( I)*PS1( I-1+(IX-1), J) ) * DY
1 + ( HY( J-1)*PS1( I,J+1) + HY( J)*PS1( I,J-1) ) * DX
C
2 + 0.5D0 * DX * DY * WP1( I,J) )/ DIXY
C
PS1( I,J) = PS1( I,J) + Q*( RDL - PS1( I,J) )
GO TO 300
C
C ----- Continuative left ---
C
18 PS1( I,J)
1 = ( ( HX( I-1)*PS1( I,J) + HX( I)*PS1( I,J) ) *DY
2 + ( HY( J-1)*PS1( I,J+1) + HY( J)*PS1( I,J-1) )*DX
3 + 0.5D0*DX*DY* WP1( I,J) )/ DIXY
GO TO 300
C
22 IF ( PPT( I,J,4)) GO TO 28
C
PS1( I,J) = PS1( I,J-(IY-1) )
GO TO 300
C
28 PS1( I,J) = ( HX( I)* PS1( I-1,J) + HX( I-1)* PS1( I+1,J)
1 + 0.5D0* DX* WP1( I,J) )/ ( HX( I-1) + HX( I) )
GO TO 300
C
32 IF ( PPT( I,J,4)) GO TO 38 ! Continuative
C
DELP = ( PS1( I,J+1) + PS1( I,J-1+( IY-1) )
1 + ( PS1( I+1,J) + PS1( I-1,J) + AFS* WP1( I,J) )
2 / FS )* CFI / FS
C
PS1( I,J) = PS1( I,J) + Q*( DELP - PS1( I,J) )
GO TO 300
C
38 PS1( I,J) = ( HX( I)* PS1( I-1,J) + HX( I-1)* PS1( I+1,J)
1 + 0.5D0* DX* WP1( I,J) )/ ( HX( I-1) + HX( I) )
C
GO TO 300
C
300 CONTINUE
C -------------------------------------------------
IF ( DRT) GO TO 10
C
200 CONTINUE
C ------------------------------------------------------------
C
IF ( DABS( RUD).LE.1.D-6) GO TO 101
C
PSS = 0.D0
C -------------------------------------------------
DO 205 I = 1, IX
C
PSS = PSS + PS1( I,2) + 0.5D0* AS* WP1( I,1)
C
205 CONTINUE
C -------------------------------------------------
C
UL = 0.D0
C
PSS = ( PSS/ IX)
C
C -------------------------------------------------
DO 207 I = 1, IX
C
PS1( I,1) = PSS
C
207 CONTINUE
C -------------------------------------------------
C
101 IF ( DABS( RUU ).LE.1.D-6) GO TO 102
C
PSS = 0.D0
C ----------------------------------------------------
DO 209 I = 1, IX
C
PSS = PSS + PS1( I,IY-1) + 0.5D0* AS* WP1( I,IY)
C
209 CONTINUE
C -----------------------------------------------------
C
UU = 0.D0
DO 211 I = 1, IX
211 PS1( I,IY) = PSS
C
102 CONTINUE
IF ( IDC.GT.0) GO TO 105
GO TO 777
C
98 CONTINUE
WRITE (6,2000) LOP
2000 FORMAT(/' *** Diverged ! ITN >2000 in PSCL*Incompl. sol.of PS
1I:* STOP at LOP =',I4)
C
STOP
C
777 CONTINUE
C
C --------------------------
DO 120 I = 1, IX
DO 120 J = 1, IY
C
PS( I,J) = PS1( I,J)
C
120 CONTINUE
C --------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE PRSB ( PPT, WPT, EPT, PS1, WP1, PS, FF, W,
1 U, V, UB, VB, WD, WZ, HH, BX,
2 BY, XZ, YZ, RAF, RAS, HX, HY, IX,
3 IY, IX1, IY1 )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
LOGICAL*1 PPT( IX,IY,4), WPT( IX,IY,4), EPT( IX,IY,4)
C
DIMENSION PS1( IX,IY), WP1( IX,IY), PS( IX,IY), FF( IX,IY),
1 W ( IX,IY), U ( IX,IY), V ( IX,IY), UB( IX,IY),
2 VB ( IX,IY), WD ( IX,IY), WZ( IX,IY), HH( IX,IY),
3 BX ( IX,IY), BY ( IX,IY), XZ( IX), YZ( IY)
C
DIMENSION RAF( IX1), RAS( IY1), HX( IX1), HY( IY1)
C
C
CALL TTLE
C
C ----- Data Input -------
C
CALL CDATA ( IY )
C
C ----- Initialization ---
C
CALL SBINT( PPT, WPT, EPT, PS1, WP1, PS, FF, W, U, V,
1 UB, VB, WD, WZ, HH, BX, BY, XZ, YZ, IX,
2 IY )
C
C ----- Mesh generation ---
C
CALL GMSH ( XZ, YZ, HX, HY, IX, IY, IX1, IY1 )
C
CALL DINP ( HX, HY, RAF, RAS, IX, IY, IX1, IY1 )
C
C ----- Geometry Data Read ---
C
CALL GINP ( WP1, PS1, WPT, PPT, EPT, IX, IY )
C
CALL IGSS ( PS1, PPT, XZ, YZ, IX, IY )
C
CALL PSCF( IX, IY, HX, HY, IX1, IY1 )
C
CALL TMCOF( HX, HY, IX, IY, IX1, IY1 )
C
RETURN
END
C **********************************************************************
C
SUBROUTINE CDATA ( IY )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
COMMON /CTRL/ MSH, IHS, WEI
COMMON /IDAT/ DT, LOPE, REN, LOP, NPL
COMMON /UDAT/ AS, ITN
C
C -------------------------------------------
C
C MSH = 1 : Mesh of Equal length
C = 2 : Mesh of Unequal length
C
C IHS = 1 : With Correction
C WEI : Weighting parameter
C
C -------------------------------------------
C
WRITE(6,*) ' '
WRITE(6,*) ' ** MSH = 1 : Mesh of Equal length 2 : Mesh of Une
1qual length'
C
** READ (5,*) MSH
C --------------------
MSH = 2
C
IHS = 1
C
WEI = 0.5D0 ! Better than = 0 for the equal length meshes ---
C
WRITE(6,2233) MSH, IHS, WEI
WRITE(8,2233) MSH, IHS, WEI
2233 FORMAT(/,' * MSH =', I2,' IHS =', I2,' WEI =',F5.2)
C
WRITE(6,*) ' '
WRITE(6,*) ' ** Re = ?'
** READ (5, *) REN
REN = 400.D0
C =================
C
WRITE(6,*) ' '
WRITE(6,1333) REN
WRITE(8,1333) REN
1333 FORMAT(' *** Re =', F8.2)
C
VSCF = 1.D0/ REN
C ------ Attention -----
C
WRITE(6,*) ' '
WRITE(6,*) ' ** DT = ? & LOP = ?'
*** READ(5,*) DT, LOPE
DT = 0.02D0
LOPE = 1200
C
WRITE(6,1125) DT, LOPE
WRITE(8,1125) DT, LOPE
1125 FORMAT(/,' *** Cavity Flow Problem * DT =',F8.5,' LOP =',I5 )
C
C ----- Attention : Approx. For equal length of meshes --- A: DY ---
C
A = 1.D0/ DFLOAT( IY-1 )
AS = A* A
C
RETURN
END
C **********************************************************************
C
SUBROUTINE GMSH ( XZ, YZ, HX, HY, IX, IY, IX1, IY1 )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION XZ( IX), YZ ( IY), HX( IX1), HY( IY1)
C
COMMON /CTRL/ MSH, IHS, WEI
C
C ----- To set positions of NDEs ---
C
C MSH = 1 : Meshes of equal length
C
C MSH = 2 : Meshes of unequal length
C ( AA : Init. coeff., RR : Ratio )
C
WRITE(6,*) ' '
C
IF ( MSH.EQ.1) GO TO 200
IF ( MSH.EQ.2) GO TO 230
C *******************************
C
200 CONTINUE
C
WRITE(6,*) '* Meshes of equal length '
WRITE(8,*) '* Meshes of equal length '
C
DO 25 I = 1, IX
C
XZ( I) = DFLOAT( I-1)/ DFLOAT( IX-1)
C ------------------------------------------
25 CONTINUE
C
DO 26 J = 1, IY
C
YZ( J) = DFLOAT( J-1)/ DFLOAT( IY-1)
C ------------------------------------------
26 CONTINUE
C
GO TO 34
C
* 220 CONTINUE
C
* WRITE(6,*) '* Meshes of unequal length ( 1) ( Equal ratio ) '
* WRITE(8,*) '* Meshes of unequal length ( 1) ( Equal ratio ) '
* WRITE(6,*) ' '
C
* AA = 0.01926D0
* RR = 1.2D0
C
* MDX = ( IX + 1)/ 2
* MX1 = MDX - 1
* MDY = ( IY + 1)/ 2
* MY1 = MDY - 1
* WRITE(6,*) ' MDX, MX1, MDY, MY1=', MDX, MX1, MDY, MY1
C
* XZ( 1) = 0.D0
* XZ( MDX) = 0.5D0
* XZ( IX) = 1.D0
* YZ( 1) = 0.D0
* YZ( MDY) = 0.5D0
* YZ( IY) = 1.D0
C
* DO 35 I = 2, MX1
C
* XZ( I) = AA* ( RR** ( I-1) - 1.D0 )/ ( RR - 1.D0 )
* XZ( IX+1-I) = 1.D0 - XZ( I)
C
* 35 CONTINUE
C
* WRITE(6,*) ' * XZ=', ( XZ( I), I=1, IX )
* STOP
C
* DO 36 I = 2, MY1
C
* YZ( I) = AA* ( RR**( I-1) - 1.D0)/ ( RR - 1.D0 )
* YZ( IY+1-I) = 1.D0 - YZ ( I)
C
* 36 CONTINUE
C
* GO TO 34
C
230 CONTINUE
C
WRITE(6,*) '*** Meshes of unequal length ' ! Ref: JSME - Paper
WRITE(8,*) '*** Meshes of unequal length '
C
AA = 2.D0
DX = 2.D0 / DFLOAT( IX - 1)
DY = 2.D0 / DFLOAT( IY - 1)
C
MDX = ( IX+1)/2
MX1 = MDX-1
MDY = ( IY+1)/2
MY1 = MDY-1
C
XZ( MDX) = 0.5D0
YZ( MDY) = 0.5D0
C
DO 37 I = 1, MX1
C
XZ( I) = ( DEXP( AA* DFLOAT( I-1)* DX ) - 1.D0 )
1 / ( 2.D0* ( DEXP( AA ) - 1.D0 ) )
C
XZ( IX+1-I) = 1.D0 - XZ( I)
C
37 CONTINUE
C
*** WRITE(6,*) ' * XZ=', ( XZ( I), I=1, IX )
C
DO 38 J = 1, MY1
C
YZ( J) = ( DEXP( AA* DFLOAT( J-1)* DY) - 1.D0 )
1 / ( 2.D0* ( DEXP( AA) - 1.D0 ) )
C
YZ( IY+1-J) = 1.D0 - YZ( J)
C
38 CONTINUE
C
C -------------------------------------------
C
34 CONTINUE
C
DO 27 I = 1, IX-1
C
HX( I) = DABS( XZ( I+1) - XZ( I) )
27 CONTINUE
C
DO 28 J = 1, IY-1
C
HY( J) = DABS( YZ( J+1) - YZ( J) )
28 CONTINUE
C -------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE SBINT( PPT, WPT, EPT, PS1, WP1, PS, FF, W, U, V, UB,
1 VB, WD, WZ, HH, BX, BY, XZ, YZ, IX, IY )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
LOGICAL*1 PPT( IX,IY,4), WPT( IX,IY,4), EPT( IX,IY,4)
C
DIMENSION PS1( IX,IY), WP1( IX,IY), PS( IX,IY), FF( IX,IY),
1 W ( IX,IY), U ( IX,IY), V ( IX,IY), UB( IX,IY),
2 VB ( IX,IY), WD ( IX,IY), WZ( IX,IY), HH( IX,IY),
3 BX ( IX,IY), BY ( IX,IY), XZ( IX), YZ( IY)
C
C *********************************
DO 100 I = 1, IX
C
DO 100 J = 1, IY
C
DO 50 K = 1, 4
C
PPT( I,J,K) = .FALSE.
WPT( I,J,K) = .FALSE.
EPT( I,J,K) = .FALSE.
C
50 CONTINUE
C
PS1( I,J) = 0.D0
WP1( 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
C
WD ( I,J) = 0.D0
WZ ( I,J) = 0.D0
HH ( I,J) = 0.D0
BX ( I,J) = 0.D0
BY ( I,J) = 0.D0
C
XZ ( I) = 0.D0
YZ ( J) = 0.D0
C
100 CONTINUE
C ****************************
C
RETURN
END
C **********************************************************************
C
SUBROUTINE VARC ( WP1, WB, WD, WZ, VRW, II, JJ, IX, IY )
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION WP1( IX,IY), WB( IX,IY), WD( IX,IY), WZ( IX,IY)
C
DO 100 I = 1, IX
DO 100 J = 1, IY
C
WD( I,J) = 0.D0
WZ( I,J) = 0.D0
C
100 CONTINUE
C
WL1 = 0.D0
WL2 = 0.D0
WDEL = 0.D0
VRW = 0.D0
C
DO 200 I = 1, IX
DO 200 J = 1, IY
C
WZ ( I,J) = DABS( WP1( I,J) )
IF ( WZ( I,J).GT.WL1) WL1 = WZ( I,J)
C
200 CONTINUE
C
WL2 = WL1*0.01D0
C
DO 300 I = 1, IX
DO 300 J = 1, IY
C
ABSW = DABS ( WP1( I,J) )
C
IF ( ABSW.LE.WL2) GO TO 300
WD( I,J) = DABS ( ( WP1( I,J)-WB( I,J) )/ WP1( I,J) )
C
IF ( WD( I,J).GT.WDEL ) THEN
WDEL = WD( I,J)
II = I
JJ = J
END IF
C
300 CONTINUE
C
VRW = WDEL*100.D0
C
DO 400 I = 1, IX
DO 400 J = 1, IY
C
WB( I,J) = WP1( I,J)
C
400 CONTINUE
C
RETURN
END
C **********************************************************************
C
SUBROUTINE CKUV ( U, UB, DT, DUX, II, JJ, IX, IY )
C
C ----- Unsteady term ---
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION U( IX,IY), UB( IX,IY)
C
DUX = 0.D0
C
DO 100 I = 1, IX
DO 100 J = 1, IY
C
DIFU = DABS ( ( U( I,J)-UB( I,J) )/ DT )
C
IF ( DIFU.LT.DUX) GO TO 100
C
DUX = DIFU
II = I
JJ = J
C
100 CONTINUE
C
RETURN
END
C **********************************************************************
C
SUBROUTINE TMCL ( WP1, FF, WPT, EPT, U, V, HH, HX, HY,
1 BX, BY, IX, IY, IX1, IY1 )
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION WP1( IX,IY), FF( IX,IY), U ( IX,IY),
1 V ( IX,IY), HH( IX,IY), BX( IX,IY),
2 BY ( IX,IY), HX( IX1), HY( IY1)
C
COMMON /CM2/ TM1( 51), TM2( 51), TM3( 51), TY1( 51), TY2( 51),
1 TY3( 51)
C
LOGICAL*1 WPT( IX,IY,4), EPT( IX,IY,4)
C
COMMON /IDAT/ DT, LOPE, REN, LOP, NPL
COMMON /UDAT/ AS, ITN
COMMON /CTRL/ MSH, IHS, WEI
C
C ***************
FS = 1.D0
CF = 1.D0
C ***************
C
C ----- Attention : Approx. For equal length of meshes ---
C
DO 185 I = 1, IX
DO 185 J = 1, IY
C
BX( I,J) = U( I,J)* DT* DFLOAT( IX1)
BY( I,J) = V( I,J)* DT* DFLOAT( IY1)
C
185 CONTINUE
C
RH = 1.D0 / REN ! For Vorticity Equation ==> Prandtle Number
C ******************
C
IF ( IHS.NE.1) GO TO 81
C
C ----- Attention : Approx. For equal length meshes ---
C
RX = RH* DT* ( IX1**2)
RY = RH* DT* ( IY1**2)
C
C ---------------------------------------------
C
DO 280 I = 1, IX
DO 280 J = 1, IY
C
HH( I,J) = ( 1.D0 + ( BX( I,J)* BY( I,J) )/ ( RX+RY) )* RH* DT
C
280 CONTINUE
C
GO TO 85
C ---------------------------------
C
81 CONTINUE
C
DO 82 I = 1, IX
DO 82 J = 1, IY
C
HH( I,J) = RH*DT
C
82 CONTINUE
C
85 CONTINUE
C
DO 90 I = 1, IX
DO 90 J = 1, IY
C
IF ( WPT( I,J,1) ) GO TO 90
WP1( I,J) = FF( I,J)
C
90 CONTINUE
C
C ***************************************************
DO 700 J = 1, IY
C
C *****************************************
DO 500 I = 1, IX
C
IF ( WPT( I,J,1) ) GO TO 1
IF ( WPT( I,J,2) ) GO TO 101
IF ( WPT( I,J,3) ) GO TO 200
C
C ----- Not Fluid, Not conducting, Not insulated ---
C
GO TO 101
C
1 CONTINUE
C
IF ( EPT( I,J,3) ) GO TO 2
IF ( EPT( I,J,4) ) GO TO 12
IF ( EPT( I,J,1) ) GO TO 22
IF ( EPT( I,J,2) ) GO TO 32
C
C ----- Normal interior to fluid point ---
C
WP1( I,J) = WP1( I,J) +2.D0*HH( I,J)
1 *( TM1( I)*FF( I-1,J) - TM2( I)*FF( I,J)
C
2 + TM3( I)*FF( I+1,J) + TY1( J)*FF( I,J-1)
3 - TY2( J)*FF( I,J) + TY3( J)*FF( I,J+1) )
C
GO TO 100
C
2 CONTINUE
C
IF ( WPT( I,J,4) ) GO TO 8
C
C ----- Periodic right ---
C
A1 = 1.D0 / ( HX( I-2)*( HX( I-2) + HX( I-1) ) )
B1 = 1.D0 / ( HX( I-2)*HX( I-1) )
C1 = 1.D0 / ( HX( I-1)*( HX( I-2) + HX( I-1) ) )
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
WP1( I,J) = WP1( I,J) + 2.D0* HH( I,J)
1 * ( A1*FF( I-2,J) - B1*FF( I-1,J) + C1*FF( I,J)
2 + A2*FF( I,J-1) - B2*FF( I,J) + C2*FF( I,J+1) )
C
GO TO 100
C
C ----- Continuative right ---
C
8 CONTINUE
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
WP1( I,J) = WP1( I,J) +2.D0*HH( I,J)
1 *( A2*FF( I,J-1) - B2*FF( I,J) + C2*FF( I,J+1) )
C
GO TO 100
C
12 CONTINUE
C
IF ( WPT( I,J,4) ) GO TO 18
C
C ----- Periodic left ---
C
A1 = 1.D0 / ( HX( I)*( HX( I) + HX( I+1) ) )
B1 = 1.D0 / ( HX( I)*HX( I+1) )
C1 = 1.D0 / ( HX( I+1)*( HX( I) + HX( I+1) ) )
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
WP1( I,J) = WP1( I,J) +2.D0*HH( I,J)
1 *( A1*FF( I,J) - B1*FF( I+1,J) + C1*FF( I+2,J)
2 + A2*FF( I,J-1) - B2*FF( I,J) + C2*FF( I,J+1) )
C
GO TO 100
C
C ----- Continuative left ---
C
18 CONTINUE
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
WP1( I,J) = WP1( I,J) + 2.D0*HH( I,J)
1 *( A2*FF( I,J-1) - B2*FF( I,J) + C2*FF( I,J+1) )
C
GO TO 100
C
22 CONTINUE
C
IF ( WPT( I,J,4) ) GO TO 28
C
C ----- Periodic Upper ---
C
GO TO 100
C
C ----- Continuative Upper ---
C
28 CONTINUE
C
IF ( EPT( I-1,J,4).OR.EPT( I+1,J,3) ) GO TO 29
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) )
1 *( HX( I-1) + HX( I) + HX( I+1) )
C
HD = HX( I)*HX( I+1)
HE = HX( I) + HX( I+1)
C
A1 = -2.D0*( HX( I-1)*(2.D0*HX( I) + HX( I+1) )
1 - HX( I) *( HX( I) + HX( I+1) ) )/ HA
C
B1 = ( -A1*( 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 + 2.D0*( 2.D0*HX( I) + HX( I+1) ) )/ HB
C
D1 = ( A1*( HX( I-2) + HX( I-1) )
1 *( HX( I-2) + HX( I-1) + HX( I) + HX( I+1) )
2 + B1*HX( I-1)*( HX( I-1) + HX( I) + HX( I+1) )
3 - 2.D0 )/ HD
C
E1 = ( A1*( HX( I-2) + HX( I-1) ) + B1*HX( I-1) - D1*HX( I) )/ HE
C
C1 = - ( A1 + B1 + D1 + E1)
C
A2 = 2.D0 / ( HY( J-2)*( HY( J-2) + HY( J-1) ) )
C
B2 = 2.D0 / ( HY( J-2)*HY( J-1) )
C
C2 = 2.D0 / ( HY( J-1)*( HY( J-2) + HY( J-1) ) )
C
WP1( I,J) = WP1( I,J)
1 + HH( I,J)*( A1*FF( I-2,J) + B1*FF( I-1,J) + C1*FF( I,J)
C
2 + D1*FF( I+1,J) + E1*FF( I+2,J) + A2*FF( I,J-2)
3 - B2*FF( I,J-1) + C2*FF( I,J) )
C
GO TO 100
C
29 CONTINUE
C
A1 = 1.D0 / ( HX( I-1)*( HX( I-1) + HX( I) ) )
B1 = 1.D0 / ( HX( I-1)*HX( I) )
C1 = 1.D0 / ( HX( I) *( HX( I-1) + HX( I) ) )
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
WP1( I,J) = WP1( I,J)
1 +2.D0*HH( I,J)*( A1*FF( I-1,J) - B1*FF( I,J) + C1*FF( I+1,J) )
C
GO TO 100
C
32 CONTINUE
C
IF ( WPT( I,J,4) ) GO TO 38
C
C ----- Periodic Under ---
C
A1 = 1.D0 / ( HX( I-1)*( HX( I-1) + HX( I) ) )
B1 = 1.D0 / ( HX( I-1)*HX( I) )
C1 = 1.D0 / ( HX( I)*( HX( I-1) + HX( I) ) )
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
WP1( I,J) = WP1( I,J) +2.D0*HH( I,J)
1 *( A1*FF( I-1,J) - B1*FF( I,J) + C1*FF( I+1,J)
2 + A2*FF( I,J) - B2*FF( I,J+1) + C2*FF( I,J+2) )
C
GO TO 100
C
C ----- Continuative Under ---
C
38 CONTINUE
C
IF ( EPT( I-1,J,4) .OR. EPT( I+1,J,3) ) GO TO 39
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) )
1 *( HX( I-1) + HX( I) + HX( I+1) )
C
HD = HX( I)*HX( I+1)
C
HE = HX( I) + HX( I+1)
C
A1 = - 2.D0*( HX( I-1)*(2.D0 *HX( I) + HX( I+1) )
1 - HX( I)*( HX( I) + HX( I+1) ) )/ HA
C
B1 =( -A1*( 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 +2.D0*( 2.D0*HX( I) + HX( I+1) ) )/ HB
C
D1 =( A1*( HX( I-2) + HX( I-1) )
1 *( HX( I-2) + HX( I-1) + HX( I) + HX( I+1) )
C
2 + B1*HX( I-1)*( HX( I-1) + HX( I) + HX( I+1) ) - 2.D0 )/ HD
C
E1 = ( A1*( HX( I-2) + HX( I-1) ) + B1*HX( I-1) - D1*HX( I) )/ HE
C
C1 = - ( A1 + B1 + D1 + E1)
C
A2 = 2.D0 / ( HY( J)*( HY( J) + HY( J+1) ) )
B2 = 2.D0 / ( HY( J)*HY( J+1) )
C2 = 2.D0 / ( HY( J+1)*( HY( J) + HY( J+1) ) )
C
WP1( I,J) = WP1( I,J) + HH( I,J)
1 *( A1*FF( I-2,J) + B1*FF( I-1,J) + C1*FF( I,J)
C
2 + D1*FF( I+1,J) + E1*FF( I+2,J)
3 + A2*FF( I,J) - B2*FF( I,J+1) + C2*FF( I,J+2) )
C
GO TO 100
C
39 CONTINUE
C
A1 = 1.D0 / ( HX( I-1)*( HX( I-1) + HX( I) ) )
B1 = 1.D0 / ( HX( I-1)*HX( I) )
C1 = 1.D0 / ( HX( I) *( HX( I-1) + HX( I) ) )
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
WP1( I,J) = WP1( I,J) +2.D0*HH( I,J)
1 *( A1*FF( I-1,J) - B1*FF( I,J) + C1*FF( I+1,J) )
C
GO TO 100
C
C ----- Insulated wall ---
C
200 CONTINUE
C
IF ( I.EQ.1 .AND. J.EQ.1 ) GO TO 403
IF ( I.EQ.IX .AND. J.EQ.1 ) GO TO 404
IF ( I.EQ.1 .AND. J.EQ.IY) GO TO 405
IF ( I.EQ.IX .AND. J.EQ.IY) GO TO 406
C
IF ( EPT( I,J,3) ) GO TO 203
IF ( EPT( I,J,4) ) GO TO 204
IF ( EPT( I,J,1) ) GO TO 205
IF ( EPT( I,J,2) ) GO TO 206
C
403 WP1( I,J) = WP1( I,J) + HH( I,J)*( 2.D0*FF( I,J+1)
1 + ( 2.D0*FF( I+1,J) - CF*FF( I,J) )/ FS )
C
GO TO 100
C
404 WP1( I,J) = WP1( I,J) + HH( I,J)*( 2.D0*FF( I,J+1)
1 + ( 2.D0*FF( I-1,J) - CF*FF( I,J) )/ FS )
C
GO TO 100
C
405 WP1( I,J) = WP1( I,J) + HH( I,J)*( 2.D0*FF( I,J-1)
1 + ( 2.D0*FF( I+1,J) - CF*FF( I,J) )/ FS )
C
GO TO 100

406 WP1( I,J) = WP1( I,J) + HH( I,J)*( 2.D0*FF( I,J-1)
1 + ( 2.D0*FF( I-1,J) - CF*FF( I,J) )/ FS )
C
GO TO 100
C
203 WP1( I,J) = WP1( I,J) + HH( I,J)*(FF( I,J+1) +FF( I,J-1)
1 + ( FF( I+1,J) + FF( I-1,J) - CF*T( I,J) )/ FS )
GO TO 100
C
204 WP1( I,J) = WP1( I,J) + HH( I,J)*(FF( I,J+1) +FF( I,J-1)
1 + ( FF( I-1,J) + FF( I+1,J) - CF*FF( I,J) )/ FS )
GO TO 100
C
205 WP1( I,J) = WP1( I,J) + HH( I,J)*(FF( I,J-1) +FF( I,J-1)
1 + ( FF( I+1,J) + FF( I-1,J) - CF*FF( I,J) )/ FS )
C
GO TO 100
C
206 WP1( I,J) = WP1( I,J)
1 + HH( I,J) * ( FF(I,J+1) + FF(I,J+1)
2 + ( FF(I+1,J) + FF(I-1,J) - CF* FF(I,J) )/ FS )
C
GO TO 100
C
C ----- Perfectly conducting wall ---
C
101 GO TO 100
C
100 CONTINUE
C **********************************
C
500 CONTINUE
C ******************************************
C
700 CONTINUE
C *****************************************************
C
RETURN
END
C **********************************************************************
C
SUBROUTINE WFCA( WP1, PS1, WPT, HX, HY, RAS, RAF, IX, IY, IX1,
1 IY1 )
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION WP1( IX,IY), PS1( IX,IY), HX( IX1), HY( IY1),
1 RAF( IX1), RAS( IY1)
C
LOGICAL *1 WPT( IX,IY,4)
C
COMMON /UDAT/ AS, ITN
COMMON /CTRL/ MSH, IHS, WEI
C
C ----- No-slip surfaces : Updated ---
C
DO 100 I = 1, IX
DO 100 J = 1, IY
C
IF ( WPT( I,J,1) ) GO TO 100 ! Fluid
IF ( WPT( I,J,2) ) GO TO 50 ! Slip
IF ( WPT( I,J,3) ) GO TO 200 ! Moving
C
C ----- Attention ---
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 ----- Interior No-Slip Surfaces ---
C
IF ( .NOT. WPT( I+1,J, 1) ) GO TO 1 ! Boundary
WA = ( PS1( I,J) - PS1( I+1,J) )/ RAF( I)
C
13 IF ( .NOT. WPT( I-1,J, 1) ) GO TO 3 ! Boundary
WB = ( PS1( I,J) - PS1( I-1,J) )/ RAF( I-1)
C
12 IF ( .NOT. WPT( I,J+1, 1) ) GO TO 2 ! Boundary
WC = ( PS1( I,J) - PS1( I,J+1) )/ RAS( J)
C
14 IF ( .NOT. WPT( I,J-1, 1) ) GO TO 4 ! Boundary
C
WD = ( PS1( I,J) - PS1( I,J-1) )/ RAS( J-1)
C
114 WP1( I,J) = WA + WB + WC + WD
C
GO TO 100
C
C ----- Boundary ---
C
1 WA = ( PS1(I,J) - PS1(I-1,J) )/ RAF(I-1)
C
GO TO 13
C
3 WB = ( PS1(I,J) - PS1(I+1,J) )/ RAF(I)
C
GO TO 12
C
2 WC = ( PS1(I,J) - PS1(I,J-1) )/ RAS( J-1)
C
GO TO 14
C
4 WD = ( PS1(I,J) - PS1(I,J+1) )/ RAS( J)
C
GO TO 114
C
C ----- No-slip surfaces at edges ---
C
C ( Upper Boundary )
C
10 CONTINUE
C
C ----- Cavity Flow Problem ---
C
C ( Moving Boundary)
C **************************
C
WP1( I,J) = 3.D0*( PS1( I,J) - PS1( I,J-1) - HY( J-1) )
C ----- Attention --- **************
C
1 / RAS( J-1) - 0.5D0* WP1( I,J-1)
C
IF ( I.EQ.1.OR.I.EQ.IX) WP1(I,J) = 0.D0
C
GO TO 100
C
20 WP1( I,J) = 3.D0*( PS1( I,J) - PS1( I,J+1) )/ RAS( J)
1 - 0.5D0* WP1( I,J+1)
C
IF ( I.EQ.1.OR.I.EQ.IX) WP1(I,J) = 0.D0
C
GO TO 100
C
30 WP1( I,J) = 3.D0*( PS1( I,J) - PS1( I-1,J) )
1 / RAF( I-1) - 0.5D0* WP1( I-1,J)
GO TO 100
C
40 WP1( I,J) = 3.D0*( PS1( I,J) - PS1( I+1,J) )/ RAF( I)
1 - 0.5D0* WP1( I+1,J)
GO TO 100
C
50 WP1(I,J) = 0.D0
C
GO TO 100
C
200 CONTINUE
C
IF ( J.EQ.IY) GO TO 10
IF ( J.EQ.1 ) GO TO 20
C
100 CONTINUE
C
RETURN
END
C **********************************************************************
C
SUBROUTINE STND ( XX, YY, NDE, NNP, NEL, XZ, YZ, IX, IY )
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION XX( NNP), YY( NNP), XZ( IX), YZ( IY), NDE( NEL,4)
C
IF ( IX.EQ.IY ) GO TO 1000
C
WRITE(6,*) '*IX.NE.IY*STOP*'
WRITE(8,*) '*IX.NE.IY*STOP*'
C
STOP
C
1000 NN = 0
C
C ----------------------------
DO 100 J = 1, IY
DO 100 I = 1, IX
C
NN = NN + 1
XX( NN) = XZ( I)
YY( NN) = YZ( J)
C
100 CONTINUE
C -----------------------------
C
*** WRITE(6,*) ' * XX=', ( XX(I), I=1, NNP )
C
C ------------------------------------------------
DO 200 J = 1, IY-1
DO 200 I = 1, IX-1
C ------------------------------------------------
C
NUM = ( J-1)*( IX-1) + I
C
NDE( NUM,1) = ( J-1)*( IX) + I
NDE( NUM,2) = ( J-1)*( IX) + I + 1
NDE( NUM,3) = ( J-1)*( IX) + I + ( IX + 1)
NDE( NUM,4) = ( J-1)*( IX) + I + ( IX)
C
200 CONTINUE
C ------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE STMP ( IX, IY, DIM1, DIM2, NNP )
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION DIM1( IX,IY), DIM2( NNP)
C
C --------------------------------------
NN = 0
C
DO 100 J = 1, IY
DO 100 I = 1, IX
C
NN = NN + 1
DIM2( NN) = DIM1( I,IY-J+1)
C
100 CONTINUE
C --------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE PSCF ( IX, IY, HX, HY, IX1, IY1 )
C
IMPLICIT REAL*8(A-H,O-Z)
C
COMMON /CM1/ PA1( 51,51), PA2( 51,51), PA3( 51,51),
1 PB1( 51,51), PB2( 51,51), PB3( 51,51)
C
DIMENSION HX( IX1), HY( IY1)
C
C -------------------------
DO 100 I = 1, IX
DO 100 J = 1, IY
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 --------------------------
C
C -------------------------------------------------------------
DO 200 I = 2, IX-1
DO 200 J = 2, IY-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 WUCF ( HX, HY, UF1, UF2, UF3, UNA, UNB, UNC, UND,
1 UNE, UAP, UBP, UCP, UDP, UEP, UAM, UBM, UCM,
2 UDM, UEM, IX, IY, IX1, IY1 )
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION HX ( IX1), HY ( IY1), AA( 4), BB ( 4),
1 CC ( 4), DD ( 4), EE ( 4), UF1( IY),
2 UF2( IY), UF3( IY), UNA( IX,4), UNB( IX,4),
3 UNC( IX,4), UND( IX,4), UNE( IX,4), UAP( IX,4),
4 UBP( IX,4), UCP( IX,4), UDP( IX,4), UEP( IX,4),
5 UAM( IX,4), UBM( IX,4), UCM( IX,4), UDM( IX,4),
6 UEM( IX,4)
C
C --------------------------------------
DO 10 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
IF ( K.EQ.1) BB( K) = 1.D0
IF ( K.EQ.2) CC( K) = 2.D0
IF ( K.EQ.3) DD( K) = 6.D0
IF ( K.EQ.4) EE( K) = 24.D0
C
10 CONTINUE
C --------------------------------------
DO 20 I = 1, IX
DO 20 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
20 CONTINUE
C --------------------------------------
DO 30 J = 1, IY
C
UF1( J) = 0.D0
UF2( J) = 0.D0
UF3( J) = 0.D0
C
30 CONTINUE
C ----------------------------------------------------------
DO 40 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
40 CONTINUE
C ----------------------------------------------------------
DO 50 I = 3, IX-2
C
HA = HX(I-2)*( HX( I-2)+HX( I-1) )*( HX( I-2)+HX( I-1)+HX( I) )
1 *( 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)
C
HE = HX( I) + HX( I+1)
C
C ------------------------------------------------
DO 70 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) ) )
C
3 + DD( K)*( HX(I-1) - 2.D0*HX(I) - 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) )
C
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) )
C
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) + BB( K) )/ HE
C
UNC( I,K) = AA( K)
1 - ( UNA( I,K) + UNB( I,K) + UND( I,K) + UNE( I,K) )
C
70 CONTINUE
C ------------------------------------------------
C
50 CONTINUE
C ----------------------------------------------------------
DO 80 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) )
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 90 K = 1, 4

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) + ( HX( I-1)**2)
2 - HX( I-2)*HX( I) - 2.D0*HX( I-1)*HX( I) )
C
3 + DD(K)*( HX(I-2) +2.D0*HX(I-1) - HX(I) ) + EE(K) )/ HA
C
UBP( I,K) = ( - UAP( I,K)*( 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) )
3 + BB( K)*HX( I-1)*HX( I)
4 - CC( K)*( HX( I-1) - HX( I) ) - 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)
1 - ( UAP( I,K) + UBP( I,K) + UCP( I,K) + UEP( I,K) )
C
90 CONTINUE
C ------------------------------------------------
C
80 CONTINUE
C ----------------------------------------------------------
DO 100 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 110 K = 1, 4
C
UEM( I,K) = ( BB( K)*HX( I-1)*HX( I)*( HX( I) + HX( I+1) )
1 + CC( K)*(( HX( I)**2) + HX( I)*HX( I+1)
2 - HX( I-1)*HX( I+1) - 2.D0*HX( I-1)*HX( I) )
C
3 + DD( K)*( HX(I-1) - 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) + HX( I+2) )
C
3 - BB(K)*HX(I-1)*HX(I) - 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) )*( HX( I-1)+HX( I)+HX( I+1) )
C
3 + 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)
1 - ( UAM( I,K) + UCM( I,K) + UDM( I,K) + UEM( I,K) )
C
110 CONTINUE
C ------------------------------------------------
C
100 CONTINUE
C ----------------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE WUCL ( WP1, WM1, PS1, WPT, EPT, RUF, HX, HY,
1 UF1, UF2, UF3, UNA, UNB, UNC, UND, UNE,
2 UAP, UBP, UCP, UDP, UEP, UAM, UBM, UCM,
3 UDM, UEM, IX, IY, IX1, IY1 )
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION WP1( IX,IY), PS1( IX,IY), WM1( IX,IY), RUF( IX,IY),
1 HX ( IX1), HY ( IY1),
2 UF1( IY), UF2( IY), UF3( IY),
5 UNA( IX,4), UNB( IX,4), UNC( IX,4), UND( IX,4),
6 UNE( IX,4), UAP( IX,4), UBP( IX,4), UCP( IX,4),
7 UDP( IX,4), UEP( IX,4), UAM( IX,4), UBM( IX,4),
8 UCM( IX,4), UDM( IX,4), UEM( IX,4)
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
COMMON /UDAT/ AS, ITN
COMMON /CTRL/ MSH, IHS, WEI
COMMON /IDAT/ DT, LOPE, REN, LOP, NPL
C
C ------------------------------------------------
DO 300 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
IF ( K.EQ.1) BB( K) = 1.D0
IF ( K.EQ.2) CC( K) = 2.D0
IF ( K.EQ.3) DD( K) = 6.D0
IF ( K.EQ.4) EE( K) = 24.D0
C
300 CONTINUE
C ------------------------------------------------
C
C ----------------------------------------------------------
DO 700 I = 1, IX
DO 700 J = 1, IY
C ----------------------------------------------------------
C
IF ( EPT( I,J,3) ) GO TO 14
IF ( .NOT.EPT( I,J,1) ) GO TO 5
IF ( .NOT.WPT( I,J,4) ) GO TO 15
C
RUF( I,J) = 0.D0
GO TO 700
C
5 IF ( .NOT.EPT( I,J,2) ) GO TO 10
IF ( .NOT.WPT( I,J,4) ) GO TO 20
RUF( I,J) = 0.D0
GO TO 700
10 CONTINUE
RUF(I,J) = UF1( J)* PS1(I,J-1) + UF2( J)* PS1(I,J)
1 + UF3( J)* PS1(I,J+1)
C
14 IF ( .NOT.EPT( I,J,1) ) GO TO 112
15 IF ( .NOT.WPT( I,J,2) ) GO TO 16
C
RUF( I,J) = RUF( I,J-1)
GO TO 115
C
16 IF ( .NOT.WPT( I,J,3)) GO TO 17 ! Not Moving
UU = 1.D0
RUF( I,J) = UU
GO TO 115
C
17 RUF( I,J) = 0.D0
GO TO 115
C
112 IF ( .NOT.EPT( I,J,2)) GO TO 113
20 IF ( .NOT.WPT( I,J,2)) GO TO 21
RUF( I,J) = RUF( I,J+1)
GO TO 115
C
21 IF ( .NOT.WPT( I,J,3) ) GO TO 22
UL = 0.D0
RUF( I,J) = UL ! Not moving
GO TO 115
C
22 RUF( I,J) = 0.D0
GO TO 115
C
113 IF ( .NOT.EPT( I,J,3) ) GO TO 114
IF ( .NOT.WPT( I,J,2) ) GO TO 26
RUF( I,J) = 0.D0
GO TO 115
C
26 IF ( .NOT.WPT( I,J,3) ) GO TO 27
RUF( I,J) = 0.D0
GO TO 115
C
27 IF ( .NOT.WPT( I,J,4) ) GO TO 28
C
U1 = - HY( J)/ ( HY( J-1)*( HY( J-1) + HY( J) ) )
U2 = ( HY( J) - HY( J-1) )/ ( HY( J-1)*HY( J) )
U3 = HY( J-1)/ ( HY( J)*( HY( J-1) + HY( J) ) )
C
RUF( I,J) = U1*PS1( I,J-1) + U2*PS1( I,J) + U3*PS1( I,J+1)
GO TO 115
C
28 RUF( I,J) = 0.D0
GO TO 115
C
114 IF ( .NOT.EPT( I,J,4) ) GO TO 115
IF ( .NOT.WPT( I,J,2) ) GO TO 31
C
RUF( I,J) = 0.D0
GO TO 115
C
31 IF ( .NOT.WPT( I,J,3) ) GO TO 32
RUF( I,J) = 0.D0
GO TO 115
C
32 IF ( .NOT.WPT( I,J,4) ) GO TO 33
RUF( I,J) = RUF( I,J)
GO TO 115
C
33 CONTINUE
C
115 CONTINUE
C
UT = RUF( 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
C
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)
1 + C1* WM1( I+1,J) )
2 + ( UT**2)*( A2* WM1( I-1,J) - B2* WM1( I,J)
3 + C2* WM1( I+1,J) )
C
80 CONTINUE
C
IF ( RUF( I,J).LT.0 ) GO TO 90
IF ( EPT( I,J,4) ) GO TO 333
IF ( EPT( I,J,3) ) GO TO 130
IF ( EPT( I-1,J,4) ) GO TO 110
IF ( EPT( I-2,J,4) ) GO TO 120
C
UPWN = - 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
333 WP1( I,J) = WP1( I,J) - UT* WM1( I,J)/ HX( I)
GO TO 700
C
110 UPWN = - UT*( - WM1( I-1,J) + WM1( I,J) )/ HX( I-1)
GO TO 150
C
120 CONTINUE
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) )
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) )
C
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
UPWN = - 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
4 -( UT**3)*( UPA( 3)* WM1( I-2,J) + UPB( 3)* WM1( I-1,J)
5 + 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)*( 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
UPWN = - UT *( A1* WM1( I-2,J) - B1* WM1( I-1,J)
1 + C1* WM1( I,J) )
2 + ( UT**2)*( A2* WM1( I-2,J) - B2* WM1( I-1,J)
3 + C2* WM1( I,J) )
C
WP1( I,J) = WP1( I,J) + UPWN
C
GO TO 700
C
90 IF ( EPT( I,J,4) ) GO TO 170
IF ( EPT( I,J,3) ) GO TO 200
IF ( EPT( I+1,J,3) ) GO TO 190
IF ( EPT( I+2,J,3) ) GO TO 180
C
UPWN = - 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
C
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
UPWN = - UT*( - A1* WM1( I,J) + B1* WM1( I+1,J)
1 - C1* WM1( I+2,J) )
2 +( UT** 2)*( A2* WM1( I,J) - B2* WM1( I+1,J)
3 + C2* WM1( I+2,J) )
C
WP1( I,J) = WP1( I,J) + UPWN
GO TO 700
C
180 CONTINUE
C
DO 230 K = 1, 3
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
UPD( K) = ( - BB( K)*HX(I-1)*HX(I) - CC( K)*( HX(I) - HX(I-1) )
1 + DD( K) )/ HD
C
UPC( K)=( -UPD(K)*( HX(I)+HX(I+1) )*( HX(I-1)+HX(I)+HX(I+1) )
C
1 + BB( K)*HX( I-1) + CC( K) )/ HC
C
UPA( K) =( UPD(K)*( HX(I)+HX(I+1) ) + UPC(K)*HX(I) - BB(K) )/ HA
C
UPB( K) = AA( K) - ( UPA(K) + UPC(K) + UPD(K) )
C
230 CONTINUE
C
UPWN = - UT*( UPA( 1)* WM1( I-1,J) + UPB( 1)* WM1( I,J)
1 + UPC( 1)* WM1( I+1,J) + UPD( 1)* WM1( I+2,J) )
C
2 + ( UT**2)*( UPA( 2)* WM1( I-1,J) + UPB( 2)* WM1( I,J)
3 + UPC( 2)* WM1( I+1,J) + UPD( 2)* WM1( I+2,J) )/ 2.D0
C
4 - ( UT**3)*( UPA( 3)* WM1( I-1,J) + UPB( 3)* WM1( I,J)
5 + UPC( 3)* WM1( I+1,J) + UPD( 3)* WM1( I+2,J) )/ 6.D0
GO TO 220
C
190 UPWN = - UT*( - WM1( I,J) + WM1( I+1,J) )/ HX( I)
GO TO 220
C
200 WP1( I,J) = WP1( I,J) + UT* WM1( I,J)/ HX( I-1)
GO TO 700
C
220 CONTINUE
C
150 CONTINUE
C
WP1( I,J) = WP1( I,J) + ( 1.D0-WEI )* CENT + WEI* UPWN
C ----------------------------------------------------------
C
700 CONTINUE
C ----------------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE WVCF ( HX, HY, VF1, VF2, VF3, VNA, VNB, VNC,
1 VND, VNE, VAP, VBP, VCP, VDP, VEP, VAM,
2 VBM, VCM, VDM, VEM, IX, IY, IX1, IY1 )
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION HX ( IX1), HY ( IY1), AA ( 4), BB ( 4),
1 CC ( 4), DD ( 4), EE ( 4), VF1( IX),
2 VF2( IX), VF3( IX), VNA( IY,4), VNB( IY,4),
3 VNC( IY,4), VND( IY,4), VNE( IY,4), VAP( IY,4),
4 VBP( IY,4), VCP( IY,4), VDP( IY,4), VEP( IY,4),
5 VAM( IY,4), VBM( IY,4), VCM( IY,4), VDM( IY,4),
6 VEM( IY,4)
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
IF ( K.EQ.1) BB( K) = 1.D0
IF ( K.EQ.2) CC( K) = 2.D0
IF ( K.EQ.3) DD( K) = 6.D0
IF ( K.EQ.4) EE( K) = 24.D0
C
100 CONTINUE
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 --------------------------------------
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) )
1 *( HY( J-1) + HY( J) + HY( J+1) )
C
HD = HY( J)*HY( J+1)
C
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) )
2 + VNB( J,K)*HY( J-1)*( HY( J-1) + HY( J) + HY( J+1) )
C
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)
2 - VND( J,K)*HY( J) + BB( K) )/ HE
C
VNC( J,K)= AA( K) - ( VNA( J,K) +VNB( J,K) +VND( J,K) +VNE( J,K) )
C
700 CONTINUE
C --------------------------------
C
500 CONTINUE
C ------------------------------------------------
C
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) )
1 *( HY( J-2) + HY( J-1) + HY( J) )
C
HC = HY( J-1)*( HY( J-1) + HY( J) )
C
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)
2 + ( HY( J-1)**2) - HY( J-2)*HY( J)
3 - 2.D0*HY( J-1)*HY( J) )
C
4 + DD( K)*( HY( J-2) +2.D0*HY( J-1) - HY( J) )
5 + 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) )
5 - DD( K) )/ HB
C
VCP( J,K) = ( - VAP( J,K)*( HY( J-3) + HY( J-2) + HY( J-1) )
1 *( HY( J-3) + HY( J-2) + HY( J-1) + HY( J) )
C
2 - VBP( J,K)*( HY( J-2) + HY( J-1) )
3 *( HY( J-2) + HY( J-1) + HY( J) )
C
4 - 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) )
C
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 920 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) )
HC = HY( J)*( HY( J-1) + HY( J) )
HA = HY( J-1)
C
C -------------------------------
DO 940 K = 1, 4
C
VEM( J,K) = ( BB( K)*HY( J-1)*HY( J)*( HY( J) + HY( J+1) )
1 + CC( K)*( ( HY( J)**2) + HY( J)*HY( J+1)
2 - HY( J-1)*HY( J+1) - 2.D0*HY( J-1)*HY( J) )
C
3 + DD(K)*( HY( J-1) - 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) )
C
2 - VDM( J,K)* ( HY( J) + HY( J+1) )
3 * ( HY( J-1) + HY( J) + HY( J+1) )
C
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) - ( VAM( J,K) + VCM( J,K) + VDM( J,K)
1 + VEM( J,K) )
C
940 CONTINUE
C -----------------------------
C
920 CONTINUE
C --------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE WVCL ( WP1, WM1, PS1, WPT, EPT, RVF, HX, HY,
1 VF1, VF2, VF3, UNA, UNB, UNC, UND, UNE,
2 UAP, UBP, UCP, UDP, UEP, UAM, UBM, UCM,
3 UDM, UEM, IX, IY, IX1, IY1 )
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION WP1( IX,IY), PS1( IX,IY), WM1( IX,IY), RVF( IX,IY),
1 HX ( IX1), HY ( IY1), AA ( 4), BB ( 4),
2 CC ( 4), DD ( 4), EE ( 4), UPA( 4),
3 UPB( 4), UPC( 4), UPD( 4), VF1( IX),
4 VF2( IX), VF3( IX), UNA( IY,4), UNB( IY,4),
C
5 UNC( IY,4), UND( IY,4), UNE( IY,4),
6 UAP( IY,4), UBP( IY,4), UCP( IY,4), UDP( IY,4),
7 UEP( IY,4), UAM( IY,4), UBM( IY,4), UCM( IY,4),
8 UDM( IY,4), UEM( IY,4)
C
LOGICAL*1 WPT( IX,IY,4), EPT( IX,IY,4)
C
COMMON /UDAT/ AS, ITN
COMMON /CTRL/ MSH, IHS, WEI
COMMON /IDAT/ DT, LOPE, REN, LOP, NPL
C
VR = 0.D0
VL = 0.D0
C
C ------------------------------------------------
DO 300 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
IF ( K.EQ.1) BB( K) = 1.D0
IF ( K.EQ.2) CC( K) = 2.D0
IF ( K.EQ.3) DD( K) = 6.D0
IF ( K.EQ.4) EE( K) = 24.D0
C
300 CONTINUE
C ------------------------------------------------
C
C ----------------------------------------------------------
DO 700 J = 1, IY
DO 700 I = 1, IX
C ----------------------------------------------------------
C
IF ( EPT( I,J,1) ) GO TO 14
IF ( .NOT.EPT( I,J,3) ) GO TO 5
C
IF ( WPT( I,J,4) ) GO TO 700
C
GO TO 25
C
5 IF ( .NOT.EPT( I,J,4) ) GO TO 10
IF ( WPT( I,J,4) ) GO TO 700
C
GO TO 30
C
10 CONTINUE
C
RVF( I,J)
1 = - ( VF1(I)* PS1(I-1,J) + VF2(I)* PS1(I,J)
2 + VF3(I)* PS1(I+1,J) )
C
14 IF ( .NOT.EPT( I,J,1 ) ) GO TO 112
IF ( .NOT.WPT( I,J,2) ) GO TO 16
C
RVF( I,J) = 0.D0
GO TO 115
C
16 IF ( .NOT.WPT( I,J,3) ) GO TO 17
C
RVF( I,J) = 0.D0
GO TO 115
C
17 IF ( .NOT.WPT( I,J,4) ) GO TO 18
C
V1 = - HX( I)/ ( HX( I-1)*( HX( I-1) + HX( I) ) )
V2 = ( HX( I) - HX( I-1) )/ ( HX( I-1)*HX( I) )
V3 = HX( I-1)/ ( HX( I)*( HX( I-1) + HX( I) ) )
C
RVF( I,J) = - ( V1*PS1( I-1,J) + V2*PS1( I,J) + V3*PS1( I+1,J) )
GO TO 115
C
18 RVF( I,J) = 0.D0
GO TO 115
C
112 IF ( .NOT.EPT( I,J,2) ) GO TO 113
IF ( .NOT.WPT( I,J,2) ) GO TO 21
RVF( I,J) = 0.D0
GO TO 115
C
21 IF ( .NOT.WPT( I,J,3) ) GO TO 22
C
RVF( I,J) = 0.D0
GO TO 115
C
22 IF ( .NOT.WPT( I,J,4) ) GO TO 23
RVF( I,J) = RVF( I,J)
GO TO 115
C
23 RVF( I,J) = 0.D0
GO TO 115
C
113 IF ( .NOT.EPT( I,J,3) ) GO TO 114
25 IF ( .NOT.WPT( I,J,2) ) GO TO 26
RVF( I,J) = RVF( I-1,J)
GO TO 115
C
26 IF ( .NOT.WPT( I,J,3) ) GO TO 27
C
RVF( I,J) = VR
GO TO 115
C
27 RVF( I,J) = 0.D0
GO TO 115
C
114 IF ( .NOT.EPT( I,J,4) ) GO TO 115
30 IF ( .NOT.WPT( I,J,2) ) GO TO 31
RVF( I,J) = RVF( I+1,J)
GO TO 115
C
31 IF ( .NOT.WPT( I,J,3) ) GO TO 32
RVF( I,J) = VL
GO TO 115
C
32 RVF( I,J) = 0.D0
C
115 CONTINUE
C ---------------------------
C
VT = RVF( I,J)* DT
C
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*( UNA( J,1)* WM1( I,J-2) + UNB( J,1)* WM1( I,J-1)
1 + UNC( J,1)* WM1( I,J ) + UND( J,1)* WM1( I,J+1)
2 + UNE( J,1)* WM1( I,J+2) )
C
3 + ( VT**2)* ( UNA( J,2)* WM1( I,J-2) + UNB( J,2)* WM1( I,J-1)
4 + UNC( J,2)* WM1( I,J ) + UND( J,2)* WM1( I,J+1)
5 + UNE( J,2)* WM1( I,J+2) )/ 2.D0
C
6 - ( VT**3)* ( UNA( J,3)* WM1( I,J-2) + UNB( J,3)* WM1( I,J-1)
7 + UNC( J,3)* WM1( I,J ) + UND( J,3)* WM1( I,J+1)
8 + UNE( J,3)* WM1( I,J+2) )/ 6.D0
C
9 + ( VT**4)* ( UNA( J,4)* WM1( I,J-2) + UNB( J,4)* WM1( I,J-1)
A + UNC( J,4)* WM1( I,J ) + UND( J,4)* WM1( I,J+1)
B + UNE( J,4)* WM1( I,J+2) )/ 24.D0
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)
1 + C1* WM1( I,J+1) )
2 + ( VT**2)*( A2* WM1( I,J-1) - B2* WM1( I,J) + C2* WM1( I,J+1) )
C
80 CONTINUE
C
IF ( RVF( I,J).LT.0 ) GO TO 90
IF ( EPT( I,J,1) ) GO TO 130
IF ( EPT( I,J,2) ) GO TO 100
IF ( EPT( I,J-1,2) ) GO TO 110
IF ( EPT( I,J-2,2) ) GO TO 120
C
UPWN = - VT*( UAP( J,1)* WM1( I,J-3) + UBP( J,1)* WM1( I,J-2)
1 + UCP( J,1)* WM1( I,J-1 ) + UDP( J,1)* WM1( I,J)
2 + UEP( J,1)* WM1( I,J+1 ) )
C
3 + ( VT**2)*( UAP( J,2)* WM1( I,J-3) + UBP( J,2)* WM1( I,J-2)
4 + UCP( J,2)* WM1( I,J-1 ) + UDP( J,2)* WM1( I,J)
5 + UEP( J,2)* WM1( I,J+1 ) )/ 2.D0
C
6 - ( VT**3)*( UAP( J,3)* WM1( I,J-3) + UBP( J,3)* WM1( I,J-2)
7 + UCP( J,3)* WM1( I,J-1 ) + UDP( J,3)* WM1( I,J)
8 + UEP( J,3)* WM1( I,J+1 ) )/ 6.D0
C
9 + ( VT**4)*( UAP( J,4)* WM1( I,J-3) + UBP( J,4)* WM1( I,J-2)
A + UCP( J,4)* WM1( I,J-1 ) + UDP( J,4)* WM1( I,J)
B + UEP( J,4)* WM1( I,J+1 ) )/ 24.D0
GO TO 150
C
100 WP1( I,J) = WP1( I,J) - VT* WM1( I,J)/ HY( J)
GO TO 700
C
110 UPWN = - VT*( - WM1( I,J-1) + WM1( I,J) )/ HY( J-1)
GO TO 150
C
120 CONTINUE
C
C ----------------------------------
DO 160 K = 1, 3
C
HA = HY( J-2)*( HY( J-2)+HY( J-1) )*( HY( J-2)+HY( J-1)+HY( J) )
HB = HY( J-1)*( HY( J-1) + HY( J) )
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) )
1 + UPB( K)*HY( J-1) + BB( K) )/ HD
C
UPC( K) = AA( K) - ( UPA( K)+UPB( K)+UPD( K) )
C
160 CONTINUE
C --------------------------------
C
UPWN = - VT*( UPA( 1)* WM1( I,J-2) + UPB( 1)* WM1( I,J-1)
1 + UPC( 1)* WM1( I,J) + UPD( 1)* WM1( I,J+1) )
C
2 + ( VT**2)*( UPA( 2)* WM1( I,J-2) + UPB( 2)* WM1( I,J-1)
3 + 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
UPWN = - 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) +C2* WM1( I,J) )
C
WP1( I,J) = WP1( I,J) + UPWN
C
GO TO 700
C
90 IF ( EPT( I,J,2) ) GO TO 170
IF ( EPT( I,J,1) ) GO TO 200
IF ( EPT( I,J+1,1) ) GO TO 190
IF ( EPT( I,J+2,1) ) GO TO 180
C
UPWN = - VT*( UAM( J,1)* WM1( I,J-1) + UBM( J,1)* WM1( I,J)
1 + UCM( J,1)* WM1( I,J+1) + UDM( J,1)* WM1( I,J+2)
2 + UEM( J,1)* WM1( I,J+3) )
C
3 + ( VT**2)*( UAM( J,2)* WM1( I,J-1) + UBM( J,2)* WM1( I,J)
4 + UCM( J,2)* WM1( I,J+1) + UDM( J,2)* WM1( I,J+2)
5 + UEM( J,2)* WM1( I,J+3) )/ 2.D0
C
6 - ( VT**3)*( UAM( J,3)* WM1( I,J-1) + UBM( J,3)* WM1( I,J)
7 + UCM( J,3)* WM1( I,J+1) + UDM( J,3)* WM1( I,J+2)
8 + UEM( J,3)* WM1( I,J+3) )/ 6.D0
C
9 + ( VT**4)*( UAM( J,4)* WM1( I,J-1) + UBM( J,4)* WM1( I,J)
A + UCM( J,4)* WM1( I,J+1) + UDM( J,4)* WM1( I,J+2)
B + UEM( J,4)* WM1( I,J+3) )/ 24.D0
C
GO TO 220
C
170 CONTINUE
C
A1 = ( 2.D0*HY( J)+HY( J+1) )/ ( HY( J)*( HY( J)+HY( J+1) ) )
C
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
UPWN = - VT*( - A1* WM1( I,J) + B1* WM1( I,J+1)
1 - C1* WM1( I,J+2) )
2 + ( VT**2)*( A2* WM1( I,J) - B2* WM1( I,J+1) + C2* WM1( I,J+2) )
C
WP1( I,J) = WP1( I,J) + UPWN
GO TO 700
C
180 CONTINUE
C
C -----------------------------
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 ---------------------------
C
UPWN = -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
C
GO TO 220
C
190 UPWN = - VT*( - WM1( I,J) + WM1( I,J+1) )/ HY( J)
C
GO TO 220
C
200 WP1( I,J) = WP1( I,J) + VT* WM1( I,J)/ HY( J-1)
C
GO TO 700
C
220 CONTINUE
C
150 CONTINUE
C
WP1( I,J) = WP1( I,J) + ( 1.D0 - WEI )* CENT + WEI* UPWN
C
700 CONTINUE
C ----------------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE TMCOF ( HX, HY, IX, IY, IX1, IY1 )
C
IMPLICIT REAL*8(A-H,O-Z)
C
COMMON /CM2/ TM1( 51), TM2( 51), TM3( 51), TY1( 51), TY2( 51),
1 TY3( 51)
C
DIMENSION HX( IX1), HY( IY1)
C
C --------------------------------------
DO 100 I = 1, IX
C
TM1( I) = 0.D0
TM2( I) = 0.D0
TM3( 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
TM1( I) = 1.D0 / ( HX( I-1)*( HX( I-1) + HX( I) ) )
TM2( I) = 1.D0 / ( HX( I-1)* HX( I) )
TM3( 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 **********************************************************************