1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
|
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
|
