戻る

C ********************************************************* C
C C
C H F 2 D - FDM ( SV ) C
C C
C 2D Thermal-Fluid Analysis by the 4th order FDM C
C C
C ( Meshes of Equal length & Unequal length ) C
C C
C ( Version 3.9 ) C
C C
C Copyright : Yasuhiro MATSUDA C
C C
C ********************************************************* C
C
PROGRAM MAIN
C
IMPLICIT REAL*8 (A-H,O-Z)
C
INTEGER X, Y
C
C ***** DATA ( 1) *** Attention ( Max.: 101 X 101 ) ********************
C
** PARAMETER ( X = 21, Y = 21,
C
** PARAMETER ( X = 41, Y = 41,
C
PARAMETER ( X = 51, Y = 51,
C
** PARAMETER ( X = 61, Y = 61,
** PARAMETER ( X = 81, Y = 81,
** PARAMETER ( X = 101, Y = 101,
C
1 IX1 = X-1, IY1 = Y-1, IXY = X* Y, IXY1 = IX1* IY1 )
C
C **********************************************************************
C
C KTF = 1 : Thermal Fluid Analysis ( Natural Convection Problem )
C = 0 : Fluid Analysis ( Cavity Flow Problem )
C
C MSH = 1 : Mesh of Equal length = 2 : Mesh of Unequal length
C
C IHS = 1 : With Correction
C
C WEI : Weighting Parameter
C
C ---------------------------------------------------------------------
C
LOGICAL*1 WPT( X,Y,4), PPT( X,Y,4), TPT( X,Y,4), EPT( X,Y,4)
C
DIMENSION PS1( X,Y), PS ( X,Y), WP1( X,Y), WP( X,Y),
1 TP1( X,Y), T ( X,Y)
C
DIMENSION U ( X,Y), V ( X,Y), UB ( X,Y), VB ( X,Y),
1 HH ( X,Y), SKN( X,Y), GG ( X,Y), SIM( X),
2 ANU( X), RDP( X,Y), BX ( X,Y), BY ( X,Y),
3 XZ ( X), YZ ( Y), HX ( IX1), HY ( IY1)
C
DIMENSION UF1( Y), UF2( Y), UF3( Y), UNA( X,4),
1 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),
3 UEP( X,4), UAM( X,4), UBM( X,4), UCM( X,4),
4 UDM( X,4), UEM( X,4)
C
DIMENSION VF1( X), VF2( X), VF3( X), VNA( Y,4),
1 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),
3 VEP( Y,4), VAM( Y,4), VBM( Y,4), VCM( Y,4),
4 VDM( Y,4), VEM( Y,4)
C
DIMENSION RAF( IX1), RAS( IY1), PE ( IXY), TE( IXY),
1 XX ( IXY), YY ( IXY), NDE( IXY1,4)
C
COMMON /CTR/ KTF, MSH, IHS, WEI
COMMON /IDT/ DT, LPM, PRN, CTP, RLN, REN, LOP, NPL, NPR
C
C ----- Attention --- X & Y < = 50 x 50 ---
C
COMMON /CM1/ PA1( 101,101), PA2( 101,101), PA3( 101,101),
1 PB1( 101,101), PB2( 101,101), PB3( 101,101)
C
COMMON /CM2/ TM1( 101), TM2( 101), TM3( 101),
1 TY1( 101), TY2( 101), TY3( 101)
C
C ***** DATA ( 2) ******************************************************
C
* OPEN ( 10, FILE='F2D_FDM_20.DAT')
C ----------------------------------------------------------
*
* OPEN ( 7, FILE='F2D_FDM_E20_RE100_FIG.DAT' ) ! DT = 0.04
* OPEN ( 8, FILE='F2D_FDM_E20_RE100_OUT.DAT' )
C
* OPEN ( 7, FILE='F2D_FDM_UE20_RE100_FIG.DAT' ) ! DT = 0.009
* OPEN ( 8, FILE='F2D_FDM_UE20_RE100_OUT.DAT' )
C
C ----------------------------------------------------------
C
* OPEN ( 10, FILE='HF2D_FDM_20.DAT')
C ----------------------------------------------------------
C
* OPEN ( 7, FILE='HF2D_FDM_UE20_RA10_5_FIG.DAT' )
* OPEN ( 8, FILE='HF2D_FDM_UE20_RA10_5_OUT.DAT' )
C
* OPEN ( 7, FILE='HF2D_FDM_E20_RA10_6_FIG.DAT' )
* OPEN ( 8, FILE='HF2D_FDM_E20_RA10_6_OUT.DAT' )
C
C --------------------------------------------------------------------
C
OPEN ( 10, FILE='F2D_FDM_50.DAT')
C ----------------------------------------------------------
C
* OPEN ( 7, FILE='F2D_FDM_E50_RE1000_FIG.DAT' ) ! DT = 0.01
* OPEN ( 8, FILE='F2D_FDM_E50_RE1000_OUT.DAT' )
C
OPEN ( 7, FILE='F2D_FDM_UE50_RE1000_FIG.DAT' ) ! DT = 0.008
OPEN ( 8, FILE='F2D_FDM_UE50_RE1000_OUT.DAT' )
C
C ----------------------------------------------------------
C
* OPEN ( 10, FILE='HF2D_FDM_50.DAT')
C ----------------------------------------------------------
C
* OPEN ( 7, FILE='HF2D_FDM_E50_RA10_6_FIG.DAT' )
* OPEN ( 8, FILE='HF2D_FDM_E50_RA10_6_OUT.DAT' )
C
* OPEN ( 7, FILE='HF2D_FDM_UE50_RA10_5_FIG.DAT' )
* OPEN ( 8, FILE='HF2D_FDM_UE50_RA10_5_OUT.DAT' )
C
* OPEN ( 7, FILE='HF2D_FDM_UE50_RA10_6_FIG.DAT' )
* OPEN ( 8, FILE='HF2D_FDM_UE50_RA10_6_OUT.DAT' )
C
* OPEN ( 7, FILE='W05_HF2D_FDM_UE50_RA10_7_FIG.DAT' )
* OPEN ( 8, FILE='W05_HF2D_FDM_UE50_RA10_7_OUT.DAT' )
C
* OPEN ( 7, FILE='W1_HF2D_FDM_UE50_RA10_7_FIG.DAT' )
* OPEN ( 8, FILE='W1_HF2D_FDM_UE50_RA10_7_OUT.DAT' )
C
C --------------------------------------------------------------------
C
* OPEN ( 10, FILE='HF2D_FDM_100.DAT')
C ----------------------------------------------------------
C
* OPEN ( 7, FILE='HF2D_FDM_UE100_RA10_8_FIG.DAT' )
* OPEN ( 8, FILE='HF2D_FDM_UE100_RA10_8_OUT.DAT' )
C
C **********************************************************************

CALL PRPS ( PPT, WPT, TPT, EPT, PS1, WP1, TP1, T, PS,
1 WP, U, V, UB, VB, SKN, GG, HH, BX,
2 BY, XZ, YZ, SIM, ANU, RAF, RAS, HX, HY,
3 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 --------------------------------------------------------------------
999 LOP = LOP + 1
C
ACP = 1.7D0
C
CALL PSCL ( PS1, PS, WP1, TP1, PPT, EPT, ACP, HX, HY,
1 RDP, X, Y, IX1, IY1 )
C
CALL VARC ( PS1, PS, VRP, X, Y )
C
C ----- Vorticity on Boundary ---
C
CALL WFCA ( WP1, PS1, WPT, HX, HY, RAS, RAF, X, Y,
1 IX1, IY1 )
C
IF ( KTF.NE.1) GO TO 150
C
C ----- Temperature comp. ----
C
CALL VARC ( TP1, T, VRT, X, Y )
C
CALL WUCL ( TP1, T, PS1, TPT, 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, 1 )
C
CALL WVCL ( TP1, T, PS1, TPT, 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, 1 )
C
CALL TMCL ( TP1, T, WP1, TPT, EPT, U, V, HH,
1 HX, HY, BX, BY, X, Y, IX1, IY1,
2 1 )
C
CALL CLNU ( SKN, GG, U, V, HX, HY, SIM, ANU,
1 TP1, XZ, YZ, X, Y, IX1, IY1 )
C
150 CONTINUE
C
CALL VARC ( WP1, WP, VRW, X, Y )
C
CALL WUCL ( WP1, WP, 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, 2 )
C
CALL WVCL ( WP1, WP, 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, 2 )
C
CALL TMCL ( WP1, WP, TP1, WPT, EPT, U, V, HH,
1 HX, HY, BX, BY, X, Y, IX1, IY1,
2 2 )
C
CALL CHCK ( U, V, UB, VB, VRW, VRT, VRP, X, Y, *900 )
C
GO TO 999
C
900 CONTINUE
C
CALL POUT ( TE, XX, YY, NDE, PE, PS, T, WP,
1 UB, VB, PS1, U, V, XZ, YZ, WP1,
2 TP1, SKN, GG, SIM, ANU, HX, HY, JU,
3 IV, X, Y, IXY, IXY1, IX1, IY1 )
C
CLOSE ( 7)
CLOSE ( 10)
C
STOP
END
C *********************************************************************
C
SUBROUTINE BCUV ( U, V, IX, IY )
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION U( IX,IY), V( IX, IY)
C
C ----------------------------
DO 100 I = 1, IX
C
U( I, 1) = 0.D0
V( I, 1) = 0.D0
U( I,IY) = 0.D0
V( I,IY) = 0.D0
C
100 CONTINUE
C ----------------------------
DO 200 J = 1, IY
C
U( 1, J) = 0.D0
V( 1, J) = 0.D0
U(IX, J) = 0.D0
V(IX, J) = 0.D0
C
200 CONTINUE
C ----------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE CHCK ( U, V, UB, VB, VRW, VRT, VRP, IX, IY, * )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
COMMON /UDT/ DY2, ITN
COMMON /ABC/ DNU, RN0, RN1
COMMON /CTR/ KTF, MSH, IHS, WEI
COMMON /IDT/ DT, LPM, PRN, CTP, RLN, REN, LOP, NPL, NPR
COMMON /DEF/ VAW( 250000), VAT( 250000), VVP( 250000),
1 VRN( 250000)
C
DIMENSION U( IX,IY), V( IX,IY), UB( IX,IY), VB( IX,IY)
C
KKK = 0
IF ( LOP.EQ.1) NPL = 0
C
IF ( MOD( LOP,NPR).EQ.0) THEN
NPL = NPL + 1
VAW( NPL) = VRW
VVP( NPL) = VRP
END IF
C
C ----------------------------------------------------------
IF ( KTF.EQ.1) THEN
C
C --------------------------------------------------
IF ( MOD( LOP,NPR).EQ.0) THEN
VAT( NPL) = VRT
VRN( NPL) = DNU
WRITE(6,2000) LOP, ITN, VRW, VRT, VRP, DNU
WRITE(8,2000) LOP, ITN, VRW, VRT, VRP, DNU
2000 FORMAT(' Loop=',I6,' It=', I3,' VrW=',F7.2,'% ',' VrT=',F7.2,
1 '% ',' VrP=',F7.2,'% ',' Df_Nu =',F8.3 )
END IF
C
C --------------------------------------------------
IF ( DNU.LE.0.01D0) THEN
KKK = 1
WRITE(6,2100) RLN, DT, LOP
WRITE(8,2100) RLN, DT, LOP
2100 FORMAT(/,' *** Stationary condition reached ***',
1 /7X,' ** Ra =',F12.1,' DT =',F11.7,' Loop =',I7 )
END IF
C --------------------------------------------------
IF ( LOP.GE.LPM) THEN
KKK = 1
WRITE(6,2200) RLN, DT, LOP
WRITE(8,2200) RLN, DT, LOP
2200 FORMAT(/,' *** Computation End *** Ra =',F12.1,' DT =',F11.7,
1 ' Loop =',I7,/)
END IF
C --------------------------------------------------
IF ( KKK.EQ.0) GO TO 333
C
WRITE(6,2000) LOP, ITN, VRW, VRT, VRP, DNU
WRITE(8,2000) LOP, ITN, VRW, VRT, VRP, DNU
NPL = NPL + 1
VAW( NPL) = VRW
VVP( NPL) = VRP
VAT( NPL) = VRT
VRN( NPL) = DNU
RETURN 1
C
333 CONTINUE
C --------------------------------------------------
END IF
C ----------------------------------------------------------
IF ( KTF.EQ.0) THEN
C
C -----------------------------------------------------
IF ( MOD( LOP,NPR).EQ.0) THEN
WRITE(6,2300) LOP, ITN, VRW, VRP
WRITE(8,2300) LOP, ITN, VRW, VRP
2300 FORMAT(' Loop=', I6,' It=', I3,' VrW=',F7.2,'% ',' VrP=',F7.2,
1 '% ')
END IF
C -----------------------------------------------------
IF ( VRW.LE.0.1D0 .AND. VRP.LE.0.1D0) THEN
WRITE(6,*) ' '
WRITE(8,*) ' '
WRITE(6,2400) LOP
WRITE(8,2400) LOP
2400 FORMAT(' *** Stationary condition reached *** Loop =', I5 )
WRITE(6,2300) LOP, ITN, VRW, VRP
WRITE(8,2300) LOP, ITN, VRW, VRP
NPL = NPL + 1
VAW( NPL) = VRW
VVP( NPL) = VRP
RETURN 1
END IF
C -----------------------------------------------------
C
END IF
C ----------------------------------------------------------
IF ( LOP.GE.LPM) RETURN 1
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)' *** H F 2 D - FDM ( SV ) *** '
WRITE(6,2100)' '
WRITE(6,2100)' 2D Thermal-Fluid Analysis by the 4th order FDM '
WRITE(6,2100)' '
WRITE(6,2100)' ( Version 3.9 ) '
WRITE(6,2100)' '
WRITE(6,2100)' Copyright 2011 : Yasuhiro MATSUDA '
WRITE(6,2100)' '
WRITE(6,2100)'***************************************************'
2000 FORMAT(/)
2100 FORMAT(A60)
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 /IDT/ DT, LPM, PRN, CTP, RLN, REN, LOP, NPL, NPR
C
C --------------------------------------
DO 100 I = 1, IX-1
C
RAF( I) = HX( I)* HX( I)
100 CONTINUE
C --------------------------------------
DO 200 J = 1, IY-1
C
RAS( J) = HY( J)* HY( J)
200 CONTINUE
C --------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE GINP ( WP1, TP1, PS1, WPT, TPT, PPT, EPT, IX, IY )
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION WP1( IX,IY), TP1( IX,IY), PS1( IX,IY),
1 IDX( 5), IDY( 5), IDE( 5)
C
LOGICAL*1 WPT, TPT, PPT, EPT
DIMENSION WPT( IX,IY,4), TPT( IX,IY,4), PPT( IX,IY,4),
1 EPT( IX,IY,4)
C
NAMELIST /GDAT/ OMG, PSI, TMP, SLP, MOV, CNT, PRS
LOGICAL*1 OMG, PSI, TMP, SLP, MOV, CNT, PRS
C
C ----- SPEC Data READ ---
C
IER = 0
C
C ----------------------------
DO 100 I = 1, IX
DO 100 J = 1, IY
C
WPT( I,J,1) = .TRUE.
PPT( I,J,1) = .TRUE.
TPT( I,J,1) = .TRUE.
C
100 CONTINUE
C ----------------------------
C
C ----- All points are initially in the field ---
C
DO 200 I = 1, IX
C
EPT( I, 1,2) = .TRUE. ! Lower
EPT( I,IY,1) = .TRUE. ! Upper
C
200 CONTINUE
C ----------------------------
C
DO 300 I = 1, IY
C
EPT( 1, I,4) = .TRUE. ! Left
EPT( IX,I,3) = .TRUE. ! Right
C
300 CONTINUE
C ----------------------------
C
C ----- Edge points ---
C
333 CONTINUE
C
OMG = .FALSE.
PSI = .FALSE.
TMP = .FALSE.
SLP = .FALSE.
C
MOV = .FALSE.
CNT = .FALSE.
PRS = .FALSE.
C
IDA = 1
C
READ(10,GDAT,END=999)
C ==========================
C
IF ( OMG ) GO TO 310
IF ( PSI ) GO TO 310
IF ( TMP ) GO TO 310
C
IF ( PRS ) GO TO 310
IF ( CNT ) GO TO 310
IF ( MOV ) GO TO 310
C
310 CONTINUE
C
C ----------------------------------------------------------
DO 400 IJK = 1, IDA
C
ISLP = -10
C
C --------------------------------------
DO 450 L = 1, 5
C
IDX( L) = 0
IDY( L) = 0
IDE( L) = 0
C
450 CONTINUE
C --------------------------------------
WVL = 0.D0
PVL = 0.D0
TVL = 0.D0
RVL = 0.D0
C
READ(10,215,END=999,ERR=225) ISLP, WVL, PVL, TVL,
1 ( IDX( L), IDY( L), IDE( L), L =1,5 ), RVL
C
215 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
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,2000) ISLP, WVL, PVL, TVL, ( IDX( L), IDY( L), IDE( L),
1 L= 1, 5), RVL
2000 FORMAT(5X,I2,1X,3( F6.3,1X),5( 3I3,1X),F6.3)
GO TO 400
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
C IPT = 3, LINE : Upward. IPT = 4, Vertical
C
C ------------------------------------------------
DO 500 L = 1, 5
C
IUPR = IDE( L)
GO TO ( 240, 240, 240, 235 ) IPT
C
235 CONTINUE
ILOW = IDY( L)
GO TO 245
C
240 ILOW = IDX( L)
245 I = IDX( L)
J = IDY( L)
C
C --------------------------------------
DO 550 ILU = ILOW, IUPR
C
IF ( I) 255, 500, 250
250 IF ( J) 255, 500, 251
251 IF ( I-IX) 252, 252, 255
252 IF ( J-IY) 260, 260, 255
C
C ---- IF I or J = 0 : No Boundary specified ---
C
255 IER = 1
C
WRITE (6,2100) L, IX, IY
2100 FORMAT (' * Error in SPEC No.', I1, 2X,' IX =', I3,' IY =', I3)
WRITE ( 16,2000) ISLP, WVL, PVL, TVL,
1 ( IDX( M), IDY( M), IDE( M), M = 1, 5 )
GO TO 500
C
260 IF (.NOT.OMG) GO TO 265
WPT( I,J,1) = .FALSE.
WPT( I,J,2) = SLP
WPT( I,J,3) = MOV
C
WP1 ( I,J) = WVL
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) = PVL
C
270 IF (.NOT.TMP) GO TO 275
TPT( I,J,1) = .FALSE.
TPT( I,J,2) = SLP
TPT( I,J,3) = MOV
TP1( I,J) = TVL
C
275 IF (.NOT.CNT) GO TO 278
PPT( I,J,4) = .TRUE.
TPT( I,J,4) = .TRUE.
C
278 CONTINUE
IF (.NOT.PRS) GO TO 279
C
279 CONTINUE
IF (.NOT.MOV) GO TO 280
WPT( I,J,3) = .TRUE.
C
280 GO TO ( 286, 288, 287, 285 ), IPT
C
285 J = J + 1
GO TO 550
C
286 J = J - 1
GO TO 288
C
287 J = J + 1
288 I = I + 1
C
550 CONTINUE
C --------------------------------------
C
500 CONTINUE
C ------------------------------------------------
C
400 CONTINUE
C ----------------------------------------------------------
GO TO 333
C
999 RETURN
END
C **********************************************************************
C
SUBROUTINE IGSS ( PS1, PPT, TP1, TPT, XZ, YZ, IX, IY )
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION PS1( IX,IY), TP1( IX,IY), XZ ( IX), YZ ( IY)
C
LOGICAL*1 TPT( IX,IY,4), PPT( IX,IY,4)
C
COMMON /CTR/ KTF, MSH, IHS, WEI
COMMON /IDT/ DT, LPM, PRN, CTP, RLN, REN, LOP, NPL, NPR
C
N = IY-1
C
C ---------------------------------------
DO 100 I = 1, IX
DO 100 J = 1, IY
C
IF (.NOT.PPT( I,J,1)) GO TO 100
PS1( I,J) = PS1( IX,J)
C
100 CONTINUE
C ---------------------------------------
IF ( KTF.NE.1) RETURN ! Fluid Flow Analysis
IF ( RLN.LE.0.00001D0) RETURN
C
C ----- Natural Convection Problem ---
C
DTMP = 0.D0 ! Referenced Temperature
C ============================================
C
C --------------------------------------------
DO 200 I = 1, IX
DO 200 J = 1, IY
C
IF ( TPT( I,J,1) ) GO TO 220
IF (.NOT.TPT( I,J,4)) GO TO 240
220 TP1( I,J) = DTMP
240 CONTINUE
C
200 CONTINUE
C --------------------------------------------
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 333
C
UMX = -1000.D0
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
333 CONTINUE
IF ( IH.NE.2) GO TO 550
UMX = 0.D0
I = IU
C
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
550 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 ( TE, XX, YY, NDE, PE, PS, T, WP,
1 UB, VB, PS1, U, V, XZ, YZ, WP1,
2 TP1, SKN, GG, SIM, ANU, HX, HY, JU,
3 IV, IX, IY, IXY, IXY1, IX1, IY1 )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
COMMON /ABC/ DNU, RN0, RN1
COMMON /CTR/ KTF, MSH, IHS, WEI
COMMON /IDT/ DT, LPM, PRN, CTP, RLN, REN, LOP, NPL, NPR
C
COMMON /DEF/ VAW( 250000), VAT( 250000), VVP( 250000),
1 VRN( 250000)
C
DIMENSION TE( IXY), XX( IXY), YY ( IXY), NDE( IXY1,4),
1 PE( IXY), PS( IX,IY), T ( IX,IY), WP ( IX,IY),
2 UB( IX,IY), VB( IX,IY), PS1( IX,IY), U ( IX,IY),
3 V ( IX,IY)
C
DIMENSION XZ ( IX), YZ( IY), WP1( IX,IY), TP1( IX,IY),
1 SKN( IX,IY), GG( IX,IY), SIM( IX), ANU( IX),
2 HX ( IX1), HY( IY1)
C
NNP = IX* IY
NEL = (IX-1)* (IY-1)
C
C -----------------------------------------------------------
C
CALL STND ( XX, YY, NDE, NNP, NEL, XZ, YZ, IX, IY )
C
C -----------------------------------------------------------
C
CALL CKUV ( U, UB, DT, DUX, IPU, JPU, IX, IY )
C
C -----------------------------------------------------------
C
CALL CKUV ( V, VB, DT, DVX, IPV, JPV, IX, IY )
C
C -----------------------------------------------------------
C
CALL STMP ( IX, IY, PS1, PE, NNP )
C
C -----------------------------------------------------------
WRITE(7,2000) KTF, RLN, REN, DT, LOP, MSH, NPL
2000 FORMAT(I2,F12.1, F8.2, F11.8, I7, I2, I6)
C
WRITE(7,2100) ( PE ( I), I= 1, IXY)
WRITE(7,2100) ( VAW( I), I= 2, NPL)
WRITE(7,2100) ( VVP( I), I= 2, 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(12F10.3)
WRITE(7,2200) (( NDE( I,J), J = 1, 4), I = 1, IXY1)
2200 FORMAT(20I6)
C
C ----------------------------
IF ( KTF.EQ.1) THEN
C -------------------------------------------
C
CALL STMP ( IX, IY, TP1, TE, NNP )
C
C -------------------------------------------
WRITE(7,2100) ( TE ( I), I = 1, IXY)
WRITE(7,2100) ( VAT( I), I = 2, NPL)
WRITE(7,2100) ( VRN( I), I = 2, NPL)
END IF
C ------------------------------------------------
C
CALL DTOP ( PS1, U, V, XZ, YZ, IX, IY )
C
C ------------------------------------------------
IF ( KTF.EQ.1) THEN
C --------------------------------------------------------------------
C
CALL CMNU ( U, TP1, IX, IY, IX1, IY1, SKN, GG, SIM,
1 ANU, HX, HY, XZ, YZ )
C
C --------------------------------------------------------------------
WRITE(6,2300) DNU
WRITE(8,2300) DNU
2300 FORMAT(' * Abs. Diff. of Nu =', F7.4,' *** End ***')
END IF
C ------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE DTOP ( PS1, U, V, XZ, YZ, IX, IY )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
COMMON /CTR/ KTF, MSH, IHS, WEI
COMMON /IDT/ DT, LPM, PRN, CTP, RLN, REN, LOP, NPL, NPR
C
DIMENSION PS1( IX,IY), U( IX,IY), V( IX,IY), XZ( IX), YZ( IY)
C
C ----- Attention : Approx. for equal length of meshes ---
C
DY = 1.D0/ DFLOAT( IY-1)
C =============================
C
IF ( KTF.EQ.1) RT = CTP* DT/ DY** 2
C
RW = PRN* DT / DY** 2
BV = RLU* DT / DY
C ==========================
C
IF ( RLV.GT.RLU) BV = RLV* DT/DY
C
IF ( KTF.EQ.0) WRITE(6,2000) RW, BV
IF ( KTF.EQ.0) WRITE(8,2000) RW, BV
2000 FORMAT(' * r2, b =',2F9.4)
C
IF ( KTF.EQ.1) WRITE(6,2100) RT, RW, BV
IF ( KTF.EQ.1) WRITE(8,2100) RT, RW, BV
2100 FORMAT(' * r1, r2, b =',3F9.4)
C
C --------------------------------------
DO 100 I = 1, IX
DO 100 J = 1, IY
C
PS1( I,J) = DABS( PS1( I,J))
100 CONTINUE
C --------------------------------------
IM = ( IX + 1)/ 2
JM = ( IY + 1)/ 2
PMID = PS1( IM, JM)
C -----------------------------------------------------
C
CALL MAXS ( PS1, IX, IY, PMX, IP, JP, 1 )
C
C -----------------------------------------------------
XP = XZ( IP)
YP = YZ( JP)
WRITE(6,2200) PMID, PMX, XP, YP
WRITE(8,2200) PMID, PMX, XP, YP
2200 FORMAT(/,' * P_mid.=',F8.4,' * P_max.=',F8.4,' at Xp =',F8.4,
1 ' Yp =',F8.4 )
C -------------------------------------------------
C
CALL MAXS ( U, IX, IY, RLU, IAU, JAU, 1 )
C
C -------------------------------------------------
C
CALL MAXS ( V, IX, IY, RLV, IAV, JAV, 1 )
C
C -------------------------------------------------
XAU = XZ( IAU)
YAU = YZ( JAU)
XAV = XZ( IAV)
YAV = YZ( JAV)
C -------------------------------------------------
C
CALL MAXS ( U, IX, IY, UMX, IM, JU, 2 )
C
C -------------------------------------------------
C
CALL MAXS ( V, IX, IY, VMX, IV, JM, 3 )
C
C -------------------------------------------------
XU = XZ( IM)
YU = YZ( JU)
XV = XZ( IV)
YV = YZ( JM)
WRITE(6,2300) UMX, YU, VMX, XV
WRITE(8,2300) UMX, YU, VMX, XV
2300 FORMAT(/,' * Umax.( X = 0.5)=',F10.4,' ( Y = ',F7.4,')',/,
1 ' * Vmax.( Y = 0.5)=',F10.4,' ( X = ',F7.4,')' )
C
RETURN
END
C **********************************************************************
C
SUBROUTINE PSCL ( PS1, PS, WP1, TP1, PPT, EPT, ACP, HX,
1 HY, RDP, IX, IY, IX1, IY1 )
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION PS1( IX,IY), PS( IX,IY), WP1( IX,IY), TP1( IX,IY),
1 RDP( IX,IY), HX( IX1), HY ( IY1)
C
LOGICAL*1 PPT( IX,IY,4), EPT( IX,IY,4), DRT
C
COMMON /UDT/ DY2, ITN
COMMON /IDT/ DT, LPM, PRN, CTP, RLN, REN, LOP, NPL, NPR
COMMON /CM1/ PA1( 101,101), PA2( 101,101), PA3( 101,101),
1 PB1( 101,101), PB2( 101,101), PB3( 101,101)
C
FS = 1.D0 ! FS = ( DX/ DY)** 2
CFI = 1.D0
DXS = DY2* FS ! DXS = DX* DX DY2 = DY* DY
C
UNU = 0.D0
RND = 0.D0
C
ITN = 0
EPS = 0.00005D0
DRT = .FALSE.
JDY = IY + 1
C
105 CONTINUE
IDC = 0
ITN = ITN + 1
IF ( ITN.GT.500) GO TO 777
C
C ----------------------------------------------------------
DO 200 JJ = 1, IY
C
10 CONTINUE
DRT = .NOT. DRT
IF ( DRT ) J = JJ
IF (.NOT. DRT ) J = JDY - JJ
C
C ------------------------------------------------
DO 100 I = 1, IX
C
DXH = HX( I-1)* HX( I)* ( HX( I-1) + HX( I))
DYH = HY( J-1)* HY( J)* ( HY( J-1) + HY( J))
DXY = ( HY( J-1) + HY( J))* DXH + ( HX( I-1) + HX( I))* DYH
C
IF ( PPT( I,J,1)) GO TO 1
IF ( PPT( I,J,2)) GO TO 100
GO TO 100
C
1 CONTINUE
IF ( EPT( I,J,3)) GO TO 2
IF ( EPT( I,J,4)) GO TO 12
IF ( EPT( I,J,1)) GO TO 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))
PS1( I,J) = PS1( I,J) + ACP* ( RDP( I,J) - PS1( I,J))
C
IF ( RAB2.GT.EPS) GO TO 5
IF ( ITN.LT.10) GO TO 5
GO TO 100
C
5 IDC = 1
IXN = I
IYN = J
GO TO 100
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 100
C
C ----- Continuative right ---
C
8 CONTINUE
PS1( I,J)
1 = ( HY( J)* PS1( I,J-1) + HY( J-1)* PS1( I,J+1)
2 + 0.5D0* DYH* WP1( I,J))/( HY( J-1) + HY( J))
GO TO 100
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) + DY2* WP1( I,J)) )
2 * CFI
C
RDL = (( HX( I-1)* PS1( I+1,J)+HX( I)* PS1( I-1+( IX-1),J))* DYH
1 + ( HY( J-1)* PS1( I,J+1) + HY( J)* PS1( I,J-1))* DXH
2 + 0.5D0* DXH* DYH* WP1( I,J))/ DXY
C
PS1( I,J) = PS1( I,J) + ACP* ( RDL - PS1( I,J))
GO TO 100
C
C ----- Continuative left ---
C
18 PS1( I,J)
1 = (( HX( I-1)* PS1( I,J) + HX( I)* PS1( I,J) )* DYH
2 + ( HY( J-1)* PS1( I,J+1) + HY( J)* PS1( I,J-1))* DXH
3 + 0.5D0* DXH* DYH* WP1( I,J))/ DXY
GO TO 100
C
22 IF ( PPT( I,J,4)) GO TO 28
PS1( I,J) = PS1( I,J-(IY-1))
GO TO 100
C
28 PS1( I,J) = ( HX( I)* PS1( I-1,J) + HX( I-1)* PS1( I+1,J)
1 + 0.5D0* DXH* WP1( I,J))/( HX( I-1) + HX( I))
GO TO 100
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) + DXS* WP1( I,J))
2 / FS )* CFI/ FS
C
PS1( I,J) = PS1( I,J) + ACP* ( DELP - PS1( I,J))
GO TO 100
C
38 PS1( I,J) = ( HX( I)* PS1( I-1,J) + HX( I-1)* PS1( I+1,J)
1 + 0.5D0* DXH* WP1( I,J))/( HX( I-1) + HX( I))
GO TO 100
C
100 CONTINUE
C ------------------------------------------------
IF ( DRT) GO TO 10
C
200 CONTINUE
C ----------------------------------------------------------
IF ( DABS( RND).LE.0.000001D0) GO TO 101
PSS = 0.D0
C
C ----------------------------------------------------------
DO 300 I = 1, IX
C
PSS = PSS + PS1( I,2) + 0.5D0* DY2* WP1( I,1)
300 CONTINUE
C ----------------------------------------------------------
C
UL = 0.D0
C ==============
PSS = ( PSS/ IX)
C ----------------------------------------------------------
DO 400 I = 1, IX
C
PS1( I,1) = PSS
400 CONTINUE
C ----------------------------------------------------------
C
101 IF ( DABS( UNU).LE.0.000001D0) GO TO 102
PSS = 0.D0
C
C ----------------------------------------------------------
DO 500 I = 1, IX
C
PSS = PSS + PS1( I,IY-1) + 0.5D0* DY2* WP1( I,IY)
500 CONTINUE
C ----------------------------------------------------------
C
UU = 0.D0
C ==============
C
C ----------------------------
DO 600 I = 1, IX
C
PS1( I,IY) = PSS
600 CONTINUE
C ----------------------------
102 CONTINUE
IF ( IDC.GT.0) GO TO 105
RETURN
C
777 CONTINUE
WRITE(6,2000) LOP
2000 FORMAT (/' *** Diverged ! ITN > 500 in PSCL * Incompl. sol.of PS
1I:* STOP at Loop =',I4)
STOP
C
END
C **********************************************************************
C
SUBROUTINE PRPS ( PPT, WPT, TPT, EPT, PS1, WP1, TP1, T,
1 PS, WP, U, V, UB, VB, SKN, GG,
2 HH, BX, BY, XZ, YZ, SIM, ANU, RAF,
3 RAS, HX, HY, IX, IY, IX1, IY1 )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
LOGICAL*1 PPT( IX,IY,4), WPT( IX,IY,4), TPT( IX,IY,4),
1 EPT( IX,IY,4)
C
DIMENSION PS1( IX,IY), WP1( IX,IY), TP1( IX,IY), T ( IX,IY),
1 PS ( IX,IY), WP ( IX,IY), U ( IX,IY), V ( IX,IY),
2 UB ( IX,IY), VB ( IX,IY), SKN( IX,IY), GG( IX,IY),
3 HH ( IX,IY), BX ( IX,IY), BY ( IX,IY), XZ ( IX),
4 YZ ( IY), SIM( IX), ANU( IX), RAF( IX1),
5 RAS( IY1), HX ( IX1), HY ( IY1)
C
COMMON /CTR/ KTF, MSH, IHS, WEI
C
C ----------------------------
C
CALL TTLE
C
C ----------------------------
WRITE(6,*) ' '
WRITE(6,*) ' ** KTF = 1: Thermal & Fluid Analysis = 0 : Fluid An
1alysis ** KTF = ?'
READ(5,*) KTF
** KTF = 1
C
C -------------------------
C
CALL CDAT ( IY)
C
C -------------------------
C
C ----- Initialization -----------------------------------------------
C
CALL SINT ( PPT, WPT, TPT, EPT, PS1, WP1, TP1, T, PS,
1 WP, U, V, UB, VB, SKN, GG, HH, BX,
2 BY, XZ, YZ, SIM, ANU, IX, IY )
C
C ----- Mesh generation ----------------------------------------------
C
CALL GMSH ( XZ, YZ, HX, HY, IX, IY, IX1, IY1 )
C
C ----------------------------------------------------------
C
CALL DINP ( HX, HY, RAF, RAS, IX, IY, IX1, IY1 )
C
C ----- Geometry Data --------------------------------------
C
CALL GINP ( WP1, TP1, PS1, WPT, TPT, PPT, EPT, IX, IY )
C
C ----------------------------------------------------------
C
CALL IGSS ( PS1, PPT, TP1, TPT, XZ, YZ, IX, IY )
C
C ----------------------------------------------------------
DO 100 I = 1, IX
DO 100 J = 1, IY
C
PS( I,J) = PS1( I,J)
WP( I,J) = WP1( I,J)
T ( I,J) = TP1( I,J)
C
100 CONTINUE
C --------------------------------------------------
C
CALL PSCF ( IX, IY, HX, HY, IX1, IY1 )
C
C --------------------------------------------------
C
CALL TCOF ( HX, HY, IX, IY, IX1, IY1 )
C
C --------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE CDAT ( IY)
C
IMPLICIT REAL*8 (A-H,O-Z)
C
COMMON /UDT/ DY2, ITN
COMMON /CTR/ KTF, MSH, IHS, WEI
COMMON /IDT/ DT, LPM, PRN, CTP, RLN, REN, LOP, NPL, NPR
C
C ----------------------------------------------------------
C
C KTF = 1 : Heat & Fluid Flow ( Natural Convection)
C = 0 : Fluid Flow ( Cavity Flow )
C
C MSH = 1 : Mesh of Equal length ( EM )
C = 2 : Mesh of Unequal length ( UEM )
C
C IHS = 1 : With Correction
C
C WEI : Weighting parameter
C
C ----------------------------------------------------------
C
WRITE(6,*) ' '
WRITE(6,*) ' ** MSH = 1 : Mesh of Equal length 2 : Mesh of Unequ
1al length'
C
READ (5,*) MSH
** MSH = 2
C
IHS = 1
*** IHS = 0
C
C ----------------
C
WEI = 0.D0
WEI = 1.D0
WEI = 0.5D0 ! Better than = 0 for Equal length meshes ---
C ================================================================
C
WRITE(6,2000) KTF, MSH, IHS, WEI
WRITE(8,2000) KTF, MSH, IHS, WEI
2000 FORMAT(/,' * KTF =', I2,' MSH =', I2,' IHS =', I2,' WEI =',F5.2)
C
WRITE(6,*) ' ** NPR = ?'
** READ(5,*) NPR
C
NPR = 100
NPR = 200
NPR = 1000
IF ( KTF.EQ.0) NPR = 200
IF ( MSH.EQ.1 .AND. KTF.EQ.0 ) NPR = 10 ! Equal length
C
WRITE(6,2050) NPR
2050 FORMAT(' * NPR =',I5)
C
C ----------------------------------------------------------
IF ( KTF.EQ.1) THEN
C
PRN = 0.71D0 ! Prandtle noumber
C
CTP = 1.D0 ! = 1/( Re* Pr )
C ----------- Attention --------------------------
C
WRITE(6,*) ' '
WRITE(6,2100) PRN, CTP
2100 FORMAT(' *** Pr =', F6.3, ' CTP =',F5.2)
C
WRITE(6,*) ' '
WRITE(6,*) ' ** LOP = ?'
** READ(5,*) LPM
C
LPM = 300
LPM = 1000
LPM = 35000
C
LPM = 10000 ! 20x20, UE, Ra=100000'
C
WRITE(6,*) ' ** Ra = ?'
** READ(5,*) RLN
RLN = 1000000.D0
RLN = 100000.D0
C
C --------------------
* RLN = 1.D3
* DT = 1.D-4
C --------------------
* RLN = 1.D4
* DT = 1.D-4
C --------------------
* RLN = 1.D5
* DT = 1.D-4
C --------------------
* RLN = 1.D6
* DT = 1.D-4
C --------------------
* RLN = 1.D7
* DT = 5.D-6 ! W=0
C --------------------
* RLN = 1.D8
* DT = 3.D-7 ! W=0
C --------------------
C
WRITE(6,*) ' ** DELT= ?'
*** READ(5,*) DT
C
DT = 0.000 01 D0 ! 50x50, UE, Ra=100000
C
DT = 0.000 05 D0 ! 20x20, UE, Ra=100000
C ----------------^----
C
C ---- Attention ---
C
WRITE(6,2200) RLN, DT, LPM
WRITE(8,2200) RLN, DT, LPM
2200 FORMAT(/,' * Ra =', F12.1,' DT =',F11.7,' Final Loop =',I7 )
C
END IF
C ----------------------------------------------------------
IF ( KTF.EQ.0) THEN
C
WRITE(6,*) ' ** Re = ?'
READ (5, *) REN
** REN = 100.D0
WRITE(6,*) ' '
WRITE(6,2300) REN
2300 FORMAT(' *** Re =', F8.2)
C
PRN = 1.D0/ REN
C ====================
C
WRITE(6,*) ' '
WRITE(6,*) ' ** DT = ? & Loop = ?'
** READ(5,*) DT, LPM
DT = 0.001D0
DT = 0.005D0
DT = 0.04D0
C
IF ( MSH.NE.2) GO TO 3355 ! In case of Unqual length
IF ( DABS( REN - 100.D0).LE.0.01D0 ) DT = 0.009D0
IF ( DABS( REN - 1000.D0).LE.0.01D0 ) DT = 0.008D0
3355 CONTINUE
C
IF ( MSH.NE.1) GO TO 3357 ! In case of Equal length
IF ( DABS( REN - 100.D0).LE.0.01D0 ) DT = 0.04D0
IF ( DABS( REN - 1000.D0).LE.0.01D0 ) DT = 0.015D0
IF ( DABS( REN - 1000.D0).LE.0.01D0 ) NPR = 100
3357 CONTINUE
C
LPM = 5000
C
RLN = 0.D0
C
WRITE(6,2400) DT, LPM
WRITE(8,2400) DT, LPM
2400 FORMAT(/,' *** Cavity Flow Problem * DT =',F11.7,' Loop =',I7 )
C
END IF
C ---------------------------------------------------------
C
C ----- Attention : Approx. for equal length of meshes --- A : DY ---
C
DY = 1.D0/ DFLOAT( IY-1)
DY2 = DY* DY
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 /CTR/ KTF, MSH, IHS, WEI
C
C --------------------------------------------------
C
C MSH = 1 : Meshes of equal length ----- EM
C
C MSH = 2 : JSME - Paper Version ------- UEM
C
C --------------------------------------------------
WRITE(6,*) ' '
C
IF ( MSH.EQ.1) GO TO 150
IF ( MSH.EQ.2) GO TO 230
C
150 CONTINUE
C --------------------------------------------------------------------
WRITE(6,*) '* Meshes of equal length '
WRITE(8,*) '* Meshes of equal length '
C
C --------------------------------------
DO 100 I = 1, IX
C
XZ( I) = DFLOAT( I-1)/ DFLOAT( IX-1)
100 CONTINUE
C --------------------------------------
DO 200 J = 1, IY
C
YZ( J) = DFLOAT( J-1)/ DFLOAT( IY-1)
200 CONTINUE
C --------------------------------------
GO TO 34
C
230 CONTINUE
C -------------------------------------------------------------
WRITE(6,*) '*** Meshes of Unequal length ( UEM ) ***'
WRITE(8,*) '*** Meshes of Unequal length ( UEM ) ***'
WRITE(6,*) ' '
C
AA = 2.D0
C
DCX = 2.D0/ DFLOAT( IX-1)
DCY = 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
C --------------------------------------
DO 300 I = 1, MX1
C
XZ( I) = ( DEXP( AA* DFLOAT( I-1)* DCX) - 1.D0)
1 /( 2.D0* ( DEXP( AA) - 1.D0))
C
XZ( IX+1-I ) = 1.D0 - XZ( I)
C
300 CONTINUE
C --------------------------------------
DO 400 J = 1, MY1
C
YZ( J) = ( DEXP( AA* DFLOAT( J-1)* DCY) - 1.D0)
1 /( 2.D0* ( DEXP( AA) - 1.D0))
C
YZ( IY+1-J) = 1.D0 - YZ( J)
C
400 CONTINUE
C --------------------------------------
34 CONTINUE
C
C --------------------------------------
DO 500 I = 1, IX-1
C
HX( I) = DABS( XZ( I+1) - XZ( I))
500 CONTINUE
C --------------------------------------
DO 600 J = 1, IY-1
C
HY( J) = DABS( YZ( J+1) - YZ( J))
600 CONTINUE
C --------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE SINT ( PPT, WPT, TPT, EPT, PS1, WP1, TP1, T,
1 PS, WP, U, V, UB, VB, SKN, GG,
2 HH, BX, BY, XZ, YZ, SIM, ANU, IX,
3 IY )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
LOGICAL*1 PPT( IX,IY,4), WPT( IX,IY,4), TPT( IX,IY,4),
1 EPT( IX,IY,4)
C
DIMENSION PS1( IX,IY), WP1( IX,IY), TP1( IX,IY),
1 T (IX, IY), PS ( IX,IY), WP ( IX,IY),
2 U ( IX,IY), V ( IX,IY), UB ( IX,IY),
3 VB ( IX,IY), SKN( IX,IY), GG ( IX,IY),
4 HH ( IX,IY), BX ( IX,IY), BY ( IX,IY),
5 XZ ( IX), YZ ( IY), SIM( IX),
6 ANU( IX)
C
C ----------------------------------------------------------
DO 100 I = 1, IX
C
C ------------------------------------------------
DO 130 J = 1, IY
C
C --------------------------------------
DO 150 K = 1, 4
C
PPT ( I,J,K) = .FALSE.
WPT ( I,J,K) = .FALSE.
TPT ( I,J,K) = .FALSE.
EPT ( I,J,K) = .FALSE.
C
150 CONTINUE
C --------------------------------------
PS1( I,J) = 0.D0
WP1( I,J) = 0.D0
TP1( I,J) = 0.D0
T ( I,J) = 0.D0
C
PS ( I,J) = 0.D0
WP ( 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
SKN( I,J) = 0.D0
C
GG ( 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
130 CONTINUE
C ------------------------------------------------
SIM( I) = 0.D0
ANU( I) = 0.D0
C
100 CONTINUE
C ----------------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE VARC ( WP1, WB, VRW, IX, IY )
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION WP1( IX,IY), WB( IX,IY)
C
COMMON /IDT/ DT, LPM, PRN, CTP, RLN, REN, LOP, NPL, NPR
C
C --------------------------------------
WL1 = 0.D0
DO 100 I = 1, IX
DO 100 J = 1, IY
C
ABS1 = DABS( WP1( I,J))
IF ( ABS1.GT.WL1) WL1 = ABS1
C
100 CONTINUE
WL2 = WL1* 0.01D0
C --------------------------------------
WDL = 0.D0
DO 200 I = 1, IX
DO 200 J = 1, IY
C
ABSW = DABS( WP1( I,J))
IF ( ABSW.LE.WL2) GO TO 200
WDS = DABS((WP1( I,J) - WB( I,J))/WP1( I,J))
IF ( WDS.GT.WDL) THEN
WDL = WDS
II = I
JJ = J
END IF
C
200 CONTINUE
C --------------------------------------
VRW = WDL* 100.D0
DO 300 I = 1, IX
DO 300 J = 1, IY
C
WB( I,J) = WP1( I,J)
300 CONTINUE
C --------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE CKUV ( U, UB, DT, DUX, IP, JP, IX, IY )
C
C ----- Unsteady term ---
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION U( IX,IY), UB( IX,IY)
C
C ------------------------------------------------
DUX = 0.D0
DO 100 I = 1, IX
DO 100 J = 1, IY
C
DIFU = DABS(( U( I,J) - UB( I,J))/ DT )
IF ( DIFU.LT.DUX) GO TO 100
DUX = DIFU
IP = I
JP = J
C
100 CONTINUE
C ------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE TMCL ( TP1, T, WP1, TPT, EPT, U, V,
1 HH, HX, HY, BX, BY, IX, IY,
2 IX1, IY1, LTC )
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION TP1( IX,IY), T ( IX,IY), WP1( IX,IY), U ( IX,IY),
1 V ( IX,IY), HH( IX,IY), BX ( IX,IY), BY( IX,IY),
2 HX ( IX1), HY( IY1)
C
COMMON /CM2/ TM1( 101), TM2( 101), TM3( 101), TY1( 101),
1 TY2( 101), TY3( 101)
C
LOGICAL*1 TPT( IX,IY,4), EPT( IX,IY,4)
C
COMMON /UDT/ DY2, ITN
COMMON /CTR/ KTF, MSH, IHS, WEI
COMMON /IDT/ DT, LPM, PRN, CTP, RLN, REN, LOP, NPL, NPR
C
FS = 1.D0
CF = 1.D0
C ==============
C
C ----- Attention : Approx. for Equal length of meshes ---
C
DO 100 I = 1, IX
DO 100 J = 1, IY
C
BX( I,J) = U( I,J)* DT* IX1
BY( I,J) = V( I,J)* DT* IY1
C
100 CONTINUE
C -------------------------------------------------
GO TO ( 301, 304), LTC
C
301 RH = CTP ! Energy Equation ==> 1 /( Re*Pr )
C =============
C
GO TO 303
C
304 RH = PRN ! For Vorticity Equation ==> Prandtle Number
C =============
C
303 CONTINUE
C ----------------------------------
IF ( IHS.NE.1) GO TO 81
C
C ----- Correction Method ---
C
C ----- Attention : Approx. For equal length meshes ---
C
RX = RH* DT* ( IX1** 2)
RY = RH* DT* ( IY1** 2)
C
C ---------------------------------
DO 200 I = 1, IX
DO 200 J = 1, IY
C
HH( I,J) = ( 1.D0 + ( BX( I,J)* BY( I,J))/( RX + RY))* RH* DT
200 CONTINUE
C ---------------------------------
GO TO 85
C
81 CONTINUE
C ---------------------------------
DO 300 I = 1, IX
DO 300 J = 1, IY
C
HH( I,J) = RH* DT
300 CONTINUE
C ---------------------------------
C
85 CONTINUE
C ---------------------------------
DO 400 I = 1, IX
DO 400 J = 1, IY
C
IF ( TPT( I,J,1)) GO TO 400 ! Fluid
TP1( I,J) = T( I,J) ! Boundary
C
400 CONTINUE
C ------------------------------------------------
DO 700 J = 1, IY
C
C ----------------------------------
DO 750 I = 1, IX
C
IF ( TPT( I,J,1)) GO TO 1
IF ( TPT( I,J,2)) GO TO 101
IF ( TPT( I,J,3)) GO TO 555
C
C ----- Not Fluid, Not conducting, Not insulated ---
GO TO 101
C
1 CONTINUE
IF ( EPT( I,J,3)) GO TO 2 ! Right
IF ( EPT( I,J,4)) GO TO 12 ! Left
IF ( EPT( I,J,1)) GO TO 22 ! Upper
IF ( EPT( I,J,2)) GO TO 32 ! Lower
C
C ----- Normal interior to fluid point ---
C
TP1( I,J) = TP1( I,J) + 2.D0* HH( I,J)
C
1 * ( TM1( I)* T( I-1,J) - TM2( I)* T( I,J)
2 + TM3( I)* T( I+1,J) + TY1( J)* T( I,J-1)
3 - TY2( J)* T( I,J) + TY3( J)* T( I,J+1))
GO TO 333
C
2 CONTINUE
IF ( TPT( 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
TP1( I,J) = TP1( I,J) + 2.D0* HH( I,J)
1 * ( A1* T( I-2,J) - B1* T( I-1,J) + C1* T( I,J)
2 + A2* T( I,J-1) - B2* T( I,J) + C2* T( I,J+1))
GO TO 333
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
TP1( I,J) = TP1( I,J)
1 + 2.D0* HH( I,J)* ( A2* T( I,J-1) - B2* T( I,J) + C2* T( I,J+1))
C
GO TO 333
C
12 CONTINUE
IF ( TPT( 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
TP1( I,J) = TP1( I,J)
1 + 2.D0* HH( I,J)* ( A1* T( I,J) - B1* T( I+1,J) + C1* T( I+2,J)
2 + A2* T( I,J-1) - B2* T( I,J)
3 + C2* T( I,J+1))
GO TO 333
C
C ----- Continuative left ---
C
18 CONTINUE
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
TP1( I,J) = TP1( I,J)
1 + 2.D0* HH( I,J)* ( A2* T( I,J-1) - B2* T( I,J) + C2* T( I,J+1))
C
GO TO 333
C
22 CONTINUE
C
IF ( TPT( I,J,4)) GO TO 28
C
C ----- Periodic Upper ---
C
GO TO 333
C
C ----- Continuative Upper ---
C
28 CONTINUE
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))*( 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))* ( HX( I-2) + HX( I-1) + HX( I))
1 *( HX( I-2) + HX( I-1) + HX( I) + HX( I+1))
2 + 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)) - 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
TP1( I,J) = TP1( I,J)
C
1 + HH( I,J)* ( A1* T( I-2,J) + B1* T( I-1,J) + C1* T( I,J)
2 + D1* T( I+1,J) + E1* T( I+ 2,J) + A2* T( I,J-2)
3 - B2* T( I,J-1) + C2* T( I,J))
GO TO 333
C
29 CONTINUE
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
TP1( I,J) = TP1( I,J) + 2.D0* HH( I,J)
1 * ( A1* T( I-1,J) - B1* T( I,J) + C1* T( I+1,J))
GO TO 333
C
32 CONTINUE
IF ( TPT( 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)) )
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
TP1( I,J) = TP1( I,J) + 2.D0* HH( I,J)
1 * ( A1* T( I-1,J) - B1* T( I,J) + C1* T( I+1,J)
2 + A2* T( I,J) - B2* T( I,J+1) + C2* T( I,J+2))
GO TO 333
C
C ----- Continuative Under ---
C
38 CONTINUE
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))
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
TP1( I,J) = TP1( I,J)
1 + HH( I,J)* ( A1* T( I-2,J) + B1* T( I-1,J) + C1* T( I,J)
2 + D1* T( I+1,J) + E1* T( I+ 2,J) + A2* T( I,J)
3 - B2* T( I,J+1) + C2* T( I,J+2))
GO TO 333
C
39 CONTINUE
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
TP1( I,J) = TP1( I,J) + 2.D0* HH( I,J)
1 * ( A1* T( I-1,J) - B1* T( I,J) + C1* T( I+1,J))
GO TO 333
C
C ----- Insulated wall ---
C
555 CONTINUE
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 310
IF ( EPT( I,J,2)) GO TO 206
C
403 TP1( I,J) = TP1( I,J) + HH( I,J)* ( 2.D0* T( I,J+1)
1 + ( 2.D0* T( I+1,J) - CF* T( I,J))/ FS )
GO TO 333
C
404 TP1( I,J) = TP1( I,J) + HH( I,J)* ( 2.D0* T( I,J+1)
1 + ( 2.D0* T( I-1,J) - CF* T( I,J))/ FS )
GO TO 333
C
405 TP1( I,J) = TP1( I,J) + HH( I,J)* ( 2.D0* T( I,J-1)
1 + ( 2.D0* T( I+1,J) - CF* T( I,J))/ FS )
GO TO 333

