SUBROUTINE VELVCT (U,LU,V,LV,M,N,FLO,HI,NSET,LENGTH)
C
COMMON/VPLT/ IFLV,IVU1,IVU2,IVRT
C DECLARATIONS -
C
COMMON /VEC1/ ASH ,EXT ,ICTRFG ,ILAB ,
+ IOFFD ,IOFFM ,ISX ,ISY ,
+ RMN ,RMX ,SIDE ,SIZE ,
+ XLT ,YBT ,ZMN ,ZMX
C
COMMON /VEC2/ BIG,IX0,IX1,INCX,IY0,IY1,INCY
COMMON/LITARR/ LITARFL
COLD DATA LITARFL/1/
C
C FORCE THE BLOCK DATA ROUTINE, WHICH SETS DEFAULT VARIABLES, TO LOAD.
C
EXTERNAL VELDAT
C
C ARGUMENT DIMENSIONS.
C
DIMENSION U(LU,N) ,V(LV,N) ,SPV(2)
CHARACTER*14 LABEL
REAL WIND(4), VIEW(4), IAR(4)
COMMON/TOPOG/ CTP(2000),HTP(2000),Z0,ITOP
C
C ---------------------------------------------------------------------
C
C INTERNAL PARAMETERS OF VELVCT ARE AS FOLLOWS. THE DEFAULT VALUES OF
C THESE PARAMETERS ARE DECLARED IN THE BLOCK DATA ROUTINE VELDAT.
C
C NAME DEFAULT FUNCTION
C ---- ------- --------
C
C BIG R1MACH(2) CONSTANT USED TO INITIALIZE
C POSSIBLE SEARCH FOR HI.
C
C EXT 0.25 THE LENGTHS OF THE SIDES OF THE
C PLOT ARE PROPORTIONAL TO M AND
C N WHEN NSET IS LESS THAN OR
C EQUAL TO ZERO, EXCEPT WHEN
C MIN(M,N)/MAX(M,N) IS LESS THAN
C EXT, IN WHICH CASE A SQUARE
C GRAPH IS PLOTTED.
C
C ICTRFG 1 FLAG TO CONTROL THE POSITION OF
C THE ARROW RELATIVE TO A BASE
C POINT AT (MX,MY).
C
C ZERO - CENTER AT (MX,MY)
C
C POSITIVE - TAIL AT (MX,MY)
C
C NEGATIVE - HEAD AT (MX,MY)
C
C ILAB 0 FLAG TO CONTROL THE DRAWING OF
C LINE LABELS.
C
C ZERO - DO NOT DRAW THE LABELS
C
C NON-ZERO - DRAW THE LABELS
C
C INCX 1 X-COORDINATE STEP SIZE FOR LESS
C DENSE ARRAYS.
C
C INCY 1 Y-COORDINATE STEP SIZE.
C
C IOFFD 0 FLAG TO CONTROL NORMALIZATION
C OF LABEL NUMBERS.
C
C ZERO - INCLUDE A DECIMAL POINT
C WHEN POSSIBLE
C
C NON-ZERO - NORMALIZE ALL LABEL
C NUMBERS BY ASH
C
C IOFFM 0 FLAG TO CONTROL PLOTTING OF
C THE MESSAGE BELOW THE PLOT.
C
C ZERO - PLOT THE MESSAGE
C
C NON-ZERO - DO NOT PLOT IT
C
C RMN 160. ARROW SIZE BELOW WHICH THE
C HEAD NO LONGER SHRINKS, ON A
C 2**15 X 2**15 GRID.
C
C RMX 6400. ARROW SIZE ABOVE WHICH THE
C HEAD NO LONGER GROWS LARGER,
C ON A 2**15 X 2**15 GRID.
C
C SIDE 0.90 LENGTH OF LONGER EDGE OF PLOT.
C (SEE ALSO EXT.)
C
C SIZE 256. WIDTH OF THE CHARACTERS IN
C VECTOR LABELS, ON A 2**15 X
C 2**15 GRID.
C
C XLT 0.05 LEFT HAND EDGE OF THE PLOT.
C (0 IS THE LEFT EDGE OF THE
C FRAME, 1 THE RIGHT EDGE.)
C
C YBT 0.05 BOTTOM EDGE OF THE PLOT (0 IS
C THE BOTTOM OF THE FRAME, 1 THE
C TOP OF THE FRAME.)
C
C ---------------------------------------------------------------------
C
C INTERNAL FUNCTIONS WHICH MAY BE MODIFIED FOR DATA TRANSFORMATION -
C
C SCALE COMPUTES A SCALE FACTOR USED IN THE
C DETERMINATION OF THE LENGTH OF THE
C VECTOR TO BE DRAWN.
C
C DIST COMPUTES THE LENGTH OF A VECTOR.
C
C FX RETURNS THE X INDEX AS THE
C X-COORDINATE OF THE VECTOR BASE.
C
C MXF RETURNS THE X-COORDINATE OF THE VECTOR
C HEAD.
C
C FY RETURNS THE Y INDEX AS THE
C Y-COORDINATE OF THE VECTOR BASE.
C
C MYF RETURNS THE Y-COORDINATE OF THE VECTOR
C HEAD.
C
C VLAB THE VALUE FOR THE VECTOR LABEL WHEN
C ILAB IS NON-ZERO.
C
SAVE
FX(XX,YY) = XX
C FY(XX,YY) = YY
COLD
C FY(XX,YY)=YY+ITOP*((CTP(IFIX(XX))+(IFIX(XX)-XX)*(CTP(IFIX(XX))-
C 1CTP(IFIX(XX+1.))))*(Z0-YY))
FY(X,Y)=(1-ITOP)*Y + ITOP*(Y*( HTP(IFIX(X))
1 +(IFIX(X)-X)*(HTP(IFIX(X))-HTP(IFIX(X+1.))))
1+(CTP(IFIX(X))+(IFIX(X)-X)*(CTP(IFIX(X))-CTP(IFIX(X+1.))))*(Z0-Y))
DIST(XX,YY) = SQRT(XX*XX+YY*YY)
MXF(XX,YY,UU,VV,SFXX,SFYY,MXX,MYY) = MXX+IFIX(SFXX*UU)
MYF(XX,YY,UU,VV,SFXX,SFYY,MXX,MYY) = MYY+IFIX(SFYY*VV)
SCALEX(MM,NN,INCXX,INCYY,HAA,XX1,XX2,YY1,YY2,XX3,XX4,YY3,YY4,
1 LENN) = LENN/HAA
SCALEY(MM,NN,INCXX,INCYY,HAA,XX1,XX2,YY1,YY2,XX3,XX4,YY3,YY4,
1 LENN) = SCALEX(MM,NN,INCXX,INCYY,HAA,XX1,XX2,YY1,YY2,XX3,
2 XX4,YY3,YY4,LENN)
VLAB(UU,VV,II,JJ) = DIST(UU,VV)
C
C ---------------------------------------------------------------------
C
C THE FOLLOWING CALL IS FOR GATHERING STATISTICS ON LIBRARY USE AT NCAR.
C
C CALL Q8QST4 ('NSSL','VELVCT','VELVCT','VERSION 6')
C
C INITIALIZE AND TRANSFER SOME ARGUMENTS TO LOCAL VARIABLES.
C
LITARFL=1
BIG = -1.E-15
MX = LU
MY = LV
NX = M
NY = N
GL = FLO
HA = HI
ISP = 0
NC = 0
C
C COMPUTE CONSTANTS BASED ON THE ADDRESSABILITY OF THE PLOTTER.
C
CALL GETUSV('XF',ISX)
CALL GETUSV('YF',ISY)
ISX = 2**(15-ISX)
ISY = 2**(15-ISY)
LEN = LENGTH*ISX
C
C SET UP THE SCALING OF THE PLOT.
C
CALL GQCNTN(IERR,IOLDNT)
CALL GQNT(IOLDNT,IERR,WIND,VIEW)
X1 = VIEW(1)
X2 = VIEW(2)
Y1 = VIEW(3)
Y2 = VIEW(4)
X3 = WIND(1)
X4 = WIND(2)
Y3 = WIND(3)
Y4 = WIND(4)
CALL GETUSV('LS',IOLLS)
C
C SAVE NORMALIZATION TRANSFORMATION 1
C
CALL GQNT(1,IERR,WIND,VIEW)
C
IF (NSET) 101,102,106
C
101 X3 = 1.
X4 = FLOAT(NX)
Y3 = 1.
Y4 = FLOAT(NY)
GO TO 105
C
102 X1 = XLT
X2 = XLT+SIDE
Y1 = YBT
Y2 = YBT+SIDE
X3 = 1.
Y3 = 1.
X4 = FLOAT(NX)
Y4 = FLOAT(NY)
IF (AMIN1(X4,Y4)/AMAX1(X4,Y4) .LT. EXT) GO TO 105
C
IF (NX-NY) 103,105,104
103 X2 = XLT+SIDE*X4/Y4
GO TO 105
104 Y2 = YBT+SIDE*Y4/X4
C
105 CALL SET(X1,X2,Y1,Y2,X3,X4,Y3,Y4,1)
IF (NSET .EQ. 0) CALL PERIM (1,0,1,0)
C
C CALCULATE A LENGTH IF NONE PROVIDED.
C
106 IF (LEN .NE. 0) GO TO 107
CALL FL2INT(FX(1.,1.),FY(1.,1.),MX,MY)
CALL FL2INT(FX(FLOAT(1+INCX),FLOAT(1+INCY)),
+ FY(FLOAT(1+INCX),FLOAT(1+INCY)),LX,LY)
LEN = SQRT((FLOAT(MX-LX)**2+FLOAT(MY-LY)**2)/2.)
C
C SET UP SPECIAL VALUES.
C
107 IF (ISP .EQ. 0) GO TO 108
SPV1 = SPV(1)
SPV2 = SPV(2)
IF (ISP .EQ. 4) SPV2 = SPV(1)
C
C FIND THE MAXIMUM VECTOR LENGTH.
C
108 IF (HA .GT. 0.) GO TO 118
C
HA = BIG
IF (ISP .EQ. 0) GO TO 115
C
DO 114 J=IY0,NY-IY1,INCY
DO 113 I=IX0,NX-IX1,INCX
IF (ISP-2) 109,111,110
109 IF (U(I,J) .EQ. SPV1) GO TO 113
GO TO 112
110 IF (U(I,J) .EQ. SPV1) GO TO 113
111 IF (V(I,J) .EQ. SPV2) GO TO 113
112 HA = AMAX1(HA,DIST(U(I,J),V(I,J)))
113 CONTINUE
114 CONTINUE
GO TO 126
C
115 DO 117 J=IY0,NY-IY1,INCY
DO 116 I=IX0,NX-IX1,INCX
HA = AMAX1(HA,DIST(U(I,J),V(I,J)))
116 CONTINUE
117 CONTINUE
C
C BRANCH IF NULL VECTOR SIZE.
C
C 126 IF (HA .LE. 0.) GO TO 125
126 IF (HA .LE. 1.E-15) GO TO 125
C
C COMPUTE SCALE FACTORS.
C
118 SFX = SCALEX(M,N,INCX,INCY,HA,X1,X2,Y1,Y2,X3,X4,Y3,Y4,LEN)
SFY = SCALEY(M,N,INCX,INCY,HA,X1,X2,Y1,Y2,X3,X4,Y3,Y4,LEN)
IOFFDT = IOFFD
IF (GL.NE.0.0 .AND. (ABS(GL).LT.0.1 .OR. ABS(GL).GE.1.E5))
1 IOFFDT = 1
IF (HA.NE.0.0 .AND. (ABS(HA).LT.0.1 .OR. ABS(HA).GE.1.E5))
1 IOFFDT = 1
ASH = 1.0
IF (IOFFDT .NE. 0)
1 ASH = 10.**(3-IFIX(ALOG10(AMAX1(ABS(GL),ABS(HA)))-500.)-500)
IZFLG = 0
C
C COMPUTE ZMN AND ZMX, WHICH ARE USED IN DRWVEC.
C
ZMN = LEN*(GL/HA)
ZMX = FLOAT(LEN)+.01
C
C DRAW THE VECTORS.
C
DO 123 J=IY0,NY-IY1,INCY
DO 122 I=IX0,NX-IX1,INCX
UI = U(I,J)
VI = V(I,J)
IF (ISP-1) 121,119,120
119 IF (UI-SPV1) 121,122,121
120 IF (VI .EQ. SPV2) GO TO 122
IF (ISP .GE. 3) GO TO 119
121 X = I
Y = J
CALL FL2INT(FX(X,Y),FY(X,Y),MX,MY)
LX = MAX0(1,MXF(X,Y,UI,VI,SFX,SFY,MX,MY))
LY = MAX0(1,MYF(X,Y,UI,VI,SFX,SFY,MX,MY))
IZFLG = 1
IF (ILAB .NE. 0) CALL ENCD(VLAB(UI,VI,I,J),ASH,LABEL,NC,
+ IOFFDT)
CALL DRWVEC (MX,MY,LX,LY,LABEL,NC)
122 CONTINUE
123 CONTINUE
C
IF (IZFLG .EQ. 0) GO TO 125
C
IF (IOFFM .NE. 0) GO TO 200
IF(IVRT.EQ.0) THEN
write (label,904) ha
904 FORMAT(F7.2,5H M/S)
ENDIF
IF(IVRT.EQ.1) THEN
write (label,9041) ha
9041 FORMAT(E10.3,4H vrt)
ENDIF
C
C TURN OFF CLIPPING SO ARROW CAN BE DRAWN
C
CALL GQCLIP(IER,ICLP,IAR)
CALL GSCLIP(0)
IF(LITARFL.EQ.1) THEN
#if (COLORPL == 1)
CALL DRWVEC (24768,3608,24768+LEN,3608,LABEL,14)
#else
CALL DRWVEC (28768,384,28768+LEN,384,LABEL,14)
#endif
ENDIF
C RESTORE CLIPPING
C
CALL GSCLIP(ICLP)
IX = 1+(28768+LEN/2)/ISX
IY = 1+(608-(5*ISX*MAX0(256/ISX,8))/4)/ISY
CALL GQCNTN(IER,ICN)
CALL GSELNT(0)
C XC = CPUX(IX)
C YC = CPUY(IY)
C CALL WTSTR (XC,YC,
C + 'MAXIMUM VECTOR',MAX0(256/ISX,8),0,0)
CALL GSELNT(ICN)
C
C DONE.
C
GOTO 200
C
C ZERO-FIELD ACTION.
C
125 IX = 1+16384/ISX
IY = 1+16384/ISY
CALL GQCNTN(IER,ICN)
CALL GSELNT(0)
XC = CPUX(IX)
YC = CPUY(IY)
CALL WTSTR (XC,YC,
+ 'ZERO FIELD',MAX0(960/ISX,8),0,0)
CALL GSELNT(ICN)
C
C RESTORE TRANS 1 AND LOG SCALING AND ORIGINAL TRANS NUMBER
C
200 CONTINUE
IF (NSET .LE. 0) THEN
CALL SET(VIEW(1),VIEW(2),VIEW(3),VIEW(4),
- WIND(1),WIND(2),WIND(3),WIND(4),IOLLS)
ENDIF
CALL GSELNT(IOLDNT)
RETURN
END
SUBROUTINE DRWVEC (M1,M2,M3,M4,LABEL,NC)
C
C THIS ROUTINE IS CALLED TO DRAW A SINGLE ARROW. IT HAS ARGUMENTS AS
C FOLLOWS -
C
C (M1,M2) - COORDINATE OF ARROW BASE, ON A 2**15 X 2**15 GRID.
C (M3,M4) - COORDINATE OF ARROW HEAD, ON A 2**15 X 2**15 GRID.
C LABEL - CHARACTER LABEL TO BE PUT ABOVE ARROW.
C NC - NUMBER OF CHARACTERS IN LABEL.
C
SAVE
C
C
COMMON /VEC1/ ASH ,EXT ,ICTRFG ,ILAB ,
+ IOFFD ,IOFFM ,ISX ,ISY ,
+ RMN ,RMX ,SIDE ,SIZE ,
+ XLT ,YBT ,ZMN ,ZMX
CHARACTER*14 LABEL
C
C SOME LOCAL PARAMETERS ARE THE FOLLOWING -
C
C CL - ARROW HEAD LENGTH SCALE FACTOR - EACH SIDE OF THE ARROW
C HEAD IS THIS LONG RELATIVE TO THE LENGTH OF THE ARROW
C ST,CT - SIN AND COS OF THE ARROW HEAD ANGLE
C PI - THE CONSTANT PI
C TWOPI - TWO TIMES PI
C OHOPI - ONE HALF OF PI
C FHOPI - FIVE HALVES OF PI
C
DATA CL / .25 /
DATA ST / .382683432365090 /
DATA CT / .923879532511287 /
DATA PI / 3.14159265358979 /
DATA TWOPI / 6.28318530717959 /
DATA OHOPI / 1.57079632679489 /
DATA FHOPI / 7.85398163397448 /
C
DIST(X,Y) = SQRT(X*X+Y*Y)
C
C TRANSFER ARGUMENTS TO LOCAL VARIABLES AND COMPUTE THE VECTOR LENGTH.
C
N1 = M1
N2 = M2
N3 = M3
N4 = M4
DX = N3-N1
DY = N4-N2
R = DIST(DX,DY)
C
C SORT OUT POSSIBLE CASES, DEPENDING ON VECTOR LENGTH.
C
IF (R .LE. ZMN) RETURN
C
c IF (R .LE. ZMX) GO TO 101
GO TO 101
C
C PLOT A POINT FOR VECTORS WHICH ARE TOO LONG.
C
c CALL PLOTIT (N1,N2,0)
c CALL PLOTIT (N1,N2,1)
c CALL PLOTIT (N1,N2,0)
c RETURN
C
C ADJUST THE COORDINATES OF THE VECTOR ENDPOINTS AS IMPLIED BY THE
C CENTERING OPTION.
C
101 IF (ICTRFG) 102,103,104
C
102 N3 = N1
N4 = N2
N1 = FLOAT(N1)-DX
N2 = FLOAT(N2)-DY
GO TO 104
C
103 N1 = FLOAT(N1)-.5*DX
N2 = FLOAT(N2)-.5*DY
N3 = FLOAT(N3)-.5*DX
N4 = FLOAT(N4)-.5*DY
C
C DETERMINE THE COORDINATES OF THE POINTS USED TO DRAW THE ARROWHEAD.
C
104 C1 = CL
C
C SHORT ARROWS HAVE HEADS OF A FIXED MINIMUM SIZE.
C
IF (R .LT. RMN) C1 = RMN*CL/R
C
C LONG ARROWS HAVE HEADS OF A FIXED MAXIMUM SIZE.
C
IF (R .GT. RMX) C1 = RMX*CL/R
C
C COMPUTE THE COORDINATES OF THE HEAD.
C
N5 = FLOAT(N3)-C1*(CT*DX-ST*DY)
N6 = FLOAT(N4)-C1*(CT*DY+ST*DX)
N7 = FLOAT(N3)-C1*(CT*DX+ST*DY)
N8 = FLOAT(N4)-C1*(CT*DY-ST*DX)
C
C PLOT THE ARROW.
C
CALL PLOTIT (N1,N2,0)
CALL PLOTIT (N3,N4,1)
CALL PLOTIT (N5,N6,0)
CALL PLOTIT (N3,N4,1)
CALL PLOTIT (N7,N8,1)
CALL PLOTIT (0,0,0)
C
C IF REQUESTED, PUT THE VECTOR MAGNITUDE ABOVE THE ARROW.
C
IF (NC .EQ. 0) RETURN
PHI = ATAN2(DY,DX)
IF (AMOD(PHI+FHOPI,TWOPI) .GT. PI) PHI = PHI+PI
IX = 1+IFIX(.5*FLOAT(N1+N3)+1.25*
+ FLOAT(ISX*MAX0(IFIX(SIZE)/ISX,8))*COS(PHI+OHOPI))/ISX
IY = 1+IFIX(.5*FLOAT(N2+N4)+1.25*
+ FLOAT(ISX*MAX0(IFIX(SIZE)/ISX,8))*SIN(PHI+OHOPI))/ISY
CALL GQCNTN(IER,ICN)
CALL GSELNT(0)
XC = CPUX(IX)
YC = CPUY(IY)
CALL WTSTR(XC,YC,
+ LABEL,MAX0(IFIX(SIZE)/ISX,8),
+ IFIX(57.2957795130823*PHI),0)
CALL GSELNT(ICN)
RETURN
END
SUBROUTINE VELVEC (U,LU,V,LV,M,N,FLO,HI,NSET)
C
C THIS ROUTINE SUPPORTS USERS OF THE OLD VERSION OF THIS PACKAGE.
C
DIMENSION U(LU,N) ,V(LV,N) ,SPV(2)
C
SAVE
C
C THE FOLLOWING CALL IS FOR GATHERING STATISTICS ON LIBRARY USE AT NCAR.
C
C CALL Q8QST4 ('CRAYLIB','VELVCT','VELVEC','VERSION 4')
CALL VELVCT (U,LU,V,LV,M,N,FLO,HI,NSET,0)
RETURN
END
BLOCK DATA VELDAT
C
C THIS 'ROUTINE' DEFINES THE DEFAULT VALUES OF THE VELVCT PARAMETERS.
C
SAVE
C
COMMON /VEC1/ ASH ,EXT ,ICTRFG ,ILAB ,
+ IOFFD ,IOFFM ,ISX ,ISY ,
+ RMN ,RMX ,SIDE ,SIZE ,
+ XLT ,YBT ,ZMN ,ZMX
C
COMMON /VEC2/ BIG,IX0,IX1,INCX,IY0,IY1,INCY
C
DATA EXT / 0.25 /
DATA ICTRFG / 1 /
DATA ILAB / 0 /
DATA IOFFD / 0 /
DATA IOFFM / 0 /
DATA RMN / 160.00 /
DATA RMX / 6400.00 /
DATA SIDE / 0.90 /
DATA SIZE / 256.00 /
DATA XLT / 0.05 /
DATA YBT / 0.05 /
DATA ZMX / 0.00 /
DATA IX0,IX1,INCX /1,0, 1 /
DATA IY0,IY1,INCY /1,0, 1 /
C
C REVISION HISTORY ----------------------------------------------------
C
C FEBRUARY, 1979 ADDED REVISION HISTORY
C MODIFIED CODE TO CONFORM TO FORTRAN 66 STANDARD
C
C JULY, 1979 FIXED HI VECTOR TRAP AND MESSAGE INDICATING
C MAXIMUM VECTOR PLOTTED.
C
C DECEMBER, 1979 CHANGED THE STATISTICS CALL FROM CRAYLIB TO NSSL
C
C MARCH, 1981 FIXED SOME FRINGE-CASE ERRORS, CHANGED THE CODE TO
C USE FL2INTT AND PLOTIT INSTEAD OF MXMY, FRSTPT, AND
C VECTOR, AND MADE THE ARROWHEADS NARROWER (45 DEGREES
C APART, RATHER THAN 60 DEGREES APART)
C
C FEBRUARY, 1984 PROVIDED A DIMENSION STATEMENT FOR A VARIABLE INTO
C WHICH A TEN-CHARACTER STRING WAS BEING ENCODED. ON
C THE CRAY, WHEN THE ENCODE WAS DONE, A WORD FOLLOWING
C THE VARIABLE WAS CLOBBERED, BUT THIS APPARENTLY MADE
C NO DIFFERENCE. ON AT LEAST ONE OTHER MACHINE, THE
C CODE BLEW UP. (ERROR REPORTED BY GREG WOODS)
C
C JULY, 1984 CONVERTED TO FORTRAN77 AND GKS.
C
C ---------------------------------------------------------------------
END
SUBROUTINE STRMLN (U,V,WORK,IMAX,IPTSX,JPTSY,NSET,IER)
C
C +-----------------------------------------------------------------+
C | |
C | Copyright (C) 1989 by UCAR |
C | University Corporation for Atmospheric Research |
C | All Rights Reserved |
C | |
C | NCARGRAPHICS Version 3.00 |
C | |
C +-----------------------------------------------------------------+
C
C SUBROUTINE STRMLN (U,V,WORK,IMAX,IPTSX,JPTSY,NSET,IER)
C
C DIMENSION OF U(IMAX,JPTSY) , V(IMAX,JPTSY) ,
C ARGUMENTS WORK(2*IMAX*JPTSY)
C
C PURPOSE STRMLN draws a streamline representation of
C the flow field. The representation is
C independent of the flow speed.
C
C USAGE If the following assumptions are met, use
C
C CALL EZSTRM (U,V,WORK,IMAX,JMAX)
C
C Assumptions:
C --The whole array is to be processed.
C --The arrays are dimensioned
C U(IMAX,JMAX) , V(IMAX,JMAX) and
C WORK(2*IMAX*JMAX).
C --Window and viewport are to be chosen
C by STRMLN.
C --PERIM is to be called.
C
C If these assumptions are not met, use
C
C CALL STRMLN (U,V,WORK,IMAX,IPTSX,JPTSY,
C NSET,IER)
C
C The user must call FRAME in the calling
C routine.
C
C The user may change various internal
C parameters via common blocks. See below.
C
C ARGUMENTS
C
C ON INPUT U, V
C Two dimensional arrays containing the
C velocity fields to be plotted.
C
C Note: If the U AND V components
C are, for example, defined in Cartesian
C coordinates and the user wishes to plot them
C on a different projection (i.e., stereo-
C graphic), then the appropriate
C transformation must be made to the U and V
C components via the functions FU and FV
C (located in DRWSTR).
C
C WORK
C User provided work array. The dimension
C of this array must be .GE. 2*IMAX*JPTSY.
C
C Caution: This routine does not check the
C size of the work array.
C
C IMAX
C The first dimension of U and V in the
C calling program. (X-direction)
C
C IPTSX
C The number of points to be plotted in the
C first subscript direction. (X-direction)
C
C JPTSY
C The number of points to be plotted in the
C second subscript direction. (Y-direction)
C
C NSET
C Flag to control scaling
C > 0 STRMLN assumes that the window
C and viewport have been set by the
C user in such a way as to properly
C scale the plotting instructions
C generated by STRMLN. PERIM is not
C called.
C = 0 STRMLN will establish the window and
C viewport to properly scale the
C plotting instructions to the standard
C configuration. PERIM is called to draw
C the border.
C < 0 STRMLN establishes the window
C and viewport so as to place the
C streamlines within the limits
C of the user's window. PERIM is
C not called.
C
C ON OUTPUT Only the IER argument may be changed. All
C other arguments are unchanged.
C
C
C IER
C = 0 when no errors are detected
C = -1 when the routine is called with ICYC
C .NE. 0 and the data are not cyclic
C (ICYC is an internal parameter
C described below); in this case the
C routine will draw the
C streamlines with the non-cyclic
C interpolation formulas.
C
C ENTRY POINTS STRMLN, DRWSTR, EZSTRM, GNEWPT, CHKCYC
C
C COMMON BLOCKS STR01, STR02, STR03, STR04
C
C REQUIRED LIBRARY GRIDAL, GBYTES, and the SPPS
C ROUTINES
C
C REQUIRED GKS LEVEL 0A
C
C I/O None
C
C PRECISION Single
C
C LANGUAGE FORTRAN 77
C
C HISTORY Written and standardized in November 1973.
C
C Converted to FORTRAN 77 and GKS in June, 1984.
C
C
C PORTABILITY FORTRAN 77
C
C ALGORITHM Wind components are normalized to the value
C of DISPL. The least significant two
C bits of the work array are
C utilized as flags for each grid box. Flag 1
C indicates whether any streamline has
C previously passed through this box. Flag 2
C indicates whether a directional arrow has
C already appeared in a box. Judicious use
C of these flags prevents overcrowding of
C streamlines and directional arrows.
C Experience indicates that a final pleasing
C picture is produced when streamlines are
C initiated in the center of a grid box. The
C streamlines are drawn in one direction then
C in the opposite direction.
C
C REFERENCE The techniques utilized here are described
C in an article by Thomas Whittaker (U. of
C Wisconsin) which appeared in the notes and
C correspondence section of Monthly Weather
C Review, June 1977.
C
C TIMING Highly variable
C It depends on the complexity of the
C flow field and the parameters: DISPL,
C DISPC , CSTOP , INITA , INITB , ITERC ,
C and IGFLG. (See below for a discussion
C of these parameters.) If all values
C are default, then a simple linear
C flow field for a 40 x 40 grid will
C take about 0.4 seconds on the CRAY1-A;
C a fairly complex flow field will take about
C 1.5 seconds on the CRAY1-A.
C
C
C INTERNAL PARAMETERS
C
C NAME DEFAULT FUNCTION
C ---- ------- --------
C
C EXT 0.25 Lengths of the sides of the
C plot are proportional to
C IPTSX and JPTSY except in
C the case when MIN(IPTSX,JPT
C / MAX(IPTSX,JPTSY) .LT. EXT
C in that case a square
C graph is plotted.
C
C SIDE 0.90 Length of longer edge of
C plot. (See also EXT.)
C
C XLT 0.05 Left hand edge of the plot.
C (0.0 = left edge of frame)
C (1.0 = right edge of frame)
C
C YBT 0.05 Bottom edge of the plot.
C (0.0 = bottom ; 1.0 = top)
C
C (YBT+SIDE and XLT+SIDE must
C be .LE. 1. )
C
C INITA 2 Used to precondition grid
C boxes to be eligible to
C start a streamline.
C For example, a value of 4
C means that every fourth
C grid box is eligible ; a
C value of 2 means that every
C other grid box is eligible.
C (see INITB)
C
C INITB 2 Used to precondition grid
C boxes to be eligible for
C direction arrows.
C If the user changes the
C default values of INITA
C and/or INITB, it should
C be done such that
C MOD(INITA,INITB) = 0 .
C For a dense grid try
C INITA=4 and INITB=2 to
C reduce the CPU time.
C
C AROWL 0.33 Length of direction arrow.
C For example, 0.33 means
C each directional arrow will
C take up a third of a grid
C box.
C
C ITERP 35 Every 'ITERP' iterations
C the streamline progress
C is checked.
C
C ITERC -99 The default value of this
C parameter is such that
C it has no effect on the
C code. When set to some
C positive value, the program
C will check for streamline
C crossover every 'ITERC'
C iterations. (The routine
C currently does this every
C time it enters a new grid
C box.)
C Caution: When this
C parameter is activated,
C CPU time will increase.
C
C IGFLG 0 A value of zero means that
C the sixteen point Bessel
C Interpolation Formula will
C be utilized where possible;
C when near the grid edges,
C quadratic and bi-linear
C interpolation will be
C used. This mixing of
C interpolation schemes can
C sometimes cause slight
C raggedness near the edges
C of the plot. If IGFLG.NE.0
C then only the bilinear
C interpolation formula
C is used; this will generall
C result in slightly faster
C plot times but a less
C pleasing plot.
C
C IMSG 0 If zero, then no missing
C U and V components are
C present.
C If .NE. 0, STRMLN will
C utilize the
C bi-linear interpolation
C scheme and terminate if
C any data points are missing
C
C UVMSG 1.E+36 Value assigned to a missing
C point.
C
C ICYC 0 Zero means the data are
C non-cyclic in the X
C direction.
C If .NE 0, the
C cyclic interpolation
C formulas will be used.
C (Note: Even if the data
C are cyclic in X, leaving
C ICYC = 0 will do no harm.)
C
C DISPL 0.33 The wind speed is
C normalized to this value.
C (See the discussion below.)
C
C DISPC 0.67 The critical displacement.
C If after 'ITERP' iterations
C the streamline has not
C moved this distance, the
C streamline will be
C terminated.
C
C CSTOP 0.50 This parameter controls
C the spacing between
C streamlines. The checking
C is done when a new grid
C box is entered.
C
C DISCUSSION OF Assume a value of 0.33 for DISPL. This
C DISPL,DISPC means that it will take three steps to move
C AND CSTOP across one grid box if the flow was all in the
C X direction. If the flow is zonal, then a
C larger value of DISPL is in order.
C If the flow is highly turbulent, then
C a smaller value is in order. The smaller
C DISPL, the more the CPU time. A value
C of 2 to 4 times DISPL is a reasonable value
C for DISPC. DISPC should always be greater
C than DISPL. A value of 0.33 for CSTOP would
C mean that a maximum of three stream-
C lines will be drawn per grid box. This max
C will normally only occur in areas of singular
C points.
C
C ***************************
C Any or all of the above
C parameters may be changed
C by utilizing common blocks
C STR02 and/or STR03
C ***************************
C
C UXSML A number which is small
C compared to the average
C normalized u component.
C Set automatically.
C
C NCHK 750 This parameter is located
C in DRWSTR. It specifies the
C length of the circular
C lists used for checking
C for STRMLN crossovers.
C For most plots this number
C may be reduced to 500
C or less and the plots will
C not be altered.
C
C ISKIP Number of bits to be
C skipped to get to the
C least two significant bits
C in a floating point number.
C The default value is set to
C I1MACH(5) - 2 . This value
C may have to be changed
C depending on the target
C computer; see subroutine
C DRWSTR.
C
C
C
DIMENSION U(IMAX,JPTSY) ,V(IMAX,JPTSY) ,
1 WORK(1)
DIMENSION WNDW(4) ,VWPRT(4)
C
COMMON /STR01/ IS ,IEND ,JS ,JEND
1 , IEND1 ,JEND1 ,I ,J
2 , X ,Y ,DELX ,DELY
3 , ICYC1 ,IMSG1 ,IGFL1
COMMON /STR02/ EXT , SIDE , XLT , YBT
CIBM8 common /str03/ arowl,uvmsg,displ,dispc,cstop,
CIBM81 inita,initb,iterp,iterc,igflg,imsg,icyc
COMMON /STR03/ INITA , INITB , AROWL , ITERP , ITERC , IGFLG
1 , IMSG , UVMSG , ICYC , DISPL , DISPC , CSTOP
C
SAVE
C
EXT = 0.25
SIDE = 0.90
XLT = 0.05
YBT = 0.05
C
C INITA = 2
C INITB = 2
AROWL = 0.40
C AROWL = 0.33
ITERP = 35
ITERC = -99
IGFLG = 0
ICYC = 0
IMSG = 0
UVMSG = 1.E+36
DISPL = 0.33
DISPC = 0.67
c CSTOP = 0.50
C
C THE FOLLOWING CALL IS FOR MONITORING LIBRARY USE AT NCAR
C
C CALL Q8QST4 ( 'GRAPHX', 'STRMLN', 'STRMLN', 'VERSION 01')
C
IER = 0
C
C LOAD THE COMMUNICATION COMMON BLOCK WITH PARAMETERS
C
IS = 1
IEND = IPTSX
JS = 1
JEND = JPTSY
IEND1 = IEND-1
JEND1 = JEND-1
IEND2 = IEND-2
JEND2 = JEND-2
XNX = FLOAT(IEND-IS+1)
XNY = FLOAT(JEND-JS+1)
ICYC1 = ICYC
IGFL1 = IGFLG
IMSG1 = 0
C
C IF ICYC .NE. 0 THEN CHECK TO MAKE SURE THE CYCLIC CONDITION EXISTS.
C
IF (ICYC1.NE.0) CALL CHKCYC (U,V,IMAX,JPTSY,IER)
C
C SAVE ORIGINAL NORMALIZATION TRANSFORMATION NUMBER
C
CALL GQCNTN ( IERR,NTORIG )
C
C SET UP SCALING
C
IF (NSET) 10 , 20 , 60
10 CALL GETUSV ( 'LS' , ITYPE )
CALL GQNT ( NTORIG,IERR,WNDW,VWPRT )
CALL GETUSV('LS',IOLLS)
X1 = VWPRT(1)
X2 = VWPRT(2)
Y1 = VWPRT(3)
Y2 = VWPRT(4)
X3 = IS
X4 = IEND
Y3 = JS
Y4 = JEND
GO TO 55
C
20 ITYPE = 1
X1 = XLT
X2 = (XLT+SIDE)
Y1 = YBT
Y2 = (YBT+SIDE)
X3 = IS
X4 = IEND
Y3 = JS
Y4 = JEND
IF (AMIN1(XNX,XNY)/AMAX1(XNX,XNY).LT.EXT) GO TO 50
IF (XNX-XNY) 30, 50, 40
30 X2 = (SIDE*(XNX/XNY) + XLT)
GO TO 50
40 Y2 = (SIDE*(XNY/XNX) + YBT)
50 CONTINUE
C
C CENTER THE PLOT
C
DX = 0.25*( 1. - (X2-X1) )
DY = 0.25*( 1. - (Y2-Y1) )
X1 = (XLT+DX)
X2 = (X2+DX )
Y1 = (YBT+DY)
Y2 = (Y2+DY )
C
55 CONTINUE
C
C SAVE NORMALIZATION TRANSFORMATION 1
C
CALL GQNT ( 1,IERR,WNDW,VWPRT )
C
C DEFINE AND SELECT NORMALIZATION TRANS, SET LOG SCALING
C
CALL SET(X1,X2,Y1,Y2,X3,X4,Y3,Y4,ITYPE)
C
IF (NSET.EQ.0) CALL PERIM (1,0,1,0)
C
60 CONTINUE
C
C DRAW THE STREAMLINES
C . BREAK THE WORK ARRAY INTO TWO PARTS. SEE DRWSTR FOR FURTHER
C . COMMENTS ON THIS.
C
CALL DRWSTR (U,V,WORK(1),WORK(IMAX*JPTSY+1),IMAX,JPTSY)
C
C RESET NORMALIATION TRANSFORMATION 1 TO ORIGINAL VALUES
C
IF (NSET .LE. 0) THEN
CALL SET(VWPRT(1),VWPRT(2),VWPRT(3),VWPRT(4),
- WNDW(1),WNDW(2),WNDW(3),WNDW(4),IOLLS)
ENDIF
CALL GSELNT (NTORIG)
C
RETURN
END
SUBROUTINE DRWSTR (U,V,UX,VY,IMAX,JPTSY)
C
PARAMETER (NCHK=750)
C
C THIS ROUTINE DRAWS THE STREAMLINES.
C . THE XCHK AND YCHK ARRAYS SERVE AS A CIRCULAR LIST. THEY
C . ARE USED TO PREVENT LINES FROM CROSSING ONE ANOTHER.
C
C THE WORK ARRAY HAS BEEN BROKEN UP INTO TWO ARRAYS FOR CLARITY. THE
C . TOP HALF OF WORK (CALLED UX) WILL HAVE THE NORMALIZED (AND
C . POSSIBLY TRANSFORMED) U COMPONENTS AND WILL BE USED FOR BOOK
C . KEEPING. THE LOWER HALF OF THE WORK ARRAY (CALLED VY) WILL
C . CONTAIN THE NORMALIZED (AND POSSIBLY TRANSFORMED) V COMPONENTS.
C
DIMENSION U(IMAX,JPTSY) ,V(IMAX,JPTSY)
1 , UX(IMAX,JPTSY) ,VY(IMAX,JPTSY)
COMMON /STR01/ IS ,IEND ,JS ,JEND
1 , IEND1 ,JEND1 ,I ,J
2 , X ,Y ,DELX ,DELY
3 , ICYC1 ,IMSG1 ,IGFL1
CIBM8 common /str03/ arowl,uvmsg,displ,dispc,cstop,
CIBM81 inita,initb,iterp,iterc,igflg,imsg,icyc
COMMON /STR03/ INITA , INITB , AROWL , ITERP , ITERC , IGFLG
1 , IMSG , UVMSG , ICYC , DISPL , DISPC , CSTOP
COMMON /STR04/ XCHK(NCHK) ,YCHK(NCHK) , UXSML, NUMCHK
COMMON/TOPOG/ CTP(2000),HTP(2000),Z0,ITOP
C
C
SAVE
C
C STATEMENT FUNCTIONS FOR SPATIAL AND VELOCITY TRANSFORMATIONS.
C . (IF THE USER WISHES OTHER TRANSFORMATIONS REPLACE THESE STATEMENT
C . FUNCTIONS WITH THE APPROPRIATE NEW ONES, OR , IF THE TRANSFORMA-
C . TIONS ARE COMPLICATED DELETE THESE STATEMENT FUNCTIONS
C . AND ADD EXTERNAL ROUTINES WITH THE SAME NAMES TO DO THE TRANS-
C . FORMING.)
C
FX(X,Y) = X
c FY(X,Y) = Y
COLD
C FY(XX,YY)=YY+ITOP*((CTP(IFIX(XX))+(IFIX(XX)-XX)*(CTP(IFIX(XX))-
C 1CTP(IFIX(XX+1.))))*(Z0-YY))
FY(X,Y)=(1-ITOP)*Y + ITOP*(Y*( HTP(IFIX(X))
1 +(IFIX(X)-X)*(HTP(IFIX(X))-HTP(IFIX(X+1.))))
1+(CTP(IFIX(X))+(IFIX(X)-X)*(CTP(IFIX(X))-CTP(IFIX(X+1.))))*(Z0-Y))
FU(X,Y) = X
FV(X,Y) = Y
C
C INITIALIZE
C
ISKIP = I1MACH(5) - 2
ISKIP1 = ISKIP + 1
UXSML = 1.E-30
C
C
NUMCHK = NCHK
LCHK = 1
ICHK = 1
XCHK(1) = 0.
YCHK(1) = 0.
KFLAG = 0
IZERO = 0
IONE = 1
ITWO = 2
C
C
C COMPUTE THE X AND Y NORMALIZED (AND POSSIBLY TRANSFORMED)
C . DISPLACEMENT COMPONENTS (UX AND VY).
C
DO 40 J=JS,JEND
DO 30 I=IS,IEND
UX(I,J) = FU(U(I,J),V(I,J))
VY(I,J) = FV(U(I,J),V(I,J))
IF (UX(I,J).NE.0. .OR. VY(I,J).NE.0.) THEN
CON = DISPL/SQRT(UX(I,J)*UX(I,J) + VY(I,J)*VY(I,J))
UX(I,J) = CON*UX(I,J)
VY(I,J) = CON*VY(I,J)
END IF
C
C BOOKKEEPING IS DONE IN THE LEAST SIGNIFICANT BITS OF THE UX ARRAY.
C . WHEN UX(I,J) IS EXACTLY ZERO THIS CAN PRESENT SOME PROBLEMS.
C . TO GET AROUND THIS PROBLEM, SET IT TO A RELATIVELY SMALL NUMBER.
C
IF(UX(I,J) .EQ. 0.) UX(I,J) = UXSML
C
C MASK OUT THE LEAST SIGNIFICANT TWO BITS AS FLAGS FOR EACH GRID BOX
C . A GRID BOX IS ANY REGION SURROUNDED BY FOUR GRID POINTS.
C . FLAG 1 INDICATES WHETHER ANY STREAMLINE HAS PREVIOUSLY PASSED
C . THROUGH THIS BOX.
C . FLAG 2 INDICATES WHETHER ANY DIRECTIONAL ARROW HAS ALREADY
C . APPEARED IN THIS BOX.
C . JUDICIOUS USE OF THESE FLAGS PREVENTS OVERCROWDING OF
C . STREAMLINES AND DIRECTIONAL ARROWS.
C
CALL SBYTES( UX(I,J) , IZERO , ISKIP , 2 , 0 , 1 )
C
IF (MOD(I,INITA).NE.0 .OR. MOD(J,INITA).NE.0)
1 CALL SBYTES( UX(I,J) , IONE , ISKIP1, 1 , 0 , 1 )
IF (MOD(I,INITB).NE.0 .OR. MOD(J,INITB).NE.0)
1 CALL SBYTES( UX(I,J) , IONE , ISKIP , 1 , 0 , 1 )
C
30 CONTINUE
40 CONTINUE
C
50 CONTINUE
C
C START A STREAMLINE. EXPERIENCE HAS SHOWN THAT A PLEASING PICTURE
C . WILL BE PRODUCED IF NEW STREAMLINES ARE STARTED ONLY IN GRID
C . BOXES THAT PREVIOUSLY HAVE NOT HAD OTHER STREAMLINES PASS THROUGH
C . THEM. AS LONG AS A REASONABLY DENSE PATTERN OF AVAILABLE BOXES
C . IS INITIALLY PRESCRIBED, THE ORDER OF SCANNING THE GRID PTS. FOR
C . AVAILABLE BOXES IS IMMATERIAL
C
C FIND AN AVAILABLE BOX FOR STARTING A STREAMLINE
C
IF (KFLAG.NE.0) GO TO 90
DO 70 J=JS,JEND1
DO 60 I=IS,IEND1
CALL GBYTES( UX(I,J) , IUX , ISKIP , 2 , 0 , 1 )
IF ( IAND( IUX , IONE ) .EQ. IZERO ) GO TO 80
60 CONTINUE
70 CONTINUE
C
C MUST BE NO AVAILABLE BOXES FOR STARTING A STREAMLINE
C
GO TO 190
80 CONTINUE
C
C INITILIZE PARAMETERS FOR STARTING A STREAMLINE
C . TURN THE BOX OFF FOR STARTING A STREAMLINE
C . CHECK TO SEE IF THIS BOX HAS MISSING DATA (IMSG.NE.0). IF SO ,
C . FIND A NEW STARTING BOX
C
CALL SBYTES( UX(I,J) , IONE , ISKIP1 , 1 , 0 , 1 )
IF ( IMSG.EQ.0) GO TO 85
IF (U(I,J).EQ.UVMSG .OR. U(I,J+1).EQ.UVMSG .OR.
1 U(I+1,J).EQ.UVMSG .OR. U(I+1,J+1).EQ.UVMSG) GO TO 50
C
85 ISAV = I
JSAV = J
KFLAG = 1
PLMN1 = +1.
GO TO 100
90 CONTINUE
C
C COME TO HERE TO DRAW IN THE OPPOSITE DIRECTION
C
KFLAG = 0
PLMN1 = -1.
I = ISAV
J = JSAV
100 CONTINUE
C
C INITIATE THE DRAWING SEQUENCE
C . START ALL STREAMLINES IN THE CENTER OF A BOX
C
NBOX = 0
ITER = 0
IF (KFLAG.NE.0) ICHKB = ICHK+1
IF (ICHKB.GT.NUMCHK) ICHKB = 1
X = FLOAT(I)+0.5
Y = FLOAT(J)+0.5
XBASE = X
YBASE = Y
CALL FL2INT (FX(X,Y),FY(X,Y),IFX,IFY)
CALL PLOTIT (IFX,IFY,0)
CALL GBYTES( UX(I,J) , IUX , ISKIP , 2 , 0 , 1 )
IF ( (KFLAG.EQ.0) .OR. (IAND( IUX , ITWO ) .NE. 0 ) ) GO TO 110
C
C GRID BOX MUST BE ELIGIBLE FOR A DIRECTIONAL ARROW
C
CALL GNEWPT (UX,VY,IMAX,JPTSY)
MFLAG = 1
GO TO 160
C
110 CONTINUE
C
C PLOT LOOP
C . CHECK TO SEE IF THE STREAMLINE HAS ENTERED A NEW GRID BOX
C
IF (I.NE.IFIX(X) .OR. J.NE.IFIX(Y)) GO TO 120
C
C MUST BE IN SAME BOX CALCULATE THE DISPLACEMENT COMPONENTS
C
CALL GNEWPT (UX,VY,IMAX,JPTSY)
C
C UPDATE THE POSITION AND DRAW THE VECTOR
C
X = X+PLMN1*DELX
Y = Y+PLMN1*DELY
CALL FL2INT (FX(X,Y),FY(X,Y),IFX,IFY)
CALL PLOTIT (IFX,IFY,1)
ITER = ITER+1
C
C CHECK STREAMLINE PROGRESS EVERY 'ITERP' OR SO ITERATIONS
C
IF (MOD(ITER,ITERP).NE.0) GO TO 115
IF (ABS(X-XBASE).LT.DISPC .AND. ABS(Y-YBASE).LT.DISPC ) GO TO 50
XBASE = X
YBASE = Y
GO TO 110
115 CONTINUE
C
C SHOULD THE CIRCULAR LISTS BE CHECKED FOR STREAMLINE CROSSOVER
C
IF ( (ITERC.LT.0) .OR. (MOD(ITER,ITERC).NE.0) ) GO TO 110
C
C MUST WANT THE CIRCULAR LIST CHECKED
C
GO TO 130
120 CONTINUE
C
C MUST HAVE ENTERED A NEW GRID BOX CHECK FOR THE FOLLOWING :
C . (1) ARE THE NEW POINTS ON THE GRID
C . (2) CHECK FOR MISSING DATA IF MSG DATA FLAG (IMSG) HAS BEEN SET.
C . (3) IS THIS BOX ELIGIBLE FOR A DIRECTIONAL ARROW
C . (4) LOCATION OF THIS ENTRY VERSUS OTHER STREAMLINE ENTRIES
C
NBOX = NBOX+1
C
C CHECK (1)
C
IF (IFIX(X).LT.IS .OR. IFIX(X).GT.IEND1) GO TO 50
IF (IFIX(Y).LT.JS .OR. IFIX(Y).GT.JEND1) GO TO 50
C
C CHECK (2)
C
IF ( IMSG.EQ.0) GO TO 125
II = IFIX(X)
JJ = IFIX(Y)
IF (U(II,JJ).EQ.UVMSG .OR. U(II,JJ+1).EQ.UVMSG .OR.
1 U(II+1,JJ).EQ.UVMSG .OR. U(II+1,JJ+1).EQ.UVMSG) GO TO 50
125 CONTINUE
C
C CHECK (3)
C
CALL GBYTES( UX(I,J) , IUX , ISKIP , 2 , 0 , 1 )
IF ( IAND( IUX , ITWO ) .NE. 0) GO TO 130
MFLAG = 2
GO TO 160
130 CONTINUE
C
C CHECK (4)
C
DO 140 LOC=1,LCHK
IF (ABS( X-XCHK(LOC) ).GT.CSTOP .OR.
1 ABS( Y-YCHK(LOC) ).GT.CSTOP) GO TO 140
LFLAG = 1
IF (ICHKB.LE.ICHK .AND. LOC.GE.ICHKB .AND. LOC.LE.ICHK) LFLAG = 2
IF (ICHKB.GE.ICHK .AND. (LOC.GE.ICHKB .OR. LOC.LE.ICHK)) LFLAG = 2
IF (LFLAG.EQ.1) GO TO 50
140 CONTINUE
LCHK = MIN0(LCHK+1,NUMCHK)
ICHK = ICHK+1
IF (ICHK.GT.NUMCHK) ICHK = 1
XCHK(ICHK) = X
YCHK(ICHK) = Y
I = IFIX(X)
J = IFIX(Y)
CALL SBYTES( UX(I,J) , IONE , ISKIP1 , 1 , 0 , 1 )
IF (NBOX.LT.5) GO TO 150
ICHKB = ICHKB+1
IF (ICHKB.GT.NUMCHK) ICHKB = 1
150 CONTINUE
GO TO 110
C
160 CONTINUE
C
C THIS SECTION DRAWS A DIRECTIONAL ARROW BASED ON THE MOST RECENT DIS-
C . PLACEMENT COMPONENTS ,DELX AND DELY, RETURNED BY GNEWPT. IN EARLIE
C . VERSIONS THIS WAS A SEPARATE SUBROUTINE (CALLED DRWDAR). IN THAT
C . CASE ,HOWEVER, FX AND FY WERE DEFINED EXTERNAL SINCE THESE
C . FUNCTIONS WERE USED BY BOTH DRWSTR AND DRWDAR. IN ORDER TO
C . MAKE ALL DEFAULT TRANSFORMATIONS STATEMENT FUNCTIONS I HAVE
C . PUT DRWDAR HERE AND I WILL USE MFLAG TO RETURN TO THE CORRECT
C . LOCATION IN THE CODE.
C
C IF ( (DELX.EQ.0.) .AND. (DELY.EQ.0.) ) GO TO 50
IF((ABS(DELX).LE.1.E-15) .AND. (ABS(DELY).LE.1.E-15)) GO TO 50
C
CALL SBYTES( UX(I,J) ,IONE , ISKIP , 1 ,0 , 1 )
D = ATAN2(-DELX,DELY)
D30 = D+0.5
170 YY = -AROWL*COS(D30)+Y
XX = +AROWL*SIN(D30)+X
CALL FL2INT (FX(XX,YY),FY(XX,YY),IFXX,IFYY)
CALL PLOTIT (IFXX,IFYY,1)
CALL FL2INT (FX(X,Y),FY(X,Y),IFX,IFY)
CALL PLOTIT (IFX,IFY,0)
IF (D30.LT.D) GO TO 180
D30 = D-0.5
GO TO 170
180 IF (MFLAG.EQ.1) GO TO 110
IF (MFLAG.EQ.2) GO TO 130
C
190 CONTINUE
C
C FLUSH PLOTIT BUFFER
C
CALL PLOTIT(0,0,0)
RETURN
END
SUBROUTINE GNEWPT (UX,VY,IMAX,JPTSY)
C
C INTERPOLATION ROUTINE TO CALCULATE THE DISPLACEMANT COMPONENTS
C . THE PHILOSPHY HERE IS TO UTILIZE AS MANY POINTS AS POSSIBLE
C . (WITHIN REASON) IN ORDER TO OBTAIN A PLEASING AND ACCURATE PLOT.
C . INTERPOLATION SCHEMES DESIRED BY OTHER USERS MAY EASILY BE
C . SUBSTITUTED IF DESIRED.
C
DIMENSION UX(IMAX,JPTSY) ,VY(IMAX,JPTSY)
COMMON /STR01/ IS ,IEND ,JS ,JEND
1 , IEND1 ,JEND1 ,I ,J
2 , X ,Y ,DELX ,DELY
3 , ICYC1 ,IMSG1 ,IGFL1
CIBM8 common /str03/ arowl,uvmsg,displ,dispc,cstop,
CIBM81 inita,initb,iterp,iterc,igflg,imsg,icyc
COMMON /STR03/ INITA , INITB , AROWL , ITERP , ITERC , IGFLG
1 , IMSG , UVMSG , ICYC , DISPL , DISPC , CSTOP
C
SAVE
C
C FDLI - DOUBLE LINEAR INTERPOLATION FORMULA
C FBESL - BESSEL 16 PT INTERPOLATION FORMULA ( MOST USED FORMULA )
C FQUAD - QUADRATIC INTERPOLATION FORMULA
C
FDLI(Z,Z1,Z2,Z3,DX,DY) = (1.-DX)*((1.-DY)*Z +DY*Z1)
1 + DX *((1.-DY)*Z2+DY*Z3)
FBESL(Z,ZP1,ZP2,ZM1,DZ)=Z+DZ*(ZP1-Z+0.25*(DZ-1.)*((ZP2-ZP1-Z+ZM1)
1 +0.666667*(DZ-0.5)*(ZP2-3.*ZP1+3.*Z-ZM1)))
FQUAD(Z,ZP1,ZM1,DZ)=Z+0.5*DZ*(ZP1-ZM1+DZ*(ZP1-2.*Z+ZM1))
C
DX = X-AINT(X)
DY = Y-AINT(Y)
C
IF( IMSG.NE.0.OR.IGFLG.NE.0) GO TO 20
C
IM1 = I-1
IP2 = I+2
C
C DETERMINE WHICH INTERPOLATION FORMULA TO USE DEPENDING ON I,J LOCATION
C . THE FIRST CHECK IS FOR I,J IN THE GRID INTERIOR.
C
IF (J.GT.JS .AND. J.LT.JEND1 .AND. I.GT.IS .AND. I.LT.IEND1)
1 GO TO 30
IF (J.EQ.JEND1 .AND. I.GT.IS .AND. I.LT.IEND1) GO TO 40
IF (J.EQ.JS) GO TO 20
C
IF (ICYC1.EQ.1) GO TO 10
C
C MUST NOT BE CYCLIC
C
IF (I.EQ.IS) GO TO 20
IF (I.EQ.IEND1) GO TO 50
GO TO 20
10 CONTINUE
C
C MUST BE CYCLIC IN THE X DIRECTION
C
IF (I.EQ.IS .AND. J.LT.JEND1) GO TO 12
IF (I.EQ.IEND1 .AND. J.LT.JEND1) GO TO 14
IF (J.EQ.JEND1 .AND. I.EQ.IS) GO TO 16
IF (J.EQ.JEND1 .AND. I.EQ.IEND1) GO TO 18
GO TO 20
12 IM1 = IEND1
GO TO 30
14 IP2 = IS+1
GO TO 30
16 IM1 = IEND1
GO TO 40
18 IP2 = IS+1
GO TO 40
C
20 CONTINUE
C
C DOUBLE LINEAR INTERPOLATION FORMULA. THIS SCHEME WORKS AT ALL POINTS
C . BUT THE RESULTING STREAMLINES ARE NOT AS PLEASING AS THOSE DRAWN
C . BY FBESL OR FQUAD. CURRENTLY THIS IS USED AT THIS IS UTILIZED
C . ONLY AT CERTAIN BOUNDARY POINTS OR IF IGFLG IS NOT EQUAL TO ZERO.
C
DELX = FDLI (UX(I,J),UX(I,J+1),UX(I+1,J),UX(I+1,J+1),DX,DY)
DELY = FDLI (VY(I,J),VY(I,J+1),VY(I+1,J),VY(I+1,J+1),DX,DY)
RETURN
30 CONTINUE
C
C USE A 16 POINT BESSEL INTERPOLATION SCHEME
C
UJM1 = FBESL (UX(I,J-1),UX(I+1,J-1),UX(IP2,J-1),UX(IM1,J-1),DX)
UJ = FBESL (UX(I,J),UX(I+1,J),UX(IP2,J),UX(IM1,J),DX)
UJP1 = FBESL (UX(I,J+1),UX(I+1,J+1),UX(IP2,J+1),UX(IM1,J+1),DX)
UJP2 = FBESL (UX(I,J+2),UX(I+1,J+2),UX(IP2,J+2),UX(IM1,J+2),DX)
DELX = FBESL (UJ,UJP1,UJP2,UJM1,DY)
VJM1 = FBESL (VY(I,J-1),VY(I+1,J-1),VY(IP2,J-1),VY(IM1,J-1),DX)
VJ = FBESL (VY(I,J),VY(I+1,J),VY(IP2,J),VY(IM1,J),DX)
VJP1 = FBESL (VY(I,J+1),VY(I+1,J+1),VY(IP2,J+1),VY(IM1,J+1),DX)
VJP2 = FBESL (VY(I,J+2),VY(I+1,J+2),VY(IP2,J+2),VY(IM1,J+2),DX)
DELY = FBESL (VJ,VJP1,VJP2,VJM1,DY)
RETURN
40 CONTINUE
C
C 12 POINT INTERPOLATION SCHEME APPLICABLE TO ONE ROW FROM TOP BOUNDARY
C
UJM1 = FBESL (UX(I,J-1),UX(I+1,J-1),UX(IP2,J-1),UX(IM1,J-1),DX)
UJ = FBESL (UX(I,J),UX(I+1,J),UX(IP2,J),UX(IM1,J),DX)
UJP1 = FBESL (UX(I,J+1),UX(I+1,J+1),UX(IP2,J+1),UX(IM1,J+1),DX)
DELX = FQUAD (UJ,UJP1,UJM1,DY)
VJM1 = FBESL (VY(I,J-1),VY(I+1,J-1),VY(IP2,J-1),VY(IM1,J-1),DX)
VJ = FBESL (VY(I,J),VY(I+1,J),VY(IP2,J),VY(IM1,J),DX)
VJP1 = FBESL (VY(I,J+1),VY(I+1,J+1),VY(IP2,J+1),VY(IM1,J+1),DX)
DELY = FQUAD (VJ,VJP1,VJM1,DY)
RETURN
50 CONTINUE
C
C 9 POINT INTERPOLATION SCHEME FOR USE IN THE NON-CYCLIC CASE
C . AT I=IEND1 ; JS.LT.J AND J.LE.JEND1
C
UJP1 = FQUAD (UX(I,J+1),UX(I+1,J+1),UX(IM1,J+1),DX)
UJ = FQUAD (UX(I,J),UX(I+1,J),UX(IM1,J),DX)
UJM1 = FQUAD (UX(I,J-1),UX(I+1,J-1),UX(IM1,J-1),DX)
DELX = FQUAD (UJ,UJP1,UJM1,DY)
VJP1 = FQUAD (VY(I,J+1),VY(I+1,J+1),VY(IM1,J+1),DX)
VJ = FQUAD (VY(I,J),VY(I+1,J),VY(IM1,J),DX)
VJM1 = FQUAD (VY(I,J-1),VY(I+1,J-1),VY(IM1,J-1),DX)
DELY = FQUAD (VJ,VJP1,VJM1,DY)
RETURN
END
SUBROUTINE CHKCYC (U,V,IMAX,JPTSY,IER)
C
C CHECK FOR CYCLIC CONDITION
C
DIMENSION U(IMAX,JPTSY) ,V(IMAX,JPTSY)
COMMON /STR01/ IS ,IEND ,JS ,JEND
1 , IEND1 ,JEND1 ,I ,J
2 , X ,Y ,DELX ,DELY
3 , ICYC1 ,IMSG1 ,IGFL1
C
SAVE
DO 10 J=JS,JEND
IF (U(IS,J).NE.U(IEND,J)) GO TO 20
IF (V(IS,J).NE.V(IEND,J)) GO TO 20
10 CONTINUE
C
C MUST BE CYCLIC
C
RETURN
20 CONTINUE
C
C MUST NOT BE CYCLIC
C . CHANGE THE PARAMETER AND SET IER = -1
C
ICYC1 = 0
IER = -1
RETURN
C
C-----------------------------------------------------------------------
C REVISION HISTORY
C
C OCTOBER, 1979 FIRST ADDED TO ULIB
C
C OCTOBER, 1980 ADDED BUGS SECTION
C
C JUNE, 1984 REMOVED STATEMENT FUNCTIONS ANDF AND ORF,
C CONVERTED TO FORTRAN77 AND GKS.
C
C MAY, 1988 CHANGED CODE (IN SUBROUTINE DRWSTR) WHICH PROTECTS
C UX ELEMENTS FROM BECOMING ZERO. THE ORIGINAL CODE
C CAUSED UNDERFLOW ON IBM MACHINES. (DJK)
C-----------------------------------------------------------------------
C
END