戻る

C
PROGRAM MAIN
C
C ******************************************************************** C
C C
C Graphics Software for F3D-FEM-R C
C C
C ( Rectangular Parallelepiped Elements ) C
C C
C ( Version: 3.8 ) C
C C
C 3D Fluid Analysis Results C
C C
C Copyright : Yasuhiro MATSUDA C BR> C C
C ******************************************************************** C
C
IMPLICIT REAL*8 (A-H,O-Z)
C
C ***** DATA ( 1 )******************************************************
C
*** PARAMETER ( NEX = 2, NEY = 2, NEZ = 2,
C
PARAMETER ( NEX = 3, NEY = 3, NEZ = 3,
C
*** PARAMETER ( NEX = 4, NEY = 4, NEZ = 4,
*** PARAMETER ( NEX = 5, NEY = 5, NEZ = 5,
*** PARAMETER ( NEX = 6, NEY = 6, NEZ = 6,
*** PARAMETER ( NEX = 10, NEY = 10, NEZ = 10,
*** PARAMETER ( NEX = 20, NEY = 20, NEZ = 20,
*** PARAMETER ( NEX = 30, NEY = 30, NEZ = 30,
C
C **********************************************************************
C
1 LX = NEX + 1, LY = NEY + 1, LZ = NEZ + 1,
2 NNP = LX* LY* LZ, NXY = NEX* NEY, NYZ = NEY* NEZ,
3 NZX = NEZ* NEX )
C
DIMENSION X( NNP), Y( NNP), Z ( NNP), U ( NNP), V( NNP),
1 W( NNP), P( NNP), XLX( LX), YLY( LY)
C
DIMENSION LXY( NXY,4), LYZ( NYZ,4), LZX( NZX,4)
C
C ***** Attention ***
C
DIMENSION VVR( 2000), PVR( 2000)
C
COMMON /PFAC/ ZFAC
COMMON /CNST/ LAYX, LAYY, MESH
C
C ***** DATA ( 2 ) *****************************************************
C
*** OPEN ( 11, FILE='E2_F3D_FEM_FIG' )
*** OPEN ( 11, FILE='E3_F3D_FEM_FIG' ) ! Re= 100 DT = 0.3
*** OPEN ( 11, FILE='E4_F3D_FEM_FIG' ) ! Re= 100 DT = 0.25
*** OPEN ( 11, FILE='E6_F3D_FEM_FIG' )
*** OPEN ( 11, FILE='E10_F3D_FEM_FIG.DAT' )
C
OPEN ( 11, FILE='UE3_F3D_FEM_FIG' ) ! Re = 100 DT = 0.3
C
*** OPEN ( 11, FILE='UE4_F3D_FEM_FIG' ) ! Re = 100 DT = 0.2
*** OPEN ( 11, FILE='UE10_F3D_FEM_FIG' ) ! Re = 100 DT = 0.05
*** OPEN ( 11, FILE='UE20_F3D_FEM_FIG' ) ! Re = 400 DT = 0.05
*** OPEN ( 11, FILE='UE30_F3D_FEM_FIG' ) ! Re = 1000 DT = 0.04
C
C **********************************************************************
C
CALL INSB ( U, V, W, P, X, Y, Z, PVR, VVR,
1 XLX, YLY, LXY, LYZ, LZX, LOP, NNP, LX, LY,
2 NXY, NYZ, NZX, NEZ, IKY )
C
CALL VGA@
C
IF ( IKY.EQ.1) CALL MSH3 ( XLX, YLY, LAYX, LAYY, NEX, MESH )
C
IF ( IKY.EQ.2) CALL PTVR ( VVR, LOP, 1, NEX )
IF ( IKY.EQ.2) CALL PTVR ( PVR, LOP, 2, NEX )
IF ( IKY.EQ.2) CALL VW3D ( X, Y, Z, NNP, U, V, W )
C
IF ( IKY.EQ.3) THEN
C
CALL PLTT ( X, Y, LXY, P, NNP, NXY, 3, LOP )
CALL PLTN ( X, Y, U, V, NNP, LXY, NXY, 3 )
C
CALL PLTT ( Y, Z, LYZ, P, NNP, NYZ, 4, LOP )
CALL PLTN ( Y, Z, V, W, NNP, LYZ, NYZ, 4 )
C
CALL PLTT ( Z, X, LZX, P, NNP, NZX, 5, LOP )
CALL PLTN ( Z, X, W, U, NNP, LZX, NZX, 5 )
C
END IF
C -----------------------------------------------------
C
IF ( IKY.EQ.4) CALL CMPLT ( U, V, X, Y, NNP )
C
C -----------------------------------------------------
C
CLOSE ( 11)
C
STOP
END
C **********************************************************************
C
SUBROUTINE INSB ( U, V, W, P, X, Y, Z, PVR,
1 VVR, XLX, YLY, LXY, LYZ, LZX, LOP, NNP,
2 LX, LY, NXY, NYZ, NZX, NEZ, IKY )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION X ( NNP), Y ( NNP), Z ( NNP), U ( NNP),
1 V ( NNP), W ( NNP), P ( NNP), XLX( LX),
2 YLY( LY), LXY( NXY,4), LYZ( NYZ,4), LZX( NZX,4)
C
DIMENSION VVR( 2000), PVR( 2000)
C
C ----- Attention : Max.No_Meshes < = 200 --------
C
COMMON /PFAC/ ZFAC
COMMON /CNST/ LAYX, LAYY, MESH
C
C -----------------
C
CALL TTLE
C
C -----------------
C
WRITE(6,*) ' '
WRITE(6,*) ' *** IKY = 1: ( Meshes) = 2: ( Var.&Vel.& Perspective
1 View )'
C
C ----- Attention --- IKY = 4 : U-V Distr. at Z = 0.5 : NEZ > = 4 ---
C
IF ( NEZ.GE.4)
1WRITE(6,*) ' = 3: ( X-Y face)( Y-Z face)( Z-X face) = 4: U-V Dist
2r. at Z = 0.5 ***'
C
C ----- Attention : Not defined for NEX = 3 etc. ---
C
READ(6,1000) IKY
1000 FORMAT(I3)
WRITE(6,2000) IKY
2000 FORMAT(/,' * IKY =',I2 )
READ(11,1200) LAYX, LAYY, LAYZ, NEX, MESH
1200 FORMAT(5I3)
WRITE(6,2100) LAYX, LAYY, LAYZ, NEX, MESH
2100 FORMAT(/,' LAYX =', I3,' LAYY =', I3,' LAYZ =', I3,' NEX =',
1 I3,' MESH =',I2 )
C
READ(11,1300) ( XLX( I), I=1, LAYX ), ( YLY( I), I=1, LAYY )
1300 FORMAT(9D14.7)
READ(11,1400) RE, DT, LOP, ZFAC
1400 FORMAT(F7.1,F7.4,I5,D11.4)
READ(11,1300) ( VVR( I), I=1, LOP )
READ(11,1300) ( PVR( I), I=1, LOP )
READ(11,1500) NNP
1500 FORMAT( I5)
WRITE(6,2200) RE, DT, LOP, ZFAC, NNP
2200 FORMAT(/,' * Re =',F7.1,' DT =',F7.4,' Loop =', I5,
1 ' Zfac =',F6.3,' NNP =', I6 )
C
C ----------------------------
C
READ(11,1300) ( X( I), Y( I), Z( I), I = 1, NNP )
READ(11,1300) ( U( I), V( I), W( I), I = 1, NNP )
READ(11,1700) NXY, NYZ, NZX
1700 FORMAT(3I5)
READ(11,1300) ( P( I), I = 1, NNP )
READ(11,2300) (( LXY( I,J), I = 1, NXY ), J = 1, 4 )
READ(11,2300) (( LYZ( I,J), I = 1, NYZ ), J = 1, 4 )
READ(11,2300) (( LZX( I,J), I = 1, NZX ), J = 1, 4 )
2300 FORMAT(16I5)
C
C ----- To set Node numbers for U-V Distribution on Z = 0.5 ---
C
CALL STND ( X, Y, LXY, NXY, NNP )
C
C -------------------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE STND ( X, Y, LXY, NXY, NNP )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION X( NNP), Y( NNP), LXY( NXY,4)
C
COMMON /CMPD/ KSY( 200), NY, KSX( 200), NX
C
C ----- To set Node numbers for U-V Distribution on Z = 0.5 ---
C ********
C ==> LXY( NXY,4) <==> NDE( NEL,4) ---
C
C ----- All points information ---------
C
NY = 0
C
DO 100 I = 1, NXY
C
NI = LXY( I,1) ! or 2
C ----- Attention ----------------------
C
IF ( DABS( X( NI)-0.5D0).LE.0.01D0) THEN
NY = NY + 1 ! NY < = 200
KSY( NY) = NI
END IF
C
100 CONTINUE
C --------------------------------------
DO 200 I = 1, NXY
C
NI = LXY( I,2) ! or 2
C ----------------------------
C
IF ( DABS( X( NI)-0.5D0).LE.0.01D0) THEN
NY = NY + 1 ! NY < = 200
KSY( NY) = NI
END IF
C
200 CONTINUE
C --------------------------------------
DO 300 I = 1, NXY
C
NI = LXY( I,3) ! or 2
C ----------------------------
C
IF ( DABS( X( NI)-0.5D0).LE.0.01D0) THEN
NY = NY + 1 ! NY < = 200
KSY( NY) = NI
END IF
C
300 CONTINUE
C --------------------------------------
DO 400 I = 1, NXY
C
NI = LXY( I,4) ! or 2
C ----------------------------
C
IF ( DABS( X( NI)-0.5D0).LE.0.01D0) THEN
NY = NY + 1 ! NY < = 200
KSY( NY) = NI
END IF
C
400 CONTINUE
C --------------------------------------
NX = 0
C
DO 500 I = 1, NXY
C
NI = LXY( I,1)
IF ( DABS( Y( NI)-0.5D0).LE.0.01D0) THEN
NX = NX + 1 ! NX < = 200
KSX( NX) = NI
END IF
C
500 CONTINUE
C ---------------------------------------
DO 600 I = 1, NXY
C
NI = LXY( I,2)
IF ( DABS( Y( NI)-0.5D0).LE.0.01D0) THEN
NX = NX + 1 ! NX < = 200
KSX( NX) = NI
END IF
C
600 CONTINUE
C --------------------------------------
DO 700 I = 1, NXY
C
NI = LXY( I,3)
IF ( DABS( Y( NI)-0.5D0).LE.0.01D0) THEN
NX = NX + 1 ! NX < = 200
KSX( NX) = NI
END IF
C
700 CONTINUE
C --------------------------------------
C
DO 800 I = 1, NXY
C
NI = LXY( I,4)
IF ( DABS( Y( NI)-0.5D0).LE.0.01D0) THEN
NX = NX + 1 ! NX < = 200
KSX( NX) = NI
END IF
C
800 CONTINUE
C --------------------------------------
C
WRITE(6,2000) NY, NX
2000 FORMAT(/,' ** NY =',I4,' NX =',I4)
C
* WRITE(6,*) ' * KSY=', ( KSY(I), I=1, NY )
* WRITE(6,*) ' * KSX=', ( KSX(I), I=1, NX )
C
RETURN
END
C **********************************************************************
C
SUBROUTINE TTLE
C
IMPLICIT REAL*8 (A-H,O-Z)
C
WRITE(6,*) ' '
WRITE(6,*) 'C ************************************************* C'
WRITE(6,*) 'C C'
WRITE(6,*) 'C FIG - F3D - FEM - R C'
WRITE(6,*) 'C C'
WRITE(6,*) 'C ( Rectangular Parallelepiped Elements ) C'
WRITE(6,*) 'C C'
WRITE(6,*) 'C ( Version : 3.8 ) C'
WRITE(6,*) 'C C'
WRITE(6,*) 'C ( Graphics Output ) C'
WRITE(6,*) 'C C'
WRITE(6,*) 'C 3D Fluid Analysis Results C'
WRITE(6,*) 'C C'
WRITE(6,*) 'C Copyright 2011 : Yasuhiro MATSUDA C'
WRITE(6,*) 'C C'
WRITE(6,*) 'C ************************************************* C'
C
RETURN
END
C **********************************************************************
C
SUBROUTINE CMPLT ( UN, VN, X, Y, NNP )
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 ICR, IX1, IY1, IX2, IY2
C
DIMENSION UN( NNP), VN( NNP), X( NNP), Y( NNP)
C
COMMON /TIMC/ DT, LPLM, LOP1, LOP, ALP, RE
C
NPL = 0
C ---------------
C
IX1 = 0
IY1 = 0
IX2 = 639
IY2 = 479
ICR = 7
C
CALL CLEAR_SCREEN_AREA@( IX1, IY1, IX2, IY2, ICR )
C
CALL RFCMP ( X, Y, UN, VN, NNP )
C
RETURN
END
C **********************************************************************
C
SUBROUTINE RFCMP ( X, Y, U, V, NNP )
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION X( NNP), Y( NNP), U( NNP), V( NNP)
C
COMMON /CMPD/ KSY( 200), NY, KSX( 200), NX
C
INTEGER*2 X0, Y0, X1, Y1, X2, Y2, IX1, IY1, IX2, IY2, ICR, IC2
C
CHARACTER*80 STR0
C
DATA STR0 /'** Results for F3D-FEM-R ** ( Velocity U-V Distributi
1on Z = 0.5 )'/
C
NPL = 0
C
ICR = 11
IC2 = 16
C
CALL SET_TEXT_ATTRIBUTE@( 3, 1., 0., 0. )
C
CALL DRAW_TEXT@( STR0, 10, 5, IC2 )
C
X1 = 345
Y1 = 50
X0 = 70
Y0 = 50
C
IWD = 400
C
X2 = IDNINT( DFLOAT( IWD)* 0.5D0)
Y2 = IDNINT( DFLOAT( IWD)* 0.5D0)
C
IX1 = X0
IY1 = Y0
IX2 = IX1 + IWD
IY2 = IY1 + IWD
C
CALL RECTANGLE@( IX1, IY1, IX2, IY2, IC2 )
C
IX1 = X0 + X2
IY1 = Y0 + IWD
C
IX2 = IX1
IY2 = Y0
C
CALL RECTANGLE@( IX1, IY1, IX2, IY2, IC2 )
C
IX1 = X0
IY1 = Y0 + IDNINT( DFLOAT( IWD)* 0.5D0)
C
IX2 = IX1 + IWD
IY2 = IY1
C
CALL DRAW_LINE@( IX1, IY1, IX2, IY2, IC2 )
C
IX1 = 40
IY1 = 80 + NPL* 20
C
IX2 = 80
IY2 = IY1
C
C -------------- Scaling -----------------------------
C
CALL SET_TEXT_ATTRIBUTE@( 3, 1.5, 0., 0. )
C
CALL DRAW_TEXT@(' 1', 450, 475, IC2 )
CALL DRAW_TEXT@('-1', 50, 475, IC2 )
CALL DRAW_TEXT@(' 0', 250, 475, IC2 )
C
CALL DRAW_TEXT@(' 1', 40, 70, IC2 )
CALL DRAW_TEXT@(' 0', 40, 270, IC2 )
C
C ---------------------------------------------------
C
** WRITE(6,*) ' * KSY =', KSY
C
C ----- Computed U ----------------------------------
C
DO 300 J = 1, NY
C
NN = KSY( J)
C
*** WRITE(6,*) ' * U =', U( NN)
C
IX1 = X0 + X2 + IDNINT( U( NN)* DFLOAT( X2))
IY1 = Y0 + IWD - IDNINT( Y( NN)* DFLOAT( IWD))
C
IX2 = 2
IY2 = 2
C
CALL ELLIPSE@( IX1, IY1, IX2, IY2, ICR )
C
NN1 = KSY( J+1)
C
IX2 = X0 + X2 + IDNINT( U( NN1)* DFLOAT( X2))
IY2 = Y0 + IWD - IDNINT( Y( NN1)* DFLOAT( IWD))
C
300 CONTINUE
C
** WRITE(6,*) ' * NX, NY=', NX, NY
C
C ----- Computed V ----------------------------------
C
DO 500 J = 1, NX
C
NN = KSX( J)
C
IX1 = X0 + IDNINT( X( NN)* DFLOAT( IWD))
IY1 = Y0 + IDNINT( DFLOAT( IWD)* 0.5 )
1 - IDNINT( V( NN)* DFLOAT( Y2))
C
IX2 = 2
IY2 = 2
C
CALL FILL_ELLIPSE@( IX1, IY1, IX2, IY2, ICR )
C
NN1 = KSX( J+1)
C
IX2 = X0 + IDNINT( X( NN1)* DFLOAT( IWD))
IY2 = Y0 + IDNINT( DFLOAT( IWD)* 0.5 )
1 - IDNINT( V( NN1)* DFLOAT( Y2))
C
IF ( J.EQ.NX ) THEN
IX2 = X0 + IWD
IY2 = Y0 + IDNINT( DFLOAT( IWD)* 0.5 )
END IF
C
500 CONTINUE
C -----------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE PLTN ( X, Y, UN, VN, NNP, NDE, NEL, IKY )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION X( NNP), Y( NNP), UN( NNP), VN( NNP), NDE( NEL,4)
C
INTEGER*2 ICR, IX1, IY1, IX2, IY2
C
COMMON /PFAC/ ZFAC
C
ICR = 15
IWD = 160
C
IX1 = 20
IY1 = 20
C
C ----------------------------
IF ( IKY.EQ.3) THEN
IX2 = IX1 + IWD
IY2 = IY1 + IWD
END IF
C ----------------------------
IF ( IKY.EQ.4) THEN
IX1 = IX1 + 200
IX2 = IX1 + IWD
IY2 = IY1 + IDNINT( DFLOAT( IWD)* ZFAC )
END IF
C ----------------------------
IF ( IKY.EQ.5) THEN
IX1 = IX1 + 400
IX2 = IX1 + IDNINT( DFLOAT( IWD)* ZFAC )
IY2 = IY1 + IWD
END IF
C -------------------------------------------------
C
CALL RECTANGLE@( IX1, IY1, IX2, IY2, ICR )
C
C -------------------------------------------------
C
ICR = 11
C
X0 = 0.1D0
Y0 = 0.1D0
C
SCX = 0.8D0
SCY = 0.8D0
C
SCL = 0.2D0
IF ( NEL.LE.75) SCL = 0.4D0 ! = < 4 x 4 x 4
C
C ----------------------------------------------------------------------
DO 100 JE = 1, NEL
C
ND = NDE( JE,1)
C
XS = X0 + X( ND)* SCX
YS = Y0 + Y( ND)* SCY
C
AL = DSQRT( UN( ND)** 2 + VN( ND)** 2 )
C
C -------------------------------------------------------
C
IF ( AL.GE.0.0001D0) THEN
C
C --------------------------------
C
XE = XS + UN( ND)* SCL
C
IF ( IKY.EQ.5) XE = XS - UN( ND)* SCL ! Direction of 'Z'
C ----- Attention -------
C
YE = YS + VN( ND)* SCL
C
C --------------------------------
C
IF ( IKY.EQ.3) THEN
IX1 = IDNINT( 200.D0* XS )
IY1 = IDNINT( 200.D0* ( 1.D0 - YS ))
IX2 = IDNINT( 200.D0* XE )
IY2 = IY1 -IDNINT( 200.D0* VN( ND)* SCL ) ! X-Y
END IF
C
IF ( IKY.EQ.4) THEN
IX1 = IDNINT( 200.D0* XS ) + 200
IY1 = IDNINT( 200.D0* ( 1.D0 - YS ))
IX2 = IDNINT( 200.D0* XE ) + 200
IY2 = IY1 + IDNINT( 200.D0* VN( ND)* SCL ) ! Y-Z
C
C ---------------- Attention: Direction of "Z" - W ----
END IF
C
IF ( IKY.EQ.5) THEN
IX1 = IDNINT( 200.D0* XS) + 400
IY1 = IDNINT( 200.D0* ( 1.D0 - YS))
IX2 = IDNINT( 200.D0* XE ) + 400
IY2 = IY1 + IDNINT( 200.D0* VN( ND)* SCL) ! Z-X
C
C ---------------- Attention: Direction of "X" - U ----
END IF
C
C -------------------------------------------------
C
CALL DRAW_LINE@( IX1, IY1, IX2, IY2, ICR )
C
C -------------------------------------------------
C
END IF
C ------------------------------------------------------------
C
100 CONTINUE
C ----------------------------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE PLTT ( X, Y, NDE, TEP, NNP, NEL, IKY, LOP )
C
C ----- Plot contour lines for rectangular meshes ---
C
C ( Linear interpolation )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION X ( NNP), Y ( NNP), NDE( NEL,4), TEP( NNP)
C
DIMENSION PX( 4), PY( 4), HVN( 501)
C
INTEGER*2 ICR, IX1, IY1, IX2, IY2, IX3, IY3
C
COMMON /PFAC/ ZFAC
C
CHARACTER*11 TSB( 3)
DATA TSB / '( Z = 0.5 )', '( X = 0.5 )', '( Y = 0.5 )'/
C
C ***** Attention **********
C
IF ( IKY.EQ.5 ) THEN
C
C ----------------------------
DO 100 I = 1, NNP
C
Y( I) = 1.D0 - Y( I)
C
100 CONTINUE
C ----------------------------
C
END IF
C
ICR = 15
C
X0 = 0.1D0
Y0 = 0.1D0
C
SCX = 0.8D0
SCY = 0.8D0
C
RMI = 100.D0
RMX = 0.D0
C
C ------------------------------------------
DO 200 I = 1, NNP
C
IF ( TEP( I).LE.RMI ) RMI = TEP( I)
IF ( TEP( I).GE.RMX ) RMX = TEP( I)
C
200 CONTINUE
C ------------------------------------------
C
WRITE(6,1231) IKY, RMI, RMX
1231 FORMAT(/,' ** IKY =',I3,' Rmin.=',F8.3,5X,' Rmax.=',F8.3 )
C
IF ( LOP.LT.500) DRM = 0.005D0
IF ( LOP.LT.500) NM = 400
C
IF ( LOP.GE.500) DRM = 0.01D0
IF ( LOP.GE.500) NM = 200
C
IF ( LOP.GE.1100) DRM = 0.005D0
IF ( LOP.GE.1100) NM = 200
C
C ------------------------------------------------
IF ( NNP.LE.125) THEN ! 4 X 4 X 4
NM = 25
DRM = ( RMX - RMI)/ DFLOAT( NM)
END IF
C ------------------------------------------------
IF ( NNP.GE.126 .AND. NNP.LE.300) THEN
NM = 400
DRM = ( RMX - RMI)/ DFLOAT( NM)
END IF
C ------------------------------------------------
WRITE(6,1133) DRM
1133 FORMAT(' ** dP =',F7.3)
C
C ----------------------------------------
DO 300 I = 1, NM+1
C
HVN( I) = RMX - DRM* DFLOAT( I-1)
C
300 CONTINUE
C ----------------------------------------
C
*** WRITE(6,*) ' * HVN=', ( HVN( I), I=1, NM+1 )
C
IWD = 160
IX1 = 20
C
IF ( IKY.EQ.4) IX1 = IX1 + 200
IF ( IKY.EQ.5) IX1 = IX1 + 400
IY1 = 220
IF ( IKY.EQ.4) IY1 = 220 + IDNINT( DFLOAT( IWD)*( 1.D0 - ZFAC))
C
IX3 = IX1 + IWD
IY3 = IY1 + IWD
IF ( IKY.EQ.4) IY3 = IY1 + IDNINT( DFLOAT( IWD)* ZFAC )
IF ( IKY.EQ.5) IX3 = IX1 + IDNINT( DFLOAT( IWD)* ZFAC )
C
CALL RECTANGLE@( IX1, IY1, IX3, IY3, ICR )
C
C----- Plot contour lines ---
C
ICR = 14
C
IWD = 200
PXV = 0.D0
PYV = 0.D0
C
C ---------------------------------------------------------
DO 400 L = 1, NEL
C
C ------------------------------------------
DO 420 IH = 1, NM + 1
C
HV = HVN( IH)
ICL = 0
C
C ------------------------------
DO 440 J = 1, 4
C
I1 = NDE ( L,J)
IF ( J .NE. 4) THEN
I2 = NDE ( L,J+1)
ELSE
I2 = NDE( L,1)
END IF
C
IF ( DABS( TEP( I1)-TEP( I2)).LE.0.000001D0) THEN
IF ( DABS( TEP( I1)-HV ).LE.0.000001D0) THEN
ICL = ICL + 1
PX( ICL) = X( I1)* SCX
PY( ICL) = Y( I1)* SCY
END IF
C
ELSE IF (( ( TEP( I1).LE. HV ) .AND. ( HV.LT.TEP( I2)))
1 .OR.(( TEP( I1).GE. HV ) .AND. ( HV.GT.TEP( I2))))
2 THEN
C
ICL = ICL + 1
RT = ( HV - TEP( I1))/( TEP( I2) - TEP( I1))
C
PX( ICL) = ( X( I1) + ( X( I2) - X( I1))* RT )* SCX
PY( ICL) = ( Y( I1) + ( Y( I2) - Y( I1))* RT )* SCY
C
END IF
C
440 CONTINUE
C ---------------------------------------
C
IF ( ICL.GE.2) THEN
C
DO 460 J = 1, ICL
C
JM = J
D2 = ( PXV - PX( J))** 2 + ( PYV - PY( J))** 2
C
C ---------------------------------
DO 480 JJ = J+1, ICL
C
D2V = ( PXV - PX( JJ))** 2 + ( PYV - PY( JJ))** 2
IF ( D2V.LT.D2 ) THEN
JM = JJ
D2 = D2V
END IF
C
480 CONTINUE
C ---------------------------------
C
PXV = PX( JM)
PYV = PY( JM )
PX( JM) = PX( J)
PY( JM) = PY( J)
PX( J) = PXV
PY( J) = PYV
C
460 CONTINUE
C
PXV = PX( ICL)
PYV = PY( ICL)
C
END IF
C
IF ( ICL.GT.0) THEN
C
C ---------------------------------------------
DO 490 J = 1, ICL-1
C
IX1 = IDNINT( 200.D0* ( X0 + PX( J)))
IF ( IKY.EQ.4) IX1 = IX1 + 200
IF ( IKY.EQ.5) IX1 = IX1 + 400
C
IY1 = 401 - IDNINT( 200.D0* ( Y0 + PY( J)))
IX2 = IDNINT( 200.D0* ( X0 + PX( J+1)))
IF ( IKY.EQ.4) IX2 = IX2 + 200
IF ( IKY.EQ.5) IX2 = IX2 + 400
IY2 = 401 - IDNINT( 200.D0* ( Y0 + PY( J+1)))
C
CALL DRAW_LINE@( IX1, IY1, IX2, IY2, ICR )
C
490 CONTINUE
C ------------------------------------------------
DO 495 J = 2, ICL
C
IX1 = IDNINT( 200.D0* ( X0 + PX( 1)))
IF ( IKY.EQ.4) IX1 = IX1 + 200
IF ( IKY.EQ.5) IX1 = IX1 + 400
IY1 = 401 - IDNINT( 200.D0* ( Y0 + PY( 1)))
IX2 = IDNINT( 200.D0* ( X0 + PX( J)))
C
IF ( IKY.EQ.4) IX2 = IX2 + 200
IF ( IKY.EQ.5) IX2 = IX2 + 400
IY2 = 401 - IDNINT( 200.D0* ( Y0 + PY( J)))
C
495 CONTINUE
C -----------------------------------------
C
END IF
C
420 CONTINUE
C ----------------------------------------------
C
400 CONTINUE
C ---------------------------------------------------------
C
ICR = 15
C
IY1 = 240 + IWD
IX1 = 10
IX2 = IX1 + 40
IY2 = IY1 - 40
C
CALL DRAW_LINE@( IX1, IY1, IX2, IY1, ICR )
CALL DRAW_LINE@( IX1, IY1, IX1, IY2, ICR )

CALL SET_TEXT_ATTRIBUTE@( 3, 1.3, 0., 0. )
C
CALL DRAW_TEXT@( 'X', 60, 450, ICR )
CALL DRAW_TEXT@( 'Y', 20, 410, ICR )
C
IX1 = 10 + IWD
IY1 = IY1 - 40
IX2 = IX1 + 40
IY2 = IY1 + 40
C
CALL DRAW_LINE@( IX1, IY1, IX2, IY1, ICR )
CALL DRAW_LINE@( IX1, IY1, IX1, IY2, ICR )
CALL DRAW_TEXT@( 'Z', IX1+10, 450, ICR )
CALL DRAW_TEXT@( 'Y', IX1+50, 410, ICR )
C
IX1 = 10 + IWD + IWD + 140
IX2 = IX1 + 40
IY2 = IY1 + 40
C
CALL DRAW_LINE@( IX1, IY1, IX2, IY1, ICR )
CALL DRAW_LINE@( IX2, IY1, IX2, IY2, ICR )
CALL DRAW_TEXT@( 'Z', IX1-20, 410, ICR )
CALL DRAW_TEXT@( 'X', IX1+20, 450, ICR )
C
C ----------------------------------------------------------
C
CALL SET_TEXT_ATTRIBUTE@( 3, 1.3, 0., 0. )
C
CALL DRAW_TEXT@( TSB( 1), 50, 210, ICR )
CALL DRAW_TEXT@( TSB( 2), 50+ IWD, 210, ICR )
CALL DRAW_TEXT@( TSB( 3), 50+2* IWD, 210, ICR )
C
C ----------------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE PTVR ( XX, LOP, KEY, NEX )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
INTEGER*2 ICR, IX0, IY0, IX1, IY1
C
C ----- Attention N.LE.2000 ----------------
C
DIMENSION XX( 2000)
C
INTEGER*2 IVR, IXX
DIMENSION IVR( 2000), IXX( 2000)
C
CHARACTER*70 STR0
CHARACTER*4 STR1( 3)
C
DATA STR0
1/'** Results of F3D-FEM-R ( Max_rel_var_Velocity & Pressure (%))
2**'/
DATA STR1/'0 ', 'LP/2', 'LOP'/
C
** WRITE(6,*) ' ** XX( LOP)=', XX( LOP)
C
IX0 = 50
IY0 = 40
C
ICR = 15
C
CALL SET_TEXT_ATTRIBUTE@( 3, 1., 0., 0. )
C
CALL DRAW_TEXT@( STR0, 90, 10, ICR )
C
CALL DRAW_TEXT@( STR1( 1), 45, 450, ICR )
CALL DRAW_TEXT@( STR1( 2), 315, 450, ICR )
CALL DRAW_TEXT@( STR1( 3), 605, 450, ICR )
C
IX1 = IX0 + 580
IY1 = IY0 + 400
C
CALL RECTANGLE@( IX0, IY0, IX1, IY1, 15 )
C
CALL DRAW_TEXT@( '1000', 15, 30, ICR )
CALL DRAW_TEXT@( ' 100', 15, 110, ICR )
CALL DRAW_TEXT@( ' 10', 15, 190, ICR )
C
CALL DRAW_TEXT@( ' 1', 15, 270, ICR )
CALL DRAW_TEXT@( ' 0.1', 15, 350, ICR )
CALL DRAW_TEXT@( '0.01', 15, 430, ICR )
C
CALL DRAW_TEXT@( '(%)', 15, 245, ICR )
C
C --------------------------------------
DO 100 I = 1, LOP
C
IF ( XX( I).LE.0.D0) GO TO 100
XX( I) = DLOG10( XX( I))
C
100 CONTINUE
C --------------------------------------
C
SCLN = 80.D0
C
IF ( NEX.GE.3. AND. NEX.LE.6) SCL = 20.D0
C
C --------------------------------------
DO 200 I = 1, LOP
C
IXX( I) = IX0 + IDNINT( DFLOAT( I)* 580.D0 / DFLOAT( LOP))
IVR( I) = 280 - IDNINT( XX( I)* SCLN )
C
200 CONTINUE
C --------------------------------------
C
** WRITE(6,*) ' ** IVR( LOP)=', IVR( LOP)
C
IF ( KEY.EQ.1) ICR = 11
IF ( KEY.EQ.2) ICR = 14
C
C -------------------------------------------
C
CALL POLYLINE@( IXX, IVR, LOP, ICR )
C
C -------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE MSH3 ( XLX, YLY, LAYX, LAYY, NEX, MESH )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION X ( 40), Y ( 40), GX ( 40), GY ( 40),
1 GXR ( 40), GYR ( 40), GZX ( 40), GZY ( 40),
2 GZXR( 40), GZYR( 40), GXRA( 40), GXRB( 40),
3 GYRA( 40), GYRB( 40), XLX ( LAYX), YLY ( LAYY)
C
COMMON / VIEW / HR, VR, DIS, DS2, OX, OY, OZ, SCA, RR
C
INTEGER*2 ICR, IX1, IY1, IX2, IY2, IX3, IY3, IX4, IY4,
1 IY5, IY6, IY7, IWK1, IWK2, IWK3, IWK4
C
C
CALL SET_TEXT_ATTRIBUTE@( 3, 1.3, 0., 0. )
C
ICR = 15
C
CALL DRAW_TEXT@(' F3D-FEM ', 60, 430, ICR )
C
IF ( NEX.EQ. 3) CALL DRAW_TEXT@( ' ( 3 X 3 X 3 )', 140, 430,ICR)
IF ( NEX.EQ. 4) CALL DRAW_TEXT@( ' ( 4 X 4 X 4 )', 140, 430,ICR)
C
IF ( NEX.EQ.10) CALL DRAW_TEXT@( ' ( 10 X 10 X 10)',
1 140, 430,ICR )
IF ( NEX.EQ.20) CALL DRAW_TEXT@( ' ( 20 X 20 X 20)',
1 140, 430,ICR )
IF ( NEX.EQ.30) CALL DRAW_TEXT@( ' ( 30 X 30 X 30)',
1 140, 430,ICR )
C
IF ( NEX.LE.10) CALL DRAW_TEXT@( 'Re = 100', 320, 430, ICR )
IF ( NEX.EQ.20) CALL DRAW_TEXT@( 'Re = 400', 320, 430, ICR )
IF ( NEX.EQ.30) CALL DRAW_TEXT@( 'Re = 1000', 320, 430, ICR )
C
C -------------------------------------
DO 100 I = 1, NEX
C
X( I) = XLX( I)
Y( I) = YLY( I)
C
100 CONTINUE
C -------------------------------------
C
ICR = 15
IWD = 200
C
IY0 = 30
C
IX1 = 100
IY1 = 100 + IY0
C
IX2 = IX1 + 67
IY2 = 50 + IY0
C
IX3 = IX2 + IWD
IY3 = IY2
C
IX4 = IX1 + IWD
IY4 = IY1
C
IY5 = IY1 + IWD
IY6 = IY3 + IWD
IY7 = IY4 + IWD
C
CALL DRAW_LINE@( IX1, IY1, IX4, IY4, ICR )
CALL DRAW_LINE@( IX1, IY1, IX1, IY5, ICR )
C
CALL DRAW_LINE@( IX1, IY5, IX4, IY5, ICR )
CALL DRAW_LINE@( IX4, IY7, IX3, IY6, ICR )
C
CALL DRAW_LINE@( IX1, IY1, IX2, IY2, ICR )
CALL DRAW_LINE@( IX2, IY2, IX3, IY3, ICR )
C
CALL DRAW_LINE@( IX3, IY3, IX3, IY6, ICR )
CALL DRAW_LINE@( IX4, IY1, IX4, IY5, ICR )
C
CALL DRAW_LINE@( IX4, IY4, IX4, IY1, ICR )
CALL DRAW_LINE@( IX4, IY1, IX3, IY2, ICR )
C
IWK1 = IY0
C
CALL DRAW_LINE@( IX2, IY2, IX2, IWK1, ICR )
C
C ----- Arrows -------------
C
IWK1 = IX2 - 5
IWK2 = 20 + IY0
IWK3 = IY0
C
CALL DRAW_LINE@( IX2, IWK3, IWK1, IWK2, ICR )
C
CALL SET_TEXT_ATTRIBUTE@( 3, 1.5, 0., 0. )
C
CALL DRAW_TEXT@( ' y', IX2-50, IY0, ICR )
C
IWK3 = IY0
IWK4 = IX2+5
C
CALL DRAW_LINE@( IX2, IWK3, IWK4, IWK2, ICR )
C
IWK3 = IX1-40
IWK4 = IY5+30
C
CALL DRAW_LINE@( IX1, IY5, IWK3, IWK4, ICR )
C
IWK1 = IY1+ IWD+30
IWK2 = IX1-33
IWK4 = IX1+ IWD+15 + IY0
C
CALL DRAW_LINE@( IWK3, IWK1, IWK2, IWK4, ICR )
C
CALL DRAW_TEXT@( ' z', IWK3-40, IWK1+20, ICR )
C
IWK2 = IX1-23
IWK4 = IY1+ IWD+28
C
CALL DRAW_LINE@( IWK3, IWK1, IWK2, IWK4, ICR )
C
IWK1 = IY3+ IWD
IWK2 = IX1+ IWD+140
IWK3 = IY3+ IWD
C
CALL DRAW_LINE@( IX3, IWK1, IWK2, IWK3, ICR )
C
IWK1 = IX1+ IWD+140
IWK2 = IY3+ IWD
IWK3 = IX1+ IWD+120
IWK4 = IY3+ IWD+5
C
CALL DRAW_LINE@( IWK1, IWK2, IWK3, IWK4, ICR )
C
IWK4 = IY3+ IWD-5
C
CALL DRAW_LINE@( IWK1, IWK2, IWK3, IWK4, ICR )
C
CALL DRAW_TEXT@( ' x', IWK1-10, IWK2+15, ICR )
C
C --------------------------------------
DO 200 I = 1, NEX
C
L = I
GX( L) = 200.D0* X( L)
GY( L) = 200.D0* Y( L)
C
200 CONTINUE
C --------------------------------------
IWK1 = IY1 + IWD
C
DO 300 L = 2, NEX
C
GXR( L) = 100.D0 + GX( L)
IX1 = IDNINT( GXR( L))
C
CALL DRAW_LINE@( IX1, IWK1, IX1, IY1, ICR )
C
300 CONTINUE
C --------------------------------------
C
IWK1 = 100 + IY0
IWK2 = 50 + IY0
C
DO 400 I = 2, NEX
C
GXR( I)= 100.D0 + GX( I)
IX1 = IDNINT( GXR( I))
IX2 = 167 + GX( I)
C
CALL DRAW_LINE@( IX1, IWK1, IX2, IWK2, ICR )
C
400 CONTINUE
C --------------------------------------
C
IWK1 = 100
IWK2 = 100 + IWD
C
DO 500 I = 2, NEX
C
GYR( I) = 100.D0 + GY( I)
IY1 = IDNINT( 100.D0 + GY( I)) + IY0
C
CALL DRAW_LINE@( IWK1, IY1, IWK2, IY1, ICR )
C
500 CONTINUE
C --------------------------------------
C
IWK1 = 300
IWK2 = 367
C
DO 600 I = 2, NEX
C
GYR( I) = 100.D0 + GY( I)
IY1 = IDNINT( 100.D0 + GY( I)) + IY0
IY2 = 47 + GY( I) + IY0
C
CALL DRAW_LINE@( IWK1, IY1, IWK2, IY2, ICR )
C
600 CONTINUE
C --------------------------------------
DO 700 I = 2, NEX
C
GZX( I) = 0.8D0* GX( I)
GZY( I) = 0.6D0* GX( I)
C
700 CONTINUE
C --------------------------------------
C
** WRITE(6,*) ' * IY0=', IY0
C
IF ( MESH.EQ.1) THEN
C
C -------------------------------------
DO 800 L = 2, NEX
C
GZXR( L) = 367.D0 - GZX( L)* 0.42D0
GYRA( L) = 250.D0 + GZY( L)* 0.42D0
GYRB( L) = 50.D0 + GZY( L)* 0.42D0
C
IX1 = GZXR( L)
IY1 = GYRA( L) + IY0
IY2 = GYRB( L) + IY0
C
CALL DRAW_LINE@( IX1, IY1, IX1, IY2, ICR )
C
800 CONTINUE
C --------------------------------------
DO 900 l = 2, NEX
C
GXRA( L) = 167.D0 - GZX( L)* 0.42D0
GXRB( L) = 367.D0 - GZX( L)* 0.42D0
GZYR( L) = 50.D0 + GZY( L)* 0.42D0
C
IX1 = IDNINT( GXRA( L))
IY1 = IDNINT( GZYR( L)) + IY0
IX2 = IDNINT( GXRB( L))
C
CALL DRAW_LINE@( IX1, IY1, IX2, IY1, ICR )
C
900 CONTINUE
C --------------------------------------
C
END IF
C
C ************************************
IF ( MESH.EQ.2 ) THEN
C
DO 930 L = 2, NEX
C
GZXR( L)= 367.D0 - GZX( L)* 0.4D0
GYRA( L)= 250.D0 + Gzy( L)* 0.4D0
GYRB( L)= 50.D0 + Gzy( L)* 0.4D0
C
IX1 = GZXR ( L)
IY1 = IDNINT( GYRA( L)) + IY0
IY2 = IDNINT( GYRB( L)) + IY0
C
CALL DRAW_LINE@( IX1, IY1, IX1, IY2, ICR )
C
930 CONTINUE
C --------------------------------------
C
DO 950 l = 2, NEX
C
GXRA( L) = 167.D0 - GZX( L)* 0.4D0
GXRB( L) = 367.D0 - GZX( L)* 0.4D0
GZYR( L) = 50.D0 + GZY( L)* 0.4D0
C
IX1 = IDNINT( GXRA( L))
IY1 = IDNINT( GZYR( L)) + IY0
IX2 = IDNINT( GXRB( L))
C
CALL DRAW_LINE@( IX1, IY1, IX2, IY1, ICR )
C
950 CONTINUE
C --------------------------------------
C
END IF
C **************************************
C
RETURN
END
C **********************************************************************
C
SUBROUTINE VW3D ( X, Y, Z, NNP, U, V, W )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
C ***** Attention : XXX ** ( Integer ) ***
C
DIMENSION X( NNP), Y( NNP), Z( NNP), U( NNP), V( NNP), W( NNP)
C
COMMON /VIEW/ HR, VR, DIS, DS2, OX, OY, OZ, SCA, RR

INTEGER*2 ICR, IX1, IY1, IX2, IY2, IX3, IY3, IX4, IY4, IX5, IY5,
1 IX6, IY6, IX7, IY7, IX8, IY8, IXS, IYS, IXE, IYE
C
SCA = 0.5D0
DIS = 2.D0
DS2 = 9.D0
C
VR = 40.D0
HR = - 30.D0
RR = - 35.D0
C
ICR = 15
C
SCL = 0.25D0
C
IF ( NNP.LE. 1000) SCL = 0.3D0
IF ( NNP.GT. 9000) SCL = 0.1D0 ! 21 X 21 X 21 = 9261
IF ( NNP.GT.29000) SCL = 0.03D0 ! 31 X 31 X 31 = 29791
C
ASCL = 0.01D0
C
XIN = X( 1)
XMX = X( 1)
C
YIN = Y( 1)
YMX = Y( 1)
C
ZIN = Z( 1)
ZMX = Z( 1)
C
C ---------------------------------
DO 100 I = 1, NNP
C
XIN = DMIN1( XIN, X( I))
XMX = DMAX1( XMX, X( I))
C
YIN = DMIN1( YIN, Y( I))
YMX = DMAX1( YMX, Y( I))
C
ZIN = DMIN1( ZIN, Z( I))
ZMX = DMAX1( ZMX, Z( I))
C
100 CONTINUE
C ---------------------------------
C
OX = ( XIN + XMX )/ 2.D0
OY = ( YIN + YMX )/ 2.D0
OZ = ( ZIN + ZMX )/ 2.D0
C
CALL VPNT ( XIN, YIN, ZIN, X1, Y1 )
CALL VPNT ( XMX, YIN, ZIN, X2, Y2 )
CALL VPNT ( XMX, YMX, ZIN, X3, Y3 )
CALL VPNT ( XIN, YMX, ZIN, X4, Y4 )
C
CALL VPNT ( XIN, YIN, ZMX, X5, Y5 )
CALL VPNT ( XMX, YIN, ZMX, X6, Y6 )
CALL VPNT ( XMX, YMX, ZMX, X7, Y7 )
CALL VPNT ( XIN, YMX, ZMX, X8, Y8 )
C
IX1 = X1
IY1 = 399-Y1
IX2 = X2
IY2 = 399-Y2
C
IX3 = X3
IY3 = 399-Y3
IX4 = X4
IY4 = 399-Y4
C
IX5 = X5
IY5 = 399-Y5
IX6 = X6
IY6 = 399-Y6
C
IX7 = X7
IY7 = 399-Y7
IX8 = X8
IY8 = 399-Y8
C
IDIF = 150
C
IX1 = IX1 + IDIF
IX2 = IX2 + IDIF
IX3 = IX3 + IDIF
IX4 = IX4 + IDIF
C
IX5 = IX5 + IDIF
IX6 = IX6 + IDIF
IX7 = IX7 + IDIF
IX8 = IX8 + IDIF
C
CALL DRAW_LINE@( IX1, IY1, IX2, IY2, ICR )
CALL DRAW_LINE@( IX2, IY2, IX3, IY3, ICR )
CALL DRAW_LINE@( IX3, IY3, IX4, IY4, ICR )
CALL DRAW_LINE@( IX4, IY4, IX1, IY1, ICR )
C
CALL DRAW_LINE@( IX1, IY1, IX5, IY5, ICR )
CALL DRAW_LINE@( IX5, IY5, IX6, IY6, ICR )
CALL DRAW_LINE@( IX6, IY6, IX7, IY7, ICR )
C
CALL DRAW_LINE@( IX7, IY7, IX8, IY8, ICR )
CALL DRAW_LINE@( IX8, IY8, IX5, IY5, ICR )
CALL DRAW_LINE@( IX8, IY8, IX4, IY4, ICR )
CALL DRAW_LINE@( IX4, IY4, IX3, IY3, ICR )
C
CALL DRAW_LINE@( IX3, IY3, IX7, IY7, ICR )
CALL DRAW_LINE@( IX7, IY7, IX6, IY6, ICR )
CALL DRAW_LINE@( IX6, IY6, IX2, IY2, ICR )
C
C ----- To plot X-Y-Z-> ---
C
XP = XIN - 0.1D0
YP = YIN - 0.1D0
ZP = ZMX
C
XC = XP
YC = YP
ZC = ZP
C
CALL VPNT ( XP, YP, ZP, X1, Y1 )
CALL VPNT ( XP+0.2D0, YP, ZP, X2, Y2 )
CALL VPNT ( XP, YP+0.2D0, ZP, X3, Y3 )
C
CALL VPNT ( XP, YP, ZP+0.2D0, X4, Y4 )
CALL VPNT ( XP, YP, ZP-0.1D0, X5, Y5 )
C
IX1 = X1
IY1 = 399-Y1
IX2 = X2
IY2 = 399-Y2
C
IX3 = X3
IY3 = 399-Y3
IX4 = X4
IY4 = 399-Y4
C
IX1 = IX1+ IDIF
IX2 = IX2+ IDIF
IX3 = IX3+ IDIF
IX4 = IX4+ IDIF
C
CALL DRAW_LINE@( IX1, IY1, IX2, IY2, ICR )
CALL DRAW_LINE@( IX1, IY1, IX3, IY3, ICR )
CALL DRAW_LINE@( IX1, IY1, IX4, IY4, ICR )
C
C ----- To plot 'X' ---
C
XP = XC + 0.25D0
YP = YC
ZP = ZC
C
CALL VPNT ( XP-0.02D0, YP+0.02D0, ZP, X1, Y1 )
CALL VPNT ( XP+0.02D0, YP-0.02D0, ZP, X2, Y2 )
CALL VPNT ( XP-0.02D0, YP-0.02D0, ZP, X3, Y3 )
CALL VPNT ( XP+0.02D0, YP+0.02D0, ZP, X4, Y4 )
C
IX1 = X1
IY1 = 399 - Y1
IX2 = X2
IY2 = 399 - Y2
IX3 = X3
IY3 = 399 - Y3
IX4 = X4
IY4 = 399 - Y4
C
IX1 = IX1 + IDIF
IX2 = IX2 + IDIF
IX3 = IX3 + IDIF
IX4 = IX4 + IDIF
C
CALL DRAW_LINE@( IX1, IY1, IX2, IY2, ICR )
CALL DRAW_LINE@( IX3, IY3, IX4, IY4, ICR )
C
C ----- To plot 'Y' ---
C
XP = XC
YP = YC + 0.25D0
ZP = ZC
C
CALL VPNT ( XP, YP, ZP, X1, Y1 )
CALL VPNT ( XP-0.02D0, YP+0.03D0, ZP, X2, Y2 )
C
CALL VPNT ( XP+0.02D0, YP+0.03D0, ZP, X3, Y3 )
CALL VPNT ( XP, YP-0.02D0, ZP, X4, Y4 )
C
IX1 = X1
IY1 = 399 - Y1
IX2 = X2
IY2 = 399 - Y2
C
IX3 = X3
IY3 = 399 - Y3
IX4 = X4
IY4 = 399 - Y4
C
IX1 = IX1 + IDIF
IX2 = IX2 + IDIF
IX3 = IX3 + IDIF
IX4 = IX4 + IDIF
C
CALL DRAW_LINE@( IX2, IY2, IX1, IY1, ICR )
CALL DRAW_LINE@( IX1, IY1, IX3, IY3, ICR )
CALL DRAW_LINE@( IX3, IY3, IX1, IY1, ICR )
CALL DRAW_LINE@( IX1, IY1, IX4, IY4, ICR )
C
C ----- To plot 'Z' ---
C
XP = XC
YP = YC
ZP = ZC + 0.25D0
C
CALL VPNT ( XP-0.02D0, YP+0.02D0, ZP, X1, Y1 )
CALL VPNT ( XP+0.02D0, YP+0.02D0, ZP, X2, Y2 )
C
CALL VPNT ( XP-0.02D0, YP-0.02D0, ZP, X3, Y3 )
CALL VPNT ( XP+0.02D0, YP-0.02D0, ZP, X4, Y4 )
C
IX1 = X1
IY1 = 399 - Y1
IX2 = X2
IY2 = 399 - Y2
C
IX3 = X3
IY3 = 399 - Y3
IX4 = X4
IY4 = 399 - Y4
C
IX1 = IX1 + IDIF
IX2 = IX2 + IDIF
IX3 = IX3 + IDIF
IX4 = IX4 + IDIF
C
CALL DRAW_LINE@( IX1, IY1, IX2, IY2, ICR )
CALL DRAW_LINE@( IX2, IY2, IX3, IY3, ICR )
CALL DRAW_LINE@( IX3, IY3, IX4, IY4, ICR )
C
S = DSIN( 0.2D0)
C = DCOS( 0.2D0)
EPS = 0.0001D0
C
ICR = 11
C ---------------------------------------------------------
DO 200 I = 1, NNP
C
AL = DSQRT( U( I)** 2 + V( I)** 2 + W( I)** 2 )
C ====================================================
C
IF ( AL.GT.EPS ) THEN
C
X1 = X( I)
Y1 = Y( I)
Z1 = Z( I)
C
CALL VPNT ( X1, Y1, Z1, XS, YS )

IXS = XS
IYS = 399 - YS
C
X1 = X1 + U( I)* SCL
Y1 = Y1 + V( I)* SCL
Z1 = Z1 + W( I)* SCL
C
CALL VPNT( X1, Y1, Z1, XE, YE )
C
IXE = XE
IYE = 399 - YE
C
IXS = IXS + IDIF
IXE = IXE + IDIF
C
CALL DRAW_LINE@( IXS, IYS, IXE, IYE, ICR )
C
AL = DSQRT(( XE-XS )** 2 + ( YE-YS )** 2 )
C ================================================
C
IF ( AL.GT.EPS) THEN
C
T = ASCL/ AL
C
IF ( T.GT.1.D0) T = 1.D0
C
X2 = XE + T* ( XS - XE )
Y2 = YE + T* ( YS - YE )
C
DX = X2 - XE
DY = Y2 - YE
C
X3 = XE + DX* C + DY* S
Y3 = YE - DX* S + DY* C
C
X4 = XE + DX* C - DY* S
Y4 = YE + DX* S + DY* C
C
IX3 = X3
IY3 = 399 - Y3
IX4 = X4
C
IY4 = 399 - Y4
IXE = XE
IYE = 399 - YE
C
IX4 = IX4 + IDIF
IXE = IXE + IDIF
IX3 = IX3 + IDIF
C
CALL DRAW_LINE@( IX4, IY4, IXE, IYE, ICR )
C
CALL DRAW_LINE@( IXE, IYE, IX3, IY3, ICR )
C
END IF
C
END IF
C
200 CONTINUE
C -----------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE VPNT ( X, Y, Z, VX, VY )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
COMMON /VIEW/ HR, VR, DIS, DS2, OX, OY, OZ, SCA, RR
C
XX = X - OX
YY = Y - OY
ZZ = Z - OZ
C
PI = 3.141592653589793 D0/ 180.D0
C
RH = PI* HR
RV = PI* VR
RL = PI* RR
C
X1 = XX* DCOS( RV) + YY* DSIN( RV)
Y2 = - XX* DSIN( RV) + YY* DCOS( RV)
X2 = X1* DCOS( RH) + ZZ* DSIN( RH)
C
Z1 = - X1* DSIN( RH) + ZZ* DCOS( RH)
X3 = X2* DCOS( RL) + Y2* DSIN( RL)
Y3 = - X2* DSIN( RL) + Y2* DCOS( RL)
C
D = DIS - Z1
C
VX = 450.D0* ( X3* DS2/( D + DS2)* SCA + 0.5D0) + 100.D0
VY = 450.D0* ( Y3* DS2/( D + DS2)* SCA + 0.5D0)
C
RETURN
END
************************************************************************