SUBROUTINE STRMLN (U,V,WORK,IMAX,IPTSX,JPTSY,NSET,IER) C C C +-----------------------------------------------------------------+ C | | C | Copyright (C) 1986 by UCAR | C | University Corporation for Atmospheric Research | C | All Rights Reserved | C | | C | NCARGRAPHICS Version 1.00 | C | | C +-----------------------------------------------------------------+ 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 LATEST REVISION JUNE 1984 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 (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 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 HISTORY WRITTEN AND STANDARDIZED IN NOVEMBER 1973. C I/O DRAWS STREAMLINES C C PRECISION SINGLE C C LANGUAGE FORTRAN C C HISTORY WRITTEN IN 1979. C CONVERTED TO FORTRAN 77 AND GKS IN JUNE 1984. 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.) CAUTION: WHEN C THIS 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 GENERALLY 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. NOTE: 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 1.E-50 THE SMALLEST REAL NUMBER C ON THE HOST COMPUTER. THIS C IS SET AUTOMATICALLY BY C R1MACH. 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 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 INITA = 2 INITB = 2 AROWL = 0.33 ITERP = 35 ITERC = -99 IGFLG = 0 ICYC = 0 IMSG = 0 C +NOAO C UVMSG = 1.E+36 uvmsg = 1.E+16 C -NOAO DISPL = 0.33 DISPC = 0.67 CSTOP = 0.50 C C THE FOLLOWING CALL IS FOR MONITORING LIBRARY USE AT NCAR 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 COMMON /STR03/ INITA , INITB , AROWL , ITERP , ITERC , IGFLG 1 , IMSG , UVMSG , ICYC , DISPL , DISPC , CSTOP COMMON /STR04/ XCHK(NCHK) ,YCHK(NCHK) , NUMCHK , UXSML 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 FY(X,Y) = Y FU(X,Y) = X FV(X,Y) = Y C C INITIALIZE C ISKIP = I1MACH(5) - 2 ISKIP1 = ISKIP + 1 UXSML = R1MACH(1) 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 IF (U(I,J).EQ.0. .AND. V(I,J).EQ.0.) GO TO 10 UX(I,J) = FU(U(I,J),V(I,J)) VY(I,J) = FV(U(I,J),V(I,J)) 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) C IF(UX(I,J) .EQ. 0.) UX(I,J) = CON*FU(UXSML,V(I,J)) C GO TO 20 10 CONTINUE 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 SOME VERY SMALL NUMBER. C UX(I,J) = FU(UXSML,0.) VY(I,J) = 0. 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 20 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 IF ( (DELX.EQ.0.) .AND. (DELY.EQ.0.) ) 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 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 EZSTRM(U,V,WORK,IMAX,JMAX) C DIMENSION U(IMAX,JMAX) ,V(IMAX,JMAX) ,WORK(1) C SAVE C C THE FOLLOWING CALL IS FOR MONITORING LIBRARY USE AT NCAR C CALL Q8QST4 ( 'GRAPHX', 'STRMLN', 'EZSTRM', 'VERSION 01') C CALL STRMLN(U,V,WORK,IMAX,IMAX,JMAX,0,IER) 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------------------------------------------------------------------- END