406 TP1( I,J) = TP1( I,J) + HH( I,J)* ( 2.D0*T( I,J-1)
1 + ( 2.D0* T( I-1,J) - CF* T( I,J))/ FS )
GO TO 333
C
203 TP1( I,J) = TP1( I,J) + HH( I,J)* ( T( I,J+1) + T( I,J-1)
1 + ( T( I+1,J) + T( I-1,J) - CF* T( I,J))/ FS )
GO TO 333
C
204 TP1( I,J) = TP1( I,J) + HH( I,J)* ( T( I,J+1) + T( I,J-1)
1 + ( T( I-1,J) + T( I+1,J) - CF* T( I,J))/ FS )
GO TO 333
C
310 TP1( I,J) = TP1( I,J) + HH( I,J)* ( T( I,J-1) + T( I,J-1)
1 + ( T( I+1,J) + T( I-1,J) - CF* T( I,J))/ FS )
GO TO 333
C
206 TP1( I,J) = TP1( I,J) + HH( I,J)* ( T( I,J+1) + T( I,J+1)
1 + ( T( I+1,J) + T( I-1,J) - CF* T( I,J))/ FS )
GO TO 333
C
C ----- Perfectly conducting wall ---
C
101 GO TO 333
C
333 CONTINUE
C ------------------------------------
C
750 CONTINUE
C --------------------------------------
C
700 CONTINUE
C ------------------------------------------------
C
IF ( LTC.NE.2) RETURN ! Energy Equation
C
IF ( RLN) 900, 900, 40
C
C ----- Natural Convection Problem ---
C
40 CONTINUE
C
C ----- Attention : Bouyancy Term comp. ---
C
C --------------------------------------
DO 800 J = 1, IY
DO 800 I = 1, IX
C
IF ( EPT( I,J,3)) GO TO 52 ! Right Boundary
IF ( EPT( I,J,4)) GO TO 62 ! Left Boundary
C
T1 = HX( I)/( HX( I-1)* ( HX( I-1) + HX( I)) )
T2 = ( HX( I) - HX( I-1))/( HX( I-1)* HX( I))
T3 = HX( I-1)/( HX( I)* ( HX( I-1) + HX( I)) )
C
TP1( I,J) = TP1( I,J) + RLN* PRN* DT
C ==============
1 * ( - T1* WP1( I-1,J) + T2* WP1( I,J) + T3* WP1( I+1,J))
C
GO TO 800
C
C ----- Right Boundary ---
C
52 IF ( TPT( I,J,4)) GO TO 800
C
C ---- Perodic right ---
C
TP1( I,J) = TP1( I,J)
C
1 + RLN* PRN* DT* ( WP1( I,J) - WP1( I-1,J))/ HX( I-1)
C
GO TO 800
C
C ---- Left Boundary ---
C
62 IF ( TPT( I,J,4)) GO TO 800 ! Continuative
C
C ----- Periodic left ---
C
TP1( I,J) = TP1( I,J)
1 + RLN* PRN* DT* ( WP1( I+1,J) - WP1( I,J))/ HX( I)
C
800 CONTINUE
C --------------------------------------
C
900 CONTINUE
C
RETURN
END
C **********************************************************************
C
SUBROUTINE WFCA ( WP1, PS1, WPT, HX, HY, RAS, RAF,
1 IX, IY, IX1, 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 /UDT/ DY2, ITN
COMMON /CTR/ KTF, MSH, IHS, WEI
C
C ----- No-slip surfaces : Updated ---
C
C -----------------------------------------------------
DO 100 I = 1, IX
DO 100 J = 1, IY
C -----------------------------------------------------
C
IF ( WPT( I,J,1)) GO TO 100 ! Fluid
IF ( WPT( I,J,2)) GO TO 50 ! Slip
IF ( WPT( I,J,3)) GO TO 200 ! Moving
C
C ----- Attention ---
C
IF ( J .EQ.IY) GO TO 111
IF ( J .EQ. 1) GO TO 222
IF ( I .EQ.IX) GO TO 333
IF ( I .EQ. 1) GO TO 444
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
WD = ( PS1( I,J) - PS1( I,J-1))/ RAS( J-1)
C
114 WP1( I,J) = WA + WB + WC + WD
GO TO 100
C
C ----- Boundary ---
C
1 WA = ( PS1( I,J) - PS1( I-1,J))/ RAF( I-1)
GO TO 13
C
3 WB = ( PS1( I,J) - PS1( I+1,J))/ RAF( I)
GO TO 12
C
2 WC = ( PS1( I,J) - PS1( I,J-1))/ RAS( J-1)
GO TO 14
C
4 WD = ( PS1( I,J) - PS1( I,J+1))/ RAS( J)
GO TO 114
C
C ----- No-slip surfaces at edges ---
C
C ( Upper Boundary )
C
111 CONTINUE
C
IF ( KTF.EQ.1) GO TO 110
C
C ----- Cavity Flow Problem --------------------------------
C
C ( Moving Boundary)
C
WP1( I,J) = 3.D0* ( PS1( I,J) - PS1( I,J-1) - HY( J-1))
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
C ----- Not-Moving Boundary --------------------------------
C
110 CONTINUE
WP1( I,J) = 3.D0* ( PS1( I,J) - PS1( I,J-1))/ RAS( J-1)
1 - 0.5D0* WP1( I,J-1)
C
IF ( I.EQ.1.OR.I.EQ.IX) WP1(I,J) = 0.D0
GO TO 100
C
222 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
GO TO 100
C
333 WP1( I,J) = 3.D0* ( PS1( I,J) - PS1( I-1,J))/ RAF( I-1)
1 - 0.5D0* WP1( I-1,J)
GO TO 100
C
444 WP1( I,J) = 3.D0* ( PS1( I,J) - PS1( I+1,J))/ RAF( I)
1 - 0.5D0* WP1( I+1,J) ! OK
GO TO 100
C
50 WP1( I,J) = 0.D0
GO TO 100
C
200 CONTINUE
C
IF ( J.EQ.IY) GO TO 111
IF ( J.EQ.1 ) GO TO 222
C
100 CONTINUE
C -----------------------------------------------------
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 999
C
WRITE(6,*) ' * IX.NE.IY * STOP * '
WRITE(8,*) ' * IX.NE.IY * STOP * '
STOP
C
999 NN = 0
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 -------------------------------------------
DO 200 J = 1, IY-1
DO 200 I = 1, IX-1
C
NM = ( J-1)* ( IX-1) + I
NDE( NM,1) = ( J-1)* IX + I
NDE( NM,2) = ( J-1)* IX + I + 1
NDE( NM,3) = ( J-1)* IX + I + IX + 1
NDE( NM,4) = ( J-1)* IX + I + IX
C
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
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 SEKB ( SKN, U, TP1, IX, IY, IX1, HX, GG, XZ, YZ )
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION SKN( IX,IY), U ( IX,IY), TP1( IX,IY), GG( IX,IY),
1 HX ( IX1), XZ( IX), YZ ( IY)
C
C ------------------------------------------------
DO 100 I = 1, IX
DO 100 J = 1, IY
C ------------------------------------------------
C
IF ( I.NE.1) GO TO 20
C
C ----- Coefficient B ( BB11) ---
C
B111 = ( TP1( I+2,J)* XZ( I)** 2 - TP1( I,J)* XZ( I+2)** 2 )
1 /( XZ( I)** 2 - XZ( I+2)** 2)
C
B122 = ( TP1( I+2,J)* XZ( I+1)** 2 - TP1( I+1,J)* XZ( I+2)** 2 )
1 /( XZ( I+1)** 2 - XZ( I+2)** 2 )
C
B133 = ( XZ( I+2)* XZ( I)** 2 - XZ( I)* XZ( I+2)** 2 )
1 /( XZ( I)** 2 - XZ( I+2)** 2 )
C
B144 = ( XZ( I+2)* XZ( I+1)** 2 - XZ( I+1)* XZ( I+2)** 2 )
1 /( XZ( I+1)** 2 - XZ( I+2)** 2 )
C
BB19 = ( B111 - B122)/( B133 - B144)
C
GG( I,J) = BB19
GO TO 40
C
20 IF ( I.NE.IX) GO TO 30
C
B112 = ( TP1( I,J)* XZ( I-2)** 2 - TP1( I-2,J)* XZ( I)** 2 )
1 /( XZ( I-2)** 2 - XZ( I)** 2 )
C
B222 = ( TP1( I,J)* XZ( I-1)** 2 - TP1( I-1,J)* XZ( I)** 2 )
1 /( XZ( I-1)** 2 - XZ( I)** 2 )
C
B132 = ( XZ( I)* XZ( I-2)** 2 - XZ( I-2)* XZ( I)** 2 )
1 /( XZ( I-2)** 2 - XZ( I)** 2 )
C
B142 = ( XZ( I)* XZ( I-1)** 2 - XZ( I-1)* XZ( I)** 2)
1 /( XZ( I-1)** 2 - XZ( I)** 2)
C
BB18 = ( B112- B222)/( B132- B142)
C
C ----------------------------------------------------------
C
C113 = ( TP1( I,J)* XZ( I-2) - TP1( I-2,J)* XZ( I))
1 /( XZ( I-2) - XZ( I))
C
C123 = ( TP1( I,J)* XZ( I-1) - TP1( I-1,J)* XZ( I))
1 /( XZ( I-1) - XZ( I))
C
C133 = ( XZ( I-2)* XZ( I)** 2 - XZ( I)* XZ( I-2)** 2)
1 /( XZ( I-2) - XZ( I))
C
C143 = ( XZ( I-1)* XZ( I)** 2 - XZ( I)* XZ( I-1)** 2)
1 /( XZ( I-1) - XZ( I))
C
CC17 = ( C113 - C123)/( C133 - C143)
C
GG( I,J) = BB18 + 2.D0 * CC17
GO TO 40
C
30 CONTINUE
C
G1 = HX( I)/( HX( I-1)* ( HX( I-1) + HX( I)) )
G2 = ( HX( I) - HX( I-1))/( HX( I-1)* HX( I))
G3 = HX( I-1)/( HX( I)* ( HX( I-1) + HX( I)) )
C
GG( I,J) = - G1* TP1( I-1,J) + G2* TP1( I,J) + G3* TP1( I+1,J)
C
40 SKN( I,J) = U( I,J)* TP1( I,J) - GG( I,J)
C
100 CONTINUE
C ------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE NUMX ( AMX, YMX, AMN, YMN, SKN, IX, IY, HX, HY,
1 IX1, IY1, I, XZ, YZ )
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION SKN( IX,IY), HX( IX1), HY( IY1), XZ( IX), YZ( IY)
C
C --------------------------------------
AMX = 0.D0
DO 100 J = 1, IY
C
IF ( SKN( I,J).LE.AMX ) GO TO 100
AMX = SKN( I,J)
JMXS = J
C
100 CONTINUE
C ---------------------------------------
AMN = AMX
DO 200 J = 1, IY
C
IF ( SKN( I,J).GE.AMN ) GO TO 200
AMN = SKN( I,J)
JMIN = J
C
200 CONTINUE
C ---------------------------------------
YMX = YZ( JMXS)
YMN = YZ( JMIN )
C
RETURN
END
C **********************************************************************
C
SUBROUTINE SIMP ( RNU, SKN, SIM, ANU, IX, IY, IX1, IY1, HX,
1 HY, XZ, YZ )
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION SKN( IX,IY), SIM( IX), ANU( IX), HX( IX-1),
1 HY ( IY-1), XZ ( IX), YZ ( IY)
C
MMX = ( IX-1)/2
C
C --------------------------------------
DO 100 I = 1, IX
DO 100 J = 1, IY-1
C
SIM( I) = ( SKN( I,J) + SKN( I,J+1))* HY( J)/ 2.D0
ANU( I) = ANU( I) + SIM( I)
C
100 CONTINUE
C --------------------------------------
SMP = 0.D0
RNU = 0.D0
DO 200 I = 1, IX-1
C
SMP = ( ANU( I) + ANU( I+1))* HX( I)/ 2.D0
RNU = RNU + SMP
C
200 CONTINUE
C --------------------------------------
C
C ----- Attention : Not used, so far ---
C
SUG = 0.D0
C
DO 300 I = 2, IX-3, 2
C
C ----- Coefficient : A ( AA11) ---
C
A11 = ( ANU( I)*YZ( I+2)*YZ( I+2) - ANU( I+2)*YZ( I)*YZ( I))
1 /( YZ( I+2)*YZ( I+2) - YZ( I+2)*YZ( I))
C
A12 = ( ANU( I)*YZ( I+1)** 2 - ANU( I+1)*YZ( I)** 2)
1 /( YZ( I+1)** 2 - YZ( I+1)*YZ( I))
C
A13 = ( YZ( I+2)** 2 - YZ( I)** 2 )
1 /( YZ( I+2)** 2 - YZ( I+2)* YZ( I))
C
A14 = ( YZ( I+1)** 2 - YZ( I)** 2)
1 /( YZ( I+1)** 2 - YZ( I+1)* YZ( I))
C
AA11 = ( A11 - A12)/( A13 - A14)
C
C ----- Coefficient : B ( BB11) ---
C
B11 = ( ANU( I+2)*YZ( I)** 2 - ANU( I)*YZ( I+2)** 2)
1 /( YZ( I)** 2 - YZ( I+2)** 2)
C
B12 = ( ANU( I+2) * YZ( I+1)** 2 - ANU( I+1) * YZ( I+2)** 2 )
1 /( YZ( I+1)** 2 - YZ( I+2)** 2)
C
B13 = ( YZ( I+2)*YZ( I)** 2 - YZ( I)*YZ( I+2)** 2)
1 /( YZ( I)** 2 - YZ( I+2)** 2)
C
B14 = ( YZ( I+2)*YZ( I+1)** 2 - YZ( I+1)*YZ( I+2)** 2 )
1 /( YZ( I+1)** 2 - YZ( I+2)** 2)
C
BB11 = ( B11- B12)/( B13- B14)
C
C ----- Coefficient : C ( CC11) ---
C
C11 = ( ANU( I+2)*YZ( I) - ANU( I)*YZ( I+2))
1 /( YZ( I) - YZ( I+2))
C
C12 = ( ANU( I+2)*YZ( I+1) - ANU( I+1)*YZ( I+2))
1 /( YZ( I+1) - YZ( I+2))
C
C13 = ( YZ( I)*YZ( I+2)** 2 - YZ( I+2)*YZ( I)** 2)
1 /( YZ( I) - YZ( I+2))
C
C14 = ( YZ( I+1)*YZ( I+2)** 2 - YZ( I+2)*YZ( I+1)** 2)
1 /( YZ( I+1) - YZ( I+2))
C
CC11 = ( C11 - C12)/( C13 - C14)
C
C ----- General Expression ---
C
HAN = AA11* YZ( I+2) + BB11* YZ( I+2)** 2/ 2.D0
1 + CC11* YZ( I+2)** 3/ 3.D0
C
2 - AA11* YZ( I) - BB11* YZ( I)** 2 / 2.D0
3 - CC11* YZ( I)** 3/ 3.D0
C
SUG = SUG + HAN
C
300 CONTINUE
C --------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE SMPS ( SKN, SIM, ANU, IX,IY, IX1, IY1, HX, HY )
C
IMPLICIT REAL*8(A-H,O-Z)
C
COMMON /ABC/ DNU, RN0, RN1
C
DIMENSION SKN( IX,IY), SIM( IX), ANU( IX), HX( IX1), HY( IY1)
C
C ------------------------------------------------
DO 100 I = 1, IX
C
C --------------------------------------
DO 150 J = 1, IY-1
C
SIM( I) = ( SKN( I,J) + SKN( I,J+1))* HY( J)/ 2.D0
ANU( I) = ANU( I) + SIM( I)
C
150 CONTINUE
C --------------------------------------
IF ( I.NE. 1) GO TO 170
RN0 = ANU( I)
170 IF ( I.NE.IX) GO TO 100
RN1 = ANU( I)
C
100 CONTINUE
C ------------------------------------------------
DNU = DABS( RN0 - RN1 )
C
RETURN
END
C **********************************************************************
C
SUBROUTINE CLNU ( SKN, GG, U, V, HX, HY, SIM, ANU, TP1,
1 XZ, YZ, IX, IY, IX1, IY1)
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION U ( IX,IY), V ( IX,IY), TP1( IX,IY), SKN( IX,IY),
1 GG( IX,IY), SIM( IX), ANU( IX), HX ( IX1),
2 HY( IY1), XZ ( IX), YZ ( IY)
C
C ----------------------------
DO 100 I = 1, IX
C
U( I, 1) = 0.D0
V( I, 1) = 0.D0
U( I,IY) = 0.D0
V( I,IY) = 0.D0
C
100 CONTINUE
C ----------------------------
DO 200 J = 1, IY
C
U( 1, J) = 0.D0
V( 1, J) = 0.D0
U( IX,J) = 0.D0
V( IX,J) = 0.D0
C
200 CONTINUE
C --------------------------------------
DO 300 I = 1, IX
C
C ----------------------------
DO 350 J = 1, IY
C
SKN( I,J) = 0.D0
GG ( I,J) = 0.D0
C
350 CONTINUE
C ----------------------------
SIM( I) = 0.D0
ANU( I) = 0.D0
C
300 CONTINUE
C --------------------------------------
C
C ----- Prep. for Numerical Integration ---------------------------
C
CALL SEKB ( SKN, U, TP1, IX, IY, IX1, HX, GG, XZ, YZ )
C
C ----- Average Nu ------------------------------------------------
C
CALL SMPS ( SKN, SIM, ANU, IX, IY, IX1, IY1, HX, HY )
C
C -----------------------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE CMNU ( U, TP1, IX, IY, IX1, IY1, SKN, GG,
1 SIM, ANU, HX, HY, XZ, YZ )
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION U ( IX,IY), TP1( IX,IY), SKN( IX,IY), GG( IX,IY),
1 SIM( IX), ANU( IX), HX ( IX1), HY( IY1),
2 XZ ( IX), YZ ( IY)
C
C ------------------------------------------------
DO 100 I = 1, IX
C
C --------------------------------------
DO 150 J = 1, IY
C
SKN( I,J) = 0.D0
GG ( I,J) = 0.D0
IF ( I.EQ.1) U( I,J) = 0.D0
C
150 CONTINUE
C --------------------------------------
SIM( I) = 0.D0
ANU( I) = 0.D0
C
100 CONTINUE
C ------------------------------------------------
C
C ----- Numerical Integration ---------------------------------------
C
CALL SEKB ( SKN, U, TP1, IX, IY, IX1, HX, GG, XZ, YZ )
C
C ----- Average Nu for all ( X,Y) -----------------------------------
C
CALL SIMP ( RNU, SKN, SIM, ANU, IX, IY, IX1, IY1, HX,
1 HY, XZ, YZ )
C
C ----- NuMax.& NuMin.on X = 0 ---------------------------------------
C
CALL NUMX ( AMX, YMX, AMN, YMN, SKN, IX, IY, HX, HY,
1 IX1, IY1, 1, XZ, YZ )
C
C --------------------------------------------------------------------
WRITE(6,2000) RNU, AMX, YMX, AMN, YMN
WRITE(8,2000) RNU, AMX, YMX, AMN, YMN
2000 FORMAT(/,' * AvNu=',F7.3,': At X=0',' N_mx=',F7.3,' at Y=',F6.3,
1 ' N_mn=',F7.3,' at Y=',F6.3)
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( 101,101), PA2( 101,101), PA3( 101,101),
1 PB1( 101,101), PB2( 101,101), PB3( 101,101)
C
DIMENSION HX( IX1), HY( IY1)
C
C --------------------------------------
DO 100 I = 1, IX
DO 100 J = 1, IY
C
PA1( I,J) = 0.D0
PA2( I,J) = 0.D0
PA3( I,J) = 0.D0
PB1( I,J) = 0.D0
PB2( I,J) = 0.D0
PB3( I,J) = 0.D0
C
100 CONTINUE
C --------------------------------------
DO 200 I = 2, IX-1
DO 200 J = 2, IY-1
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)) )
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,
1 UND, UNE, UAP, UBP, UCP, UDP, UEP, UAM,
2 UBM, UCM, 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 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
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 I = 1, IX
DO 200 K = 1, 4
C --------------------------------------
C
UNA( I,K) = 0.D0
UNB( I,K) = 0.D0
UNC( I,K) = 0.D0
UND( I,K) = 0.D0
UNE( I,K) = 0.D0
C
UAP( I,K) = 0.D0
UBP( I,K) = 0.D0
UCP( I,K) = 0.D0
UDP( I,K) = 0.D0
UEP( I,K) = 0.D0
C
UAM( I,K) = 0.D0
UBM( I,K) = 0.D0
UCM( I,K) = 0.D0
UDM( I,K) = 0.D0
UEM( I,K) = 0.D0
C
200 CONTINUE
C --------------------------------------
DO 300 J = 1, IY
C
UF1( J) = 0.D0
UF2( J) = 0.D0
UF3( J) = 0.D0
C
300 CONTINUE
C --------------------------------------
DO 400 J = 2, IY-1
C
UF1( J) = - HY( J)/( HY( J-1)* ( HY( J-1) + HY( J)) )
UF2( J) = ( HY( J) - HY( J-1))/( HY( J-1)* HY( J))
UF3( J) = HY( J-1)/( HY( J)* ( HY( J-1) + HY( J)) )
C
400 CONTINUE
C ------------------------------------------------
DO 500 I = 3, IX-2
C
HA = HX(I-2)* ( HX( I-2) + HX( I-1))
1 * ( HX( I-2) + HX( I-1) + HX( I))
2 * ( HX( I-2) + HX( I-1) + HX( I) + HX( I+1))
C
HB = HX( I-1)*( HX( I-1) + HX( I))*( HX( I-1) + HX( I) + HX( I+1))
C
HD = HX( I)* HX( I+1)
HE = HX( I) + HX( I+1)
C
C --------------------------------------
DO 550 K = 1, 4
C
UNA( I,K) = ( BB( K)* HX( I-1)* HX( I)* ( HX( I) + HX( I+1))
1 - CC( K)* ( HX( I-1)* ( 2.D0* HX( I) + HX( I+1))
2 - HX( I)* ( HX( I) + HX( I+1)) )
3 + DD( K)* ( HX( I-1) - 2.D0* HX( I) - HX( I+1))
4 + EE( K))/ HA
C
UNB( I,K) = ( - UNA( I,K)* ( HX( I-2) + HX( I-1))
1 * ( HX( I-2) + HX( I-1) + HX( I))
2 * ( HX( I-2) + HX( I-1) + HX( I) + HX( I+1))
3 - BB( K)* HX( I)* ( HX( I) + HX( I+1))
4 + CC( K)* ( 2.D0* HX( I) + HX( I+1)) - DD( K))/ HB
C
UND( I,K) = ( UNA( I,K)* ( HX( I-2) + HX( I-1))
1 * ( HX( I-2) + HX( I-1) + HX( I) + HX( I+1))
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
550 CONTINUE
C --------------------------------------
C
500 CONTINUE
C ------------------------------------------------
DO 600 I = 4, IX-1
C
HA = HX( I-3)* ( HX( I-3) + HX( I-2))
1 * ( HX( I-3) + HX( I-2) + HX( I-1))
2 * ( HX( I-3) + HX( I-2) + HX( I-1) + HX( I))
C
HB = HX( I-2)* ( HX( I-2)+ HX( I-1))*( HX( I-2)+ HX( I-1)+ HX( I))
C
HC = HX( I-1)* ( HX( I-1) + HX( I))
HE = HX( I)
C
C -------------------------------------
DO 650 K = 1, 4
C
UAP( I,K) = ( - BB( K)* HX( I-1)* HX( I)* ( HX( I-2) + HX( I-1))
1 + CC( K)* ( HX( I-2)* HX( I-1) + ( HX( I-1)** 2)
2 - HX( I-2)* HX( I) - 2.D0* HX( I-1)* HX( I))
3 + DD( K)* ( HX( I-2) + 2.D0* HX( I-1) - HX( I))
4 + 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
650 CONTINUE
C --------------------------------------
C
600 CONTINUE
C ------------------------------------------------
DO 700 I = 2, IX-3
C
HE = HX( I+2)* ( HX( I+1)+HX( I+2))* ( HX( I)+HX( I+1)+HX( I+2))
1 * ( 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))
HA = HX( I-1)
C
C --------------------------------------
DO 750 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))
3 + DD( K)* ( HX( I-1) - 2.D0* HX( I) - HX( I+1))
4 + 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)
4 - CC( K)* ( HX( I) - HX( I-1)) + DD( K))/ HD
C
UCM( I,K) = ( - UEM( I,K)* ( HX( I) + HX( I+1) + HX( I+2))
1 * ( HX( I-1) + HX( I) + HX( I+1) + HX( I+2))
C
2 - UDM( I,K) * ( HX( I)+HX( I+1)) * ( HX( I-1)+HX( I)+HX( I+1))
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
750 CONTINUE
C --------------------------------------
C
700 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, LTC )
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), AA ( 4), BB ( 4),
2 CC ( 4), DD ( 4), EE ( 4), UPA( 4),
3 UPB( 4), UPC( 4), UPD( 4), UF1( IY),
4 UF2( IY), UF3( IY), UNA( IX,4), UNB( IX,4),
5 UNC( IX,4), UND( IX,4), UNE( IX,4), UAP( IX,4),
6 UBP( IX,4), UCP( IX,4), UDP( IX,4), UEP( IX,4),
7 UAM( IX,4), UBM( IX,4), UCM( IX,4), UDM( IX,4),
8 UEM( IX,4)
C
LOGICAL *1 WPT( IX,IY,4), EPT( IX,IY,4)
C
COMMON /UDT/ DY2, ITN
COMMON /CTR/ KTF, MSH, IHS, WEI
COMMON /IDT/ DT, LPM, PRN, CTP, RLN, REN, LOP, NPL, NPR
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
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 I = 1, IX
DO 200 J = 1, IY
C ----------------------------------------------------------
C
IF ( EPT( I,J,3)) GO TO 14 ! Right
IF (.NOT.EPT( I,J,1)) GO TO 5
IF (.NOT.WPT( I,J,4)) GO TO 15 ! Periodic
RUF( I,J) = 0.D0
GO TO 200
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 200
C
10 CONTINUE
RUF(I,J)
1 = UF1( J)* PS1(I,J-1) + UF2( J)* PS1(I,J) + UF3( J)* PS1(I,J+1)
C ======================================================================
C
14 IF (.NOT.EPT( I,J,1)) GO TO 112
15 IF (.NOT.WPT( I,J,2)) GO TO 16 ! No slip
C
C ----- Slip ---
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
C ===============
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
C ==============
RUF( I,J) = UL ! Not moving
GO TO 115
C
22 RUF( I,J) = 0.D0
GO TO 115
C
113 CONTINUE
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 ! Periodic
C
C ----- Continuative ---
C
RUF( I,J) = UF1( J)* PS1( I,J-1) + UF2( J)* PS1( I,J)
1 + UF3( J)* PS1( I,J+1)
C ==========================================================
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 ! No slip
RUF( I,J) = 0.D0
GO TO 115
C
31 IF (.NOT.WPT( I,J,3)) GO TO 32 ! Not moving
RUF( I,J) = 0.D0
GO TO 115
C
32 IF (.NOT.WPT( I,J,4)) GO TO 33 ! Periodic
GO TO 115
C
33 CONTINUE
C
115 CONTINUE
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
GO TO 80
C
29 CONTINUE
A1 = HX( I)/( HX( I-1)* ( HX( I-1) + HX( I)) )
B1 = ( HX( I-1) - HX( I))/( HX( I-1)* HX( I))
C1 = HX( I-1)/( HX( I)* ( HX( I-1) + HX( I)) )
C
A2 = 1.D0 /( HX( I-1)* ( HX( I-1) + HX( I)) )
B2 = 1.D0 /( HX( I-1)* HX( I))
C2 = 1.D0 /( HX( I) * ( HX( I-1) + HX( I)) )
C
CENT = - UT* ( - A1* WM1( I-1,J) - B1* WM1( I,J) + C1* WM1(I+1,J))
1 + ( UT** 2 )*( A2* WM1( I-1,J) - B2* WM1( I,J) + C2* WM1(I+1,J))
C
80 CONTINUE
IF ( 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
UPW = - UT* ( UAP( I,1)* WM1( I-3, J) + UBP( I,1)* WM1( I-2,J)
1 + UCP( I,1)* WM1( I-1,J) + UDP( I,1)* WM1( I,J)
2 + UEP( I,1)* WM1( I+1,J))
C
3 + (UT** 2)* ( UAP( I,2)* WM1( I-3, J) + UBP( I,2)* WM1( I-2,J)
4 + UCP( I,2)* WM1( I-1,J) + UDP( I,2)* WM1( I,J)
5 + UEP( I,2)* WM1( I+1,J))/ 2.D0
C
6 - (UT** 3)* ( UAP( I,3)* WM1( I-3,J) + UBP( I,3)* WM1( I-2,J)
7 + UCP( I,3)* WM1( I-1,J) + UDP( I,3)* WM1( I,J)
8 + UEP( I,3)* WM1( I+1,J))/ 6.D0
C
9 + (UT** 4)* ( UAP( I,4)* WM1( I-3,J) + UBP( I,4)* WM1( I-2,J)
A + UCP( I,4)* WM1( I-1,J) + UDP( I,4)* WM1( I,J)
B + UEP( I,4)* WM1( I+1,J))/ 24.D0
GO TO 150
C
333 WP1( I,J) = WP1( I,J) - UT * WM1( I,J)/ HX( I)
GO TO 200
C
110 UPW = - UT * ( - WM1( I-1,J) + WM1( I,J))/ HX( I-1)
GO TO 150
C
120 CONTINUE
C
C -------------------------------------------------
DO 250 K = 1, 3
C
HA = HX( I-2)*( HX( I-2)+ HX( I-1))*( HX( I-2)+ HX( I-1)+ HX( I))
C
HB = HX( I-1)* ( HX( I-1) + HX( I))
HD = HX( I)
C
UPA( K) = ( BB( K)* HX( I-1)* HX( I)
1 + CC( K)* ( HX( I) - HX( I-1)) - DD( K))/ HA
C
UPB( K) = ( - UPA( K)* ( HX( I-2) + HX( I-1))
1 * ( HX( I-2) + HX( I-1) + HX( I))
2 - BB( K)* HX( I) + CC( K))/ HB
C
UPD( K) = ( UPA( K)* ( HX( I-2) + HX( I-1))
1 + UPB( K)* HX( I-1) + BB( K))/ HD
C
UPC( K) = AA( K) - ( UPA( K) + UPB( K) + UPD( K))
C
250 CONTINUE
C ------------------------------------------------
UPW = - UT* ( UPA( 1)* WM1( I-2,J) + UPB( 1)* WM1( I-1,J)
1 + UPC( 1)* WM1( I,J) + UPD( 1)* WM1( I+1,J))
C
2 +( UT** 2)* ( UPA( 2)* WM1( I-2,J) + UPB( 2)* WM1( I-1,J)
3 + UPC( 2)* WM1( I,J) + UPD( 2)* WM1( I+1,J))/ 2.D0
C
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
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))
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
UPW = - UT* ( A1* WM1( I-2,J) - B1* WM1( I-1,J) + C1* WM1( I,J))
1 + ( UT** 2)* ( A2* WM1( I-2,J) - B2* WM1( I-1,J) + C2* WM1( I,J))
C
WP1( I,J) = WP1( I,J) + UPW
GO TO 200
C
90 IF ( EPT( I,J,4)) GO TO 170
IF ( EPT( I,J,3)) GO TO 550
IF ( EPT( I+1,J,3)) GO TO 190
IF ( EPT( I+2,J,3)) GO TO 180
C
UPW = - UT* ( UAM( I,1)* WM1( I-1,J) + UBM( I,1)* WM1( I,J)
1 + UCM( I,1)* WM1( I+1,J) + UDM( I,1)* WM1( I+2,J)
2 + UEM( I,1)* WM1( I+3,J))
C
3 +( UT** 2)* ( UAM( I,2)* WM1( I-1,J) + UBM( I,2)* WM1( I,J)
4 + UCM( I,2)* WM1( I+1,J) + UDM( I,2)* WM1( I+2,J)
5 + UEM( I,2)* WM1( I+3, J))/ 2.D0
C
6 -( UT** 3)* ( UAM( I,3)* WM1( I-1,J) + UBM( I,3)* WM1( I,J)
7 + UCM( I,3)* WM1( I+1,J) + UDM( I,3)* WM1( I+2,J)
8 + UEM( I,3)* WM1( I+3,J))/ 6.D0
C
9 +( UT** 4)* ( UAM( I,4)* WM1( I-1,J) + UBM( I,4)* WM1( I,J)
A + UCM( I,4)* WM1( I+1,J) + UDM( I,4)* WM1( I+2,J)
B + UEM( I,4)* WM1( I+3,J))/ 24.D0
GO TO 220
C
170 A1 = ( 2.D0* HX( I) + HX( I+1))/( HX( I)* ( HX( I) + HX( I+1)) )
B1 = ( HX( I) + HX( I+1))/( HX( I)* HX( I+1))
C1 = HX( I)/( HX( I+1)* ( HX( I) + HX( I+1)) )
C
A2 = 1.D0 /( HX( I)* ( HX( I) + HX( I+1)) )
B2 = 1.D0 /( HX( I)* HX( I+1))
C2 = 1.D0 /( HX( I+1)* ( HX( I) + HX( I+1)) )
C
UPW = - UT* (- A1* WM1( I,J) + B1* WM1( I+1,J) - C1* WM1( I+2,J))
1 + ( UT** 2)*( A2* WM1( I,J) - B2* WM1( I+1,J) + C2* WM1( I+2,J))
C
WP1( I,J) = WP1( I,J) + UPW
GO TO 200
C
180 CONTINUE
C
C ------------------------------------------------
DO 230 K = 1, 3
C
HD = HX( I+1) * ( HX( I)+HX( I+1)) * ( HX( I-1)+HX( I)+HX( I+1))
HC = HX( I) * ( HX( I-1)+HX( I))
HA = HX( I-1)
C
UPD( K) = ( - BB( K)* HX( I-1)* HX( I)
1 - CC( K)* ( HX( I)-HX( I-1)) + DD( K))/ HD
C
UPC( K) = ( - UPD( K)* ( HX( I) + HX( I+1))
1 * ( HX( I-1) + HX( I) + HX( I+1))
2 + BB( K)* HX( I-1) + CC( K))/ HC
C
UPA( K) = ( UPD( K)* ( HX( I) + HX( I+1))
1 + UPC( K)* HX( I) - BB( K))/ HA
C
UPB( K) = AA( K) - ( UPA( K) + UPC( K) + UPD( K))
C
230 CONTINUE
C ------------------------------------------------
C
UPW = - UT* ( UPA( 1)* WM1( I-1,J) + 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 UPW = - UT* ( - WM1( I,J) + WM1( I+1,J))/ HX( I)
GO TO 220
C
550 WP1( I,J) = WP1( I,J) + UT* WM1( I,J)/ HX( I-1)
GO TO 200
C
220 CONTINUE
C
150 CONTINUE
C
WP1( I,J) = WP1( I,J) + ( 1.D0 - WEI )* CENT + WEI* UPW
C =============================================================
C
200 CONTINUE
C ----------------------------------------------------------
C
IF ( LTC.NE.1) RETURN
C
C ----------------------------------------------------------
DO 300 J = 1, IY
DO 300 I = 1, IX
C ----------------------------------------------------------
C
IF (.NOT.WPT( I,J,4)) GO TO 300 ! Periodic
IF ( EPT( I,J,1)) WP1( I,J) = WP1( I,J-1)
IF ( EPT( I,J,2)) WP1( I,J) = WP1( I,J+1)
IF ( EPT( I,J,3)) WP1( I,J) = WP1( I-1,J)
IF ( EPT( I,J,4)) WP1( I,J) = WP1( I+1,J)
C
300 CONTINUE
C ----------------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE WVCF ( HX, HY, VF1, VF2, VF3, UNA, UNB, UNC,
1 UND, UNE, UAP, UBP, UCP, UDP, UEP, UAM,
2 UBM, UCM, 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), VF1( IX),
2 VF2( IX), VF3( IX), UNA( IY,4), UNB( IY,4),
3 UNC( IY,4), UND( IY,4), UNE( IY,4), UAP( IY,4),
4 UBP( IY,4), UCP( IY,4), UDP( IY,4), UEP( IY,4),
5 UAM( IY,4), UBM( IY,4), UCM( IY,4), UDM( IY,4),
6 UEM( 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
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 --------------------------------------
UNA( J,K) = 0.D0
UNB( J,K) = 0.D0
UNC( J,K) = 0.D0
UND( J,K) = 0.D0
UNE( J,K) = 0.D0
C
UAP( J,K) = 0.D0
UBP( J,K) = 0.D0
UCP( J,K) = 0.D0
UDP( J,K) = 0.D0
UEP( J,K) = 0.D0
C
UAM( J,K) = 0.D0
UBM( J,K) = 0.D0
UCM( J,K) = 0.D0
UDM( J,K) = 0.D0
UEM( 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
UNA( 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)) )
3 + DD( K)* ( HY( J-1) - 2.D0* HY( J) - HY( J+1))
4 + EE( K))/ HA
C
UNB( J,K) = ( - UNA( 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))
5 - DD( K))/ HB
C
UND( J,K) = ( UNA( J,K) * ( HY( J-2) + HY( J-1))
1 * ( HY( J-2) + HY( J-1) + HY( J) + HY( J+1))
2 + UNB( J,K) * HY( J-1) * ( HY( J-1) + HY( J) + HY( J+1))
3 + BB( K) * ( HY( J) + HY( J+1)) - CC( K))/ HD
C
UNE( J,K) = ( UNA( J,K)* ( HY( J-2) + HY( J-1))
1 + UNB( J,K)* HY( J-1)
2 - UND( J,K)* HY( J) + BB( K))/ HE
C
UNC( J,K)= AA( K)
1 -( UNA( J,K) + UNB( J,K) + UND( J,K) + UNE( J,K))
C
700 CONTINUE
C --------------------------------------
C
500 CONTINUE
C ------------------------------------------------
DO 800 J = 4, IY-1
C
HA = HY( J-3)* ( HY( J-3) + HY( J-2))
1 * ( HY( J-3) + HY( J-2) + HY( J-1))
2 * ( HY( J-3) + HY( J-2) + HY( J-1) + HY( J))
C
HB = HY( J-2)* ( HY( J-2) + HY( J-1))
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
UAP( J,K) = ( - BB( K)* HY( J-1)* HY( J)* ( HY( J-2) + HY( J-1))
1 + CC( K)* ( HY( J-2)* HY( J-1) + HY( J-1)** 2
2 - HY( J-2)* HY( J) - 2.D0* HY( J-1)* HY( J))
C
3 + DD( K)* ( HY( J-2) + 2.D0* HY( J-1) - HY( J)) + EE( K))/ HA
C
UBP( J,K) = ( - UAP( J,K) * ( HY( J-3) + HY( J-2))
1 * ( HY( J-3) + HY( J-2) + HY( J-1))
2 * ( HY( J-3) + HY( J-2) + HY( J-1) + HY( J))
C
3 + BB( K)* HY( J-1)* HY( J)
4 - CC( K)* ( HY( J-1) - HY( J)) - DD( K))/ HB
C
UCP( J,K) = ( - UAP( 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 - UBP( J,K)* ( HY( J-2) + HY( J-1))
3 * ( HY( J-2) + HY( J-1) + HY( J))
4 - BB( K)* HY( J) + CC( K))/ HC
C
UEP( J,K) = ( UAP( J,K)* ( HY( J-3) + HY( J-2) + HY( J-1))
1 + UBP( J,K)* ( HY( J-2) + HY( J-1))
2 + UCP( J,K)* HY( J-1) + BB( K))/ HE
C
UDP( J,K) = AA( K)
1 - ( UAP( J,K) + UBP( J,K) + UCP( J,K) + UEP( J,K))
C
900 CONTINUE
C --------------------------------------
C
800 CONTINUE
C ------------------------------------------------
DO 950 J = 2, IY-3
C
HE = HY( J+2)* ( HY( J+1) + HY( J+2))
1 * ( HY( J) + HY( J+1) + HY( J+2))
2 * ( HY( J-1) + HY( J) + HY( J+1) + HY( J+2))
C
HD = HY( J+1)* ( HY( J) + HY( J+1))
1 * ( HY( J-1) + HY( J) + HY( J+1))
C
HC = HY( J)* ( HY( J-1) + HY( J))
HA = HY( J-1)
C
C --------------------------------------
DO 970 K = 1, 4
C
UEM( 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))
4 + EE( K))/ HE
C
UDM( J,K) = ( - UEM( 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
UCM( J,K) = ( - UEM( J,K)* ( HY( J) + HY( J+1) + HY( J+2))
1 * ( HY( J-1) + HY( J) + HY( J+1) + HY( J+2))
2 - UDM( J,K)* ( HY( J)+HY( J+1))* ( HY( J-1)+HY( J)+HY( J+1))
C
3 + BB( K) * HY( J-1) + CC( K))/ HC
C
UAM( J,K) = ( UEM( J,K)* ( HY( J) + HY( J+1) + HY( J+2))
1 + UDM( J,K)* ( HY( J) + HY( J+1))
2 + UCM( J,K)* HY( J) - BB( K))/ HA
C
UBM( J,K) = AA( K)
1 - ( UAM( J,K) + UCM( J,K) + UDM( J,K) + UEM( J,K))
C
970 CONTINUE
C --------------------------------------
C
950 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, LTC )
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),
5 UNC( IY,4), UND( IY,4), UNE( IY,4), UAP( IY,4),
6 UBP( IY,4), UCP( IY,4), UDP( IY,4), UEP( IY,4),
7 UAM( IY,4), UBM( IY,4), UCM( IY,4), UDM( IY,4),
8 UEM( IY,4)
C
LOGICAL*1 WPT( IX,IY,4), EPT( IX,IY,4)
C
COMMON /UDT/ DY2, ITN
COMMON /CTR/ KTF, MSH, IHS, WEI
COMMON /IDT/ DT, LPM, PRN, CTP, RLN, REN, LOP, NPL, NPR
C
VR = 0.D0
VL = 0.D0
C
C --------------------------------------
DO 100 K = 1, 4
C
AA( K) = 0.D0
BB( K) = 0.D0
CC( K) = 0.D0
DD( K) = 0.D0
EE( K) = 0.D0
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 --------------------------------------
C
C ----------------------------------------------------------
DO 200 J = 1, IY
DO 200 I = 1, IX
C ----------------------------------------------------------
C
IF ( EPT( I,J,1)) GO TO 14
IF (.NOT.EPT( I,J,3)) GO TO 5
IF ( WPT( I,J,4)) GO TO 200
GO TO 25
C
5 IF (.NOT.EPT( I,J,4)) GO TO 10
IF ( WPT( I,J,4)) GO TO 200
GO TO 30
C
10 CONTINUE
RVF( I,J) = - ( VF1( I)*PS1( I-1,J)
1 + VF2( I)*PS1( I,J) + VF3( I)*PS1( I+1,J))
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
RVF( I,J) = 0.D0
GO TO 115
C
17 IF (.NOT.WPT( I,J,4)) GO TO 18
RVF( I,J) = - ( VF1( I)*PS1( I-1,J)
1 + VF2( I)*PS1( I,J) + VF3( I)*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
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
C
RVF( I,J) = RVF( I-1,J)
GO TO 115
C
26 IF (.NOT.WPT( I,J,3)) GO TO 27
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
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
A1 = HY( J)/( HY( J-1)* ( HY( J-1) + HY( J)) )
B1 = ( HY( J-1)-HY( J))/( HY( J-1)* HY( J))
C1 = HY( J-1)/( HY( J)* ( HY( J-1) + HY( J)) )
C
A2 = 1.D0 /( HY( J-1)* ( HY( J-1) + HY( J)) )
B2 = 1.D0 /( HY( J-1)* HY( J))
C2 = 1.D0 /( HY( J)* ( HY( J-1) + HY( J)) )
C
CENT = - VT*( - A1* WM1( I,J-1) - B1* WM1( I,J) + C1* WM1( I,J+1))
1 + ( VT** 2)*(A2* WM1( I,J-1) - B2* WM1( I,J) + C2* WM1( I,J+1))
C
80 CONTINUE
IF ( RVF( I,J).LT.0) GO TO 90
IF ( EPT( I,J,1)) GO TO 130
IF ( EPT( I,J,2)) GO TO 333
IF ( EPT( I,J-1,2)) GO TO 110
IF ( EPT( I,J-2,2)) GO TO 120
C
UPW = - 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
333 WP1( I,J) = WP1( I,J) - VT* WM1( I,J)/ HY( J)
GO TO 200
C
110 UPW = - VT* ( - WM1( I,J-1) + WM1( I,J))/ HY( J-1)
GO TO 150
C
120 CONTINUE
C ------------------------------------------------
DO 250 K = 1, 3
C
HA = HY(J-2)* ( HY( J-2)+ HY( J-1))*( HY( J-2)+ HY( J-1) + HY( J))
C
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
250 CONTINUE
C ------------------------------------------------
C
UPW = - 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
4 - ( VT** 3)* ( UPA( 3)* WM1( I,J-2) + UPB( 3)* WM1( I,J-1)
5 + UPC( 3)* WM1( I,J) + UPD( 3)* WM1( I,J+1))/6.D0
GO TO 150
C
130 A1 = HY( J-1)/( HY( J-2)* ( HY( J-2) + HY( J-1)) )
B1 = ( HY( J-2) + HY( J-1))/( HY( J-2)* HY( J-1))
C1 = ( HY( J-2) + 2.D0* HY( J-1))
1 /( HY( J-1)* ( HY( J-2) + HY( J-1)) )
C
A2 = 1.D0/( HY( J-2)* ( HY( J-2) + HY( J-1)) )
B2 = 1.D0/( HY( J-2)* HY( J-1))
C2 = 1.D0/( HY( J-1)* ( HY( J-2) + HY( J-1)) )
C
UPW = - VT* ( A1* WM1( I,J-2) - B1* WM1( I,J-1) + C1* WM1( I,J))
1 +( VT** 2)* ( A2* WM1( I,J-2) - B2* WM1( I,J-1) + C2* WM1( I,J))
C
WP1( I,J) = WP1( I,J) + UPW
GO TO 200
C
90 IF ( EPT( I, J, 2)) GO TO 170
IF ( EPT( I, J, 1)) GO TO 350
IF ( EPT( I,J+1,1)) GO TO 190
IF ( EPT( I,J+2,1)) GO TO 180
C
UPW = - 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
GO TO 220
C
170 CONTINUE
A1 = ( 2.D0* HY( J) + HY( J+1))/( HY( J)* ( HY( J) + HY( J+1)) )
B1 = ( HY( J) + HY( J+1))/( HY( J)* HY( J+1))
C1 = HY( J)/( HY( J+1)* ( HY( J) + HY( J+1)) )
C
A2 = 1.D0/( HY( J)* ( HY( J) + HY( J+1)) )
B2 = 1.D0/( HY( J)* HY( J+1))
C2 = 1.D0/( HY( J+1)* ( HY( J) + HY( J+1)) )
C
UPW = - VT* (- A1* WM1( I,J) + B1* WM1( I,J+1) - C1* WM1( I,J+2))
1 +( VT** 2)* ( A2* WM1( I,J) - B2* WM1( I,J+1) + C2* WM1( I,J+2))
C
WP1( I,J) = WP1( I,J) + UPW
GO TO 200
C
180 CONTINUE
C
C ------------------------------------------------
DO 300 K = 1, 3
C
HD = HY( J+1)*( HY( J) + HY( J+1))*( HY( J-1) + HY( J) + HY( J+1))
C
HC = HY( J)* ( HY( J-1) + HY( J))
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
300 CONTINUE
C ------------------------------------------------
UPW = - VT* ( UPA( 1)* WM1( I,J-1) + UPB( 1)* WM1( I,J)
1 + UPC( 1)* WM1( I,J+1) + UPD( 1)* WM1( I,J+2))
C
2 + (VT** 2)* ( UPA( 2)* WM1( I,J-1) + UPB( 2)* WM1( I,J)
3 + UPC( 2)* WM1( I,J+1) + UPD( 2)* WM1( I,J+2))/ 2.D0
C
4 - (VT**3)* ( UPA( 3)* WM1( I,J-1) + UPB( 3)* WM1( I,J)
5 + UPC( 3)* WM1( I,J+1) + UPD( 3)* WM1( I,J+2))/ 6.D0
GO TO 220
C
190 UPW = - VT* ( - WM1( I,J) + WM1( I,J+1))/ HY( J)
GO TO 220
C
350 WP1( I,J) = WP1( I,J) + VT* WM1( I,J)/ HY( J-1)
GO TO 200
C
220 CONTINUE
C
150 CONTINUE
C
WP1( I,J) = WP1( I,J) + ( 1.D0 - WEI)* CENT + WEI* UPW
C
200 CONTINUE
C ----------------------------------------------------------
C
IF ( LTC.NE.1) GO TO 900
C ------------------------------------------------
DO 400 J = 1, IY
DO 400 I = 1, IX
C
IF (.NOT.WPT( I,J,4)) GO TO 400 ! Periodic
C
C ----- Continuative ---
C
IF ( EPT( I,J,1)) WP1( I,J) = WP1( I,J-1) ! Upper
IF ( EPT( I,J,2)) WP1( I,J) = WP1( I,J+1) ! Lower
IF ( EPT( I,J,3)) WP1( I,J) = WP1( I-1,J) ! Right
IF ( EPT( I,J,4)) WP1( I,J) = WP1( I+1,J) ! Left
C
400 CONTINUE
C ------------------------------------------------
C
900 CONTINUE
C
RETURN
END
C **********************************************************************
C
SUBROUTINE TCOF ( HX, HY, IX, IY, IX1, IY1)
C
IMPLICIT REAL*8(A-H,O-Z)
C
COMMON /CM2/ TM1( 101), TM2( 101), TM3( 101), TY1( 101),
1 TY2( 101), TY3( 101)
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 **********************************************************************