SUBROUTINE SUPMAP (JPROJ,POLAT,POLON,RROT,PL1,PL2,PL3,PL4,JJLTS, + JGRID,JUS,JDOT,IER) C SupMap is used to plot map outlines according to one of nine projections. C The origin and orientation of the projection are selected by the user. C Points on the Earth defined by latitude and longitude are transformed to C points in the u,v plane, the plane of projection. The u and v axes are C respectively parallel to the x and y axes of the plotter. A rectangular C frame parallel to the u and v axes is chosen and only material within the C frame is plotted. C Continental and U.S. state and county outlines are available for plotting. C In addition, a western U.S. subset of the counties has been extracted, and C up to three user files may be read. C Name Date Description C -------------- ----------- --------------------------------------- C ?? (NCAR) JAN 69 "Revised" C MAY 71 "Revised" C OCT 73 "Standardized" C JUL 74 "Revised" C AUG 76 "Revised" C M. Cairns 27 Jan 81 Added MapCol common for intensities C J. Wakefield 19 Jun 81 Added MapDas common for dash patterns C 29 Jun 81 Removed trailing blanks and modified C to read unformatted files. C 28 Sep 81 Added common block SupMpB. C 23 Mar 82 Added negative JDot function, to allow C suppression of SupMap call printing if C JUS=0. See comments below. C 16 Feb 83 Added SysDisk to data file names and C added error code for open failure. C 10 Apr 84 Added user files option. C 14 Jan 85 Fixed bug in Part < 0 option. C J. Ramer 13 Nov 85 Added common blocks necessary for using C MapInv_GC with all projections and for C recreating a map from the contents of C the common blocks. Added subroutine C RESUP which recreates the map. C J. Wakefield 19 Nov 86 Chgd SysDisk to Lib_Dev for files. C Usage C Call SupMap (JProj,PoLat,PoLong,RRot,PL1,PL2,PL3,PL4, C JJLTS,JGrid,JUS,JDot,IER) C Dimension of PL1(2),PL2(2),PL3(2),PL4(2) C arguments C On input JProj C for SupMap |JProj| defines the projection type C according to the following code: C 1 Stereographic C 2 Orthographic C 3 Lambert Conformal Conic with two standard C parallels C 4 Lambert Equal Area C 5 Gnomonic C 6 Azimuthal Equidistant C 7 Dummy -- this code is not used C 8 Cylindrical Equidistant C 9 Mercator C 10 Mollweide type C If JProj < 0, the map and grid lines are omitted. C PoLat,PoLong,RRot C If (|JProj|.ne.3) C . PoLat and PoLong define in degrees the latitude and C longitude of the point on the globe which is to C transform to the origin of the u,v plane. C -90 .le. PoLat .le. 90 C -180 .le. PoLong .le. 180 C degrees of latitude north of the Equator and C degrees of longitude east of the Greenwich Meridian C are positive. If the origin is at the North Pole, C "north" is considered to be in the direction of C (PoLong+180.). If the origin is at the South Pole, C "north" is in the direction of PoLong. C . RRot is the angle between the v axis and north at C the origin. It is measured in degrees and is taken C to be positive if the angular movement from north C to the v axis is counter-clockwise. For the C cylindrical projections (8,9,10), the axis of the C projection is parallel to the v axis. C If (|JProj|.eq.3) (Lambert Conformal Conic) C . PoLong = central meridian of projection in degrees. C . PoLat,RRot are the two standard parallels in deg. C JJLTS,PL1,PL2,PL3,PL4 C |JJLTS| can take the values 1 through 5 and specifies C one of five options on the way in which the limits of C the rectangular map are defined by the parameters C PL1, PL2, PL3, and PL4. C |JJLTS| = 1 C The maximum useful area produced by the projection C is plotted. PL1, PL2, PL3, and PL4 are not used C and may be set to zero. C |JJLTS| = 2 C In this case (PL1,PL2) and (PL3,PL4) are the C latitudes and longitudes in degrees of two points C which are to be at opposite corners of the map C (upper right and lower left, respectively). C Care must be taken when using cylindrical C projections and this option. C |JJLTS| = 3 C The minimum and maximum values of u and v are C specified by PL1 through PL4. PL1 = UMin, C PL2 = UMax, PL3 = VMin, PL4 = VMax. Knowledge of C the transformation equations is necessary for this C option to be used (see below). C |JJLTS| = 4 C Here PL1 = AUMin, PL2 = AUMax, PL3 = AVMin, C PL4 = AVMax, where C AUMin = Angular distance from origin to left C frame of map. C AUMax = Angular distance from origin to right C frame of map. C AVMin = Angular distance from origin to lower C frame. C AVMax = Angular distance from origin to upper C frame. C AUMin, AUMax, AVMin, AVMax must be positive and the C origin must be within the rectangular limits of the C map. This option is useful for polar projections. C It is not appropriate for the Lambert Conformal C with two standard parallels. An error message is C printed if an attempt is made to use JJLTS = 4 when C JProj = 3, (see below). C |JJLTS| = 5 C PL1 through PL4 are two element arrays giving the C latitudes and longitudes of four points which are C to be on the four sides of the rectangular frame. C PL1(1), PL1(2) are respectively the latitude and C longitude of a point on the left frame. Similarly, C PL2 lies on the right frame, PL3 lies on the lower C frame and PL4 lies on the upper frame. Note that C in the calling program PL1 through PL4 will be C dimensioned: C Real*4 PL1(2),PL2(2),PL3(2),PL4(2) C .If JJLTS is positive, the SupMap call is plotted C below the map. This is omitted if JJLTS is < 0. C JGrid C |JGrid| gives in degrees the interval at which lines C of latitude and longitude are to be plotted. A value C in the range 1 through 10 will usually be appropriate C but higher values are acceptable. C If JGrid < 0 the border around the map is omitted. C If JGrid = 0 no grid lines are plotted. C JUS C |JUS| C 1 World outlines C 2 U.S. state outlines C 4 U.S. counties and states C 8 PROFS-region counties and states (western U.S.) C 16,32,64 SupMap user file(s) 1,2,3 C To access a combination of datasets, use the sum of C their values (e.g. 10 = U.S. + PROFS area counties). C .If JUS is positive, the SupMap call and values of C of UMin, UMax, VMin, VMax are printed as an aid to C debugging. This is omitted if JUS is negative. C JDot C |JDot| C 0 for continuous outlines. C 1 for dotted outlines. C .If JDot is negative, the SupMap call is neither C printed nor plotted. C On output All arguments except IER are unchanged. C for SupMap C IER C Error flag with the following meanings: C If IER = C 0 Map successfully plotted. C 29 Error opening data file. C 33 Attempt to use non-existent projection. C 34 Map limits inappropriate. C 35 Angular limits too great. C 36 Map has zero area. C Entry points MaPlot, SupCon, SupFst, SupMap, SupTrp, SupVec, QCon, C QVec, VecPlt, ReSup C MaPlot C Actually draws the map. C SupCon C Once the transformation has been set up by an initial C call to SupMap, the subroutine SupCon may be called C to transform a point, (latitude, longitude) to the C corresponding point, (u, v) on the plane. Contours C may thus be readily drawn against the map background. C (See SupFst and SupVec below). C Call SupCon(RLat,RLon,U,V) C On input: C RLat,RLon are the latitude and longitude of a point C to be transformed to the u,v plane. C -90. .le. RLat .le. 90. C -180. .le. RLon .le. 180. C On output: C RLat,RLon are unchanged. C U,V are the transformed coordinates of the point. C QCon C Actually performs the above mentioned transformation. C SupFst C SupVec C To facilitate drawing lines on the map these routines C which act like the plotting routines FrstPt and C Vector are included. They are subject to the same C restrictions as SupCon above. C Call SupFst (RLat,RLon) C Call SupVec (RLat,RLon) C QVec C Decides what lines are to be drawn and where. C SupTrp C Performs interpolation to the edges of the frame. C VecPlt C Called by QVec to draw (dot) lines on the plotter. C Common blocks Name Length C SupMp1 9 FP C SupMp2 1 Int + 204 FP C SupMp3 2 Int + 5 FP C SupMp4 5 Int + 2 FP C SupMp5 1 Int + 5 FP C SupMp6 6 FP C SupMp7 3 Int + 2 FP C SupMp8 6 FP C SupMp9 3 FP C SupMpA 1 Int C SupMpB 4 FP C SupMpC 4 Int + 1 FP C SupMpD 9 FP C MapCol 4 Int C MapDas 4 Int C I/O Map plotted. Outline data is read from any of several C disk files. SupMap call printed. C Precision Single C Language FORTRAN C Algorithm The latitudes and longitudes of successive outline C points are transformed to coordinates in the plane of C projection and joined by a vector. C References Hershey, A. V., The plotting of maps on a CRT printer. C NWL Report No. 1844, 1963. C Lee, Tso-Hwa, Students Summary Reports, Work-study C program in scientific computing. NCAR 1968. C Parker, R. L., UCSD SuperMap: World Plotting Package. C Steers, J.A., An Introduction to the Study of Map C Projections. Univ. of London Press, 1962. C Accuracy The definition of the map produced is limited by the C fact that the resolution of the virtual plotter space C is 1024 units in the x and y directions. C Plotting routines PWRT, FrstPt, Vector, Point, DashLn, Perim, Set C used C Required resident ATan, Tan, Sin, Cos, ALog, SqRt, ATan2, ACos C routines COMMON/SUPMP1/PI,TOVPI,DTR,RTD,EPS,OV90,CON1,CON2,PART COMMON/SUPMP3/POLONG,CONE,RLAT,RLON,JGR,ILF,SGN COMMON/SUPMP4/IFST,IGO,IGOLD,ICROSS,IOUT,UOLD,VOLD COMMON/SUPMP5/PHIOC,SINO,COSO,SINR,COSR,IPROJ COMMON/SUPMP6/UMIN,UMAX,VMIN,VMAX,UEPS,VEPS COMMON/SUPMP7/PHIO,PHIA,IGRID,IDOT,ILTS COMMON/SUPMP8/U,V,U1,V1,U2,V2 COMMON/SUPMP9/DS,DI,DSRDI COMMON/SUPMPA/IIER COMMON/SUPMPB/X1,Y1,X2,Y2 COMMON/SUPMPC/LPROJ,ROT,JLTS,LGRID,IUS COMMON/SUPMPD/AAA,BBB,CCC,DDD,EEE,FFF,GGG,HHH,III DIMENSION PL1(2),PL2(2),PL3(2),PL4(2),LABA(20),LABB(18) REAL LAT1,LAT2,AAA,BBB,CCC,DDD,EEE,FFF,GGG,HHH,III EQUIVALENCE (PLA1,AUMIN),(PLA2,AUMAX),(PLA3,AVMIN),(PLA4,AVMAX), 1 (PHIA,LAT1),(ROT,LAT2) DATA PLTRES,RESLIM/ 1024.,10./ EXTERNAL SUPMBD SQU(X) = X*X IGDFLT = %LOC(SUPMBD) !FORCE LOAD OF BLOCK DATA from library ROT = RROT ILTS = IABS(JJLTS) JLTS = JJLTS LGRID = JGRID IGRID = IABS(LGRID) c JGR = ISIGN(1,LGRID) JGR = LGRID IUS = JUS IUSGN = ISIGN(1,IUS) PLA1 = PL1(1) PLA2 = PL2(1) PLA3 = PL3(1) PLA4 = PL4(1) LABL = 77 IDOT = IABS(JDOT) C INITIALIZATION LPROJ = JPROJ IPROJ = IABS(LPROJ) JPR = ISIGN(1,LPROJ) C JGR=JPR IOUT = ABS(IUS) !USE IOUT FOR CHOOSING DATA FILE(S) PHIA = POLAT POLONG = POLON PHIO = POLON PHIOC = 540.-PHIO ICROSS = 0 IIER = 0 ILF = 0 C COMPUTE CONSTANTS APPROPRIATE TO EACH PROJECTION IF (IPROJ .NE. 3) GO TO 30 C LAMBERT CONFORMAL CONIC SGN = SIGN(1.,0.5*(LAT1+LAT2)) CHI1 = (90.-SGN*LAT1)*DTR IF (LAT1 .EQ. LAT2) GO TO 20 CHI2 = (90.-SGN*LAT2)*DTR CONE = ALOG(SIN(CHI1)/SIN(CHI2))/ALOG(TAN(0.5*CHI1)/TAN(0.5*CHI2)) GO TO 60 20 CONE = COS(CHI1) GO TO 60 C THE OTHERS 30 XT = ROT*DTR SINR = SIN(XT) COSR = COS(XT) XT = PHIA*DTR SINO = SIN(XT) COSO = COS(XT) C CALCULATE COEFFICIENTS NECESSARY TO CONVERT DISPLACEMENT VECTOR ON EARTH'S C SURFACE FROM CARTESIAN SYSTEM WITH ORIGIN AT CENTER OF EARTH, UNIT VECTOR K C K AT THE PROJECTION CENTER, AND UNIT VECTOR J IN THE POSITIVE U DIRECTION TO C A SYSTEM WITH UNIT VECTOR K AT NORTH POLE AND UNIT VECTOR I AT INTERSECTION C OF EQUATOR AND GREENWICH MERIDIAN. NEEDED LATER FOR MAPINV_GC. AAA=COS(ROT*DTR)*COS(POLONG*DTR)*SIN(PHIA*DTR)- & SIN(ROT*DTR)*SIN(POLONG*DTR) BBB=COS(ROT*DTR)*SIN(POLONG*DTR)*SIN(PHIA*DTR)+ & SIN(ROT*DTR)*COS(POLONG*DTR) CCC=-COS(ROT*DTR)*COS(PHIA*DTR) DDD=-SIN(ROT*DTR)*COS(POLONG*DTR)*SIN(PHIA*DTR)- & COS(ROT*DTR)*SIN(POLONG*DTR) EEE=-SIN(ROT*DTR)*SIN(POLONG*DTR)*SIN(PHIA*DTR)+ & COS(ROT*DTR)*COS(POLONG*DTR) FFF=SIN(ROT*DTR)*COS(PHIA*DTR) GGG=COS(POLONG*DTR)*COS(PHIA*DTR) HHH=SIN(POLONG*DTR)*COS(PHIA*DTR) III=SIN(PHIA*DTR) IF (IPROJ-7) 60,67,40 C CYLINDRICAL PROJECTIONS [8,9,10] 40 IF (PHIA .NE. 0.0) GO TO 42 IF (ROT .EQ. 0.0) GO TO 45 IF (ABS(ROT) .EQ. 180.) GO TO 50 42 SINO1 = COSO*COSR COSO1 = SQRT(CON1-SINO1*SINO1) OVC1 = 1./COSO1 PHIO = PHIO-ATAN2(SINR*OVC1,-COSR*SINO*OVC1)*RTD PHIOC = 540.-PHIO SINR = SINR*COSO*OVC1 COSR = -SINO*OVC1 SINO = SINO1 COSO = COSO1 GO TO 60 C USE SIMPLE TRANSFORMS FOR CYLINDRICAL PROJECTIONS IF ROT = POLAT = .0 C I.E. IPROJ = 11, 12, 13 45 SINO = 1.0 IPROJ = IPROJ+3 GO TO 55 50 SINO = -1.0 PHIO = PHIO+180. PHIOC = PHIOC+180. 55 COSO = 0.0 SINR = 0.0 COSR = 1.0 ILF = 1 C ILTS = 1 THE MAXIMUM USEFUL AREA IS PLOTTED. C --------- 60 GO TO (61 ,62 ,62 ,61 ,61 ,66 ,67 ,68 ,66 ,70 , 1 68 ,66 ,70 ),IPROJ C STEREOGRAPHIC [1] C LAMBERT EQUAL AREA [4] C GNOMONIC [5] 61 UMIN = -2.0 VMIN = -2.0 UMAX = 2.0 VMAX = 2.0 GO TO 80 C ORTHOGRAPHIC [2] C LAMBERT CONFORMAL CONIC [3] 62 UMIN = -1.0 VMIN = -1.0 UMAX = 1.0 VMAX = 1.0 GO TO 80 C AZIMUTHAL EQUIDISTANT [6] C MERCATOR WITH ARBITRARY POLE [9] C MERCATOR [12] 66 UMAX = PI VMAX = PI UMIN = -PI VMIN = -PI GO TO 80 C DUMMY -- ERROR EXIT [7] 67 IIER = 33 CALL ULIBER (IIER, 1 46H SUPMAP-ATTEMPT TO USE NON-EXISTENT PROJECTION,46) GO TO 700 C CYLINDRICAL EQUIDISTANT [8,11] 68 UMAX = 180. UMIN = -180. VMAX = 90. VMIN = -90. GO TO 80 C MOLLWEIDE TYPE [10,13] 70 UMAX = 2.0 UMIN = -2.0 VMAX = 1.0 VMIN = -1.0 80 UEPS = 0.5*(UMAX-UMIN) VEPS = 0.5*(VMAX-VMIN) IF (IPROJ .EQ. 3) UEPS = 180. C COMPUTE THE APPROPRIATE MAP BOUNDARIES. GO TO (600,200,300,400,500),ILTS C ILTS = 2 POINT (PL1,PL2) IN UPPER RIGHT CORNER , (PL3,PL4) IN C --------- LOWER LEFT CORNER OF PLOT. 200 RLAT = PLA1 RLON = PLA2 CALL QCON U1 = U V1 = V RLAT = PLA3 RLON = PLA4 CALL QCON UMAX = AMAX1(U1,U) UMIN = AMIN1(U1,U) VMAX = AMAX1(V1,V) VMIN = AMIN1(V1,V) GO TO 600 C ILTS = 3 SET PLOT LIMITS DIRECTLY. C ---------- 300 UMAX = PLA2 UMIN = PLA1 VMAX = PLA4 VMIN = PLA3 GO TO 600 C ILTS = 4 USE ANGULAR DISTANCES TO SET PLOT LIMITS. C ---------- 400 COSUMI = COS(AUMIN*DTR) SINUMI = SQRT(CON1-COSUMI*COSUMI) COSUMA = COS(AUMAX*DTR) SINUMA = SQRT(CON1-COSUMA*COSUMA) COSVMI = COS(AVMIN*DTR) SINVMI = SQRT(CON1-COSVMI*COSVMI) COSVMA = COS(AVMAX*DTR) SINVMA = SQRT(CON1-COSVMA*COSVMA) GO TO (401,402,403,404,405,406,407,408,409,410, 1 408,409,410),IPROJ C STEREOGRAPHIC [1] 401 UMAX = (1.-COSUMA)/SINUMA UMIN = -(1.-COSUMI)/SINUMI VMAX = (1.-COSVMA)/SINVMA VMIN = -(1.-COSVMI)/SINVMI GO TO 600 C ORTHOGRAPHIC [2] 402 IF (AMAX1(AUMIN,AUMAX,AVMIN,AVMAX) .GT. 90.) GO TO 900 UMAX = SINUMA UMIN = -SINUMI VMAX = SINVMA VMIN = -SINVMI GO TO 600 C LAMBERT CONFORMAL CONIC [3] 403 IIER = 34 CALL ULIBER (IIER,32H SUPMAP-MAP LIMITS INAPPROPRIATE,32) GO TO 700 C LAMBERT EQUAL AREA [4] 404 UMAX = (1.+COSUMA)/SINUMA UMIN = (1.+COSUMI)/SINUMI VMAX = (1.+COSVMA)/SINVMA VMIN = (1.+COSVMI)/SINVMI UMAX = 2./SQRT(1.+UMAX*UMAX) UMIN = -2./SQRT(1.+UMIN*UMIN) VMAX = 2./SQRT(1.+VMAX*VMAX) VMIN = -2./SQRT(1.+VMIN*VMIN) GO TO 600 C GNOMONIC [5] 405 IF (AMAX1(AUMIN,AUMAX,AVMIN,AVMAX) .GE. 90.) GO TO 900 UMAX = SINUMA/COSUMA UMIN = -SINUMI/COSUMI VMAX = SINVMA/COSVMA VMIN = -SINVMI/COSVMI GO TO 600 C AZIMUTHAL EQUIDISTANT [6] 406 UMAX = AUMAX*DTR UMIN = -AUMIN*DTR VMAX = AVMAX*DTR VMIN = -AVMIN*DTR GO TO 600 C DUMMY -- ERROR EXIT [7] 407 GO TO 67 C CYLINDRICAL EQUIDISTANT [8,11] 408 UMAX = AUMAX UMIN = -AUMIN VMAX = AVMAX VMIN = -AVMIN GO TO 600 C MERCATOR [9,12] 409 IF (AMAX1(AVMIN,AVMAX) .GE. 90.) GO TO 900 UMAX = AUMAX*DTR UMIN = -AUMIN*DTR VMAX = ALOG((1.+SINVMA)/COSVMA) VMIN = -ALOG((1.+SINVMI)/COSVMI) GO TO 600 C MOLLWEIDE TYPE [10,13] 410 UMAX = AUMAX*OV90 UMIN = -AUMIN*OV90 VMAX = SINVMA VMIN = -SINVMI GO TO 600 C ILTS = 5 USE FOUR EDGE POINTS TO SET LIMITS. C ---------- 500 PLB1 = PL1(2) RLAT = PLA1 RLON = PLB1+EPS CALL QCON UMIN = U PLB2 = PL2(2) RLAT = PLA2 RLON = PLB2-EPS CALL QCON UMAX = U PLB3 = PL3(2) RLAT = PLA3 RLON = PLB3 CALL QCON VMIN = V PLB4 = PL4(2) RLAT = PLA4 RLON = PLB4 CALL QCON VMAX = V C COMPUTE MAP LIMITS FOR PLOT 600 DU = UMAX-UMIN DV = VMAX-VMIN C ERROR IF MAP HAS ZERO AREA IF (DU.EQ.0.0 .OR. DV.EQ.0.0) GO TO 905 IF(PART.LE..0)GOTO 620 !USER-SUPPLIED X1,Y1,X2,Y2 IF (DU .GT. DV) GO TO 610 Y1 = 0.5*(1.-PART) Y2 = 1.-Y1 X1 = 0.5*(1.-PART*DU/DV) X2 = 1.-X1 GO TO 620 610 X1 = 0.5*(1.-PART) X2 = 1.-X1 Y1 = 0.5*(1.-PART*DV/DU) Y2 = 1.-Y1 C ERROR IF MAP HAS ESSENTIALLY ZERO AREA 620 IF (AMIN1(X2-X1,Y2-Y1)*PLTRES .LT. RESLIM) GO TO 905 CALL SET (X1,X2,Y1,Y2,UMIN,UMAX,VMIN,VMAX,1) DS = SQU(((X2-X1)*PLTRES)/DU) DSRDI = SQRT(DS/DI) C DO WE WRITE ANYTHING? IF (JLTS.LT.0 .AND. IUSGN.LT.0) GO TO 640 IF(JDOT.LT.0)GOTO 640 C CREATE THE LABEL ENCODE (LABL,7000,LABA(1)) LPROJ,PHIA,POLONG,ROT,PLA1,PLA2,PLA3, 1 PLA4,JLTS,LGRID,IUS,IDOT 7000 FORMAT (8H SUPMAP(,I3,7(1H,,F6.1),4(1H,,I3),1H)) IF (ILTS .EQ. 5) ENCODE (61,7010,LABB(1)) PLB1,PLB2,PLB3,PLB4 7010 FORMAT (7X,1H(,24X,4(1H,,F6.1)1H)) IF (JLTS .LT. 0) GO TO 630 C WRITE SUPMAP CALL BENEATH THE MAP CALL PWRT (240,17,LABA(1),LABL,0,0) IF (ILTS .EQ. 5) CALL PWRT (240,1,LABB(1),61,0,0) 630 IF (IUSGN .LT. 0) GO TO 640 C PRINT OUT THE CALL ET AL. K = (LABL+3)/4 WRITE (6,6000) (LABA(I),I=1,K) 6000 FORMAT (30A4) IF (ILTS .EQ. 5) WRITE (6,6000) (LABB(I),I=1,18) WRITE (6,6010) UMIN,UMAX,VMIN,VMAX 6010 FORMAT (8H UMIN = ,F11.6,9H UMAX = ,F11.6,9H VMIN = , + F11.6,9H VMAX = ,F11.6) C DRAW THE MAP 640 IF (IOUT.NE.0 .OR. IGRID.NE.0 .OR. JGR.GE.0) CALL MAPLOT IDOT = 0 C RETURN IER 700 IER = IIER RETURN C ERROR RETURNS 900 IIER = 35 CALL ULIBER (IIER,32H SUPMAP-ANGULAR LIMITS TOO GREAT,32) GO TO 700 905 IIER = 36 CALL ULIBER (IIER,25H SUPMAP-MAP HAS ZERO AREA,25) GO TO 700 END C------------------------------------------------------------------------------- SUBROUTINE MAPLOT C THIS SUBROUTINE PLOTS THE CONTINENTAL AND U.S. STATE OUTLINES, C MERIDIANS, PARALLELS, LIMBS WHERE APPROPRIATE. IT LABELS KEY MERIDIANS C AND POLES, AND IT DRAWS A BORDER. COMMON/SUPMP1/PI,TOVPI,DTR,RTD,EPS,OV90,CON1,CON2,PART COMMON/SUPMP2/NPTS,MAXLAT,MINLAT,MAXLON,MINLON,PTS(200) COMMON/SUPMP3/POLONG,CONE,RLAT,RLON,JGR,ILF,SGN COMMON/SUPMP4/IFST,IGO,IGOLD,ICROSS,IOUT,UOLD,VOLD COMMON/SUPMP5/PHIOC,SINO,COSO,SINR,COSR,IPROJ COMMON/SUPMP6/UMIN,UMAX,VMIN,VMAX,UEPS,VEPS COMMON/SUPMP7/PHIO,PHIA,IGRID,IDOT,ILTS COMMON/SUPMP8/U,V,U1,V1,U2,V2 COMMON/SUPMPA/IIER COMMON/MAPCOL/MPCOL1,MPCOL2,MPCOL3,MPCOL4 COMMON/MAPDAS/LDASH1,LDASH2,LDASH3,LDASH4 DATA !Default line intensities and dash patterns 1 MPCOL1,LDash1 /255,'1777'O/, !Map lines 2 MPCOL2,LDash2 /128,'1756'O/, !Grid lines 3 MPCOL3,LDash3 /192,'1777'O/, !Limb lines 4 MPCOL4,LDash4 /255,'1777'O/ !Perimeter DIMENSION SPLAT(2) REAL MAXLAT,MINLAT,MAXLON,MINLON,MIDLAT,MIDLON CHARACTER*30 NAMFIL DATA SINLMB,COSLMB /0.017452406, 0.99984765/ DATA FLOORC / 10000. / FLOOR(X) = AINT(X+FLOORC)-FLOORC CLING(X) = FLOOR(X)+1. CALL DASHLN(LDASH1) CALL OPTN(2HIN,MPCOL1) GO TO (10,10,10,5,10,5,905,5,5,5,5,5,5),IPROJ 5 ICF=1 10 IF(IOUT.EQ.0)GOTO 100 C***Select appropriate file(s) according to bits set in IOut IBS=1 !START WITH 1ST BIT MASK=1 11 DO 12 IBIT=IBS,16 !CHECK LOWER 16 BITS IBS=IBS+1 IF((IOUT.AND.MASK).NE.0)GOTO 13 !EXIT LOOP MASK=MASK*2 12 CONTINUE GO TO 100 !FELL OUT OF LOOP SO DONE 13 MASK=MASK*2 GoTo(14,15,16,17,18,19,20)IBit GO TO 100 !OUT OF RANGE -- RETURN 14 IF((IOUT.AND.2).NE.0)THEN ! STATES WITH CONTINENTS? IBS=IBS+1 !SKIP BIT LATER MASK=MASK*2 NAMFIL='Lib_Dev:[GUDAT]CONANDSTA.DAT' !3 (SPECIAL CASE) ELSE NAMFIL='Lib_Dev:[GUDAT]CONTINENT.DAT' !1 ENDIF GO TO 25 15 NAMFIL='Lib_Dev:[GUDAT]STATE.DAT' !2 GO TO 25 16 NAMFIL='Lib_Dev:[GUDAT]USCOUNTY.DAT' !4 GO TO 25 17 NAMFIL='Lib_Dev:[GUDAT]COUNTY.DAT' !8 GO TO 25 18 NamFil='SupMap_UserFile1' !16 GoTo 25 19 NamFil='SupMap_UserFile2' !32 GoTo 25 20 NamFil='SupMap_UserFile3' !64 GoTo 25 C***Open file 25 Open(3,Name=NamFil,Type='Old',Form='UnFormatted',ReadOnly,Err=26) GoTo 30 C***Error opening file -- return 26 IIER=29 Return C***Read next line 30 Read(3,End=99)NPts,MaxLat,MinLat,MaxLon,MinLon,(Pts(M),M=1,NPts) 99 NPts=NPts/2 IF (NPTS .EQ. 0)THEN CLOSE(3) GOTO 11 !CHECK NEXT BIT ENDIF IF (NPTS .LE. 16) GO TO 70 IF (ICF .NE. 0) GO TO 70 C DOES THIS LINE INTERSECT THE SCREEN? C 1 2 3 C 4 5 C 6 7 8 MIDLAT = (MAXLAT+MINLAT)*0.5 MIDLON = (MAXLON+MINLON)*0.5 RLAT = MAXLAT RLON = MAXLON CALL QCON X3 = U Y3 = V RLON = MIDLON CALL QCON X2 = U Y2 = V RLON = MINLON CALL QCON X1 = U Y1 = V RLAT = MIDLAT CALL QCON X4 = U Y4 = V RLON = MAXLON CALL QCON X5 = U Y5 = V RLAT = MINLAT CALL QCON X8 = U Y8 = V RLON = MIDLON CALL QCON X7 = U Y7 = V RLON = MINLON CALL QCON X6 = U Y6 = V XMN = AMIN1(X1,X2,X3,X4,X5,X6,X7,X8) XMX = AMAX1(X1,X2,X3,X4,X5,X6,X7,X8) YMN = AMIN1(Y1,Y2,Y3,Y4,Y5,Y6,Y7,Y8) YMX = AMAX1(Y1,Y2,Y3,Y4,Y5,Y6,Y7,Y8) DX = AMIN1(XMX-XMN,180.) DY = AMIN1(YMX-YMN,180.) XMX = XMX+0.01*DX XMN = XMN-0.01*DX YMX = YMX+0.01*DY YMN = YMN-0.01*DY IF (XMN.GT.UMAX .OR. XMX.LT.UMIN .OR. YMN.GT.VMAX .OR. 1 YMX.LT.VMIN) GO TO 30 70 RLAT = PTS(1) RLON = PTS(2) IFST = 1 IGO = 0 CALL QVEC DO 75 J=2,NPTS RLAT = PTS(2*J-1) RLON = PTS(2*J) CALL QVEC 75 CONTINUE GO TO 30 C***DONE PLOTTING -- DRAW AND LABEL GRID, IF DESIRED 100 SPLAT(2) = 90. SPLAT(1) = -90. IF (IGRID .EQ. 0) GO TO 300 IDOTS = IDOT IDOT = 0 c if (idot .eq. 0) go to 200 C LETTER KEY MERIDIANS AND POLES C SOUTH POLE ISPF = 0 RLAT = -90. RLON = 0.0 CALL QCON IF ((U .GT. UMAX) .OR. (U .LT. UMIN) .OR. (V .GT. VMAX) .OR. 1 (V .LT. VMIN)) GO TO 110 USP = U VSP = V ISPF = 1 IPF = 1 If (phia .LT. 0) CALL PWRT (U,V,2HSP,2,1,0) C NORTH POLE 110 INPF = 0 IPF = 0 RLAT = 90. CALL QCON IF ((U .GT. UMAX) .OR. (U .LT. UMIN) .OR. (V .GT. VMAX) .OR. 1 (V .LT. VMIN)) GO TO 120 UNP = U VNP = V INPF = 1 IPF = 1 If (phia .GT. 0) CALL PWRT (U,V,2HNP,2,1,0) C EQUATOR 120 RLON = PHIO-10. RLAT = 0.0 DO 125 I=1,36 RLON = RLON+10. CALL QCON IF (U.LE.UMAX .AND. U.GE.UMIN .AND. V.LE.VMAX .AND. V.GE.VMIN) 1 GO TO 130 125 CONTINUE GO TO 140 130 CALL PWRT (U,V,2HEQ,2,1,0) C GREENWICH MERIDIAN 140 RLAT = 85. RLON = 0.0 DO 145 I=1,16 RLAT = RLAT-10. CALL QCON IF (U.LE.UMAX .AND. U.GE.UMIN .AND. V.LE.VMAX .AND. V.GE.VMIN) 1 GO TO 150 145 CONTINUE GO TO 160 150 CALL PWRT (U,V,2HGM,2,1,0) C DATE LINE 160 RLAT = 85. RLON = 180. DO 165 I=1,16 RLAT = RLAT-10. CALL QCON IF (U.LE.UMAX .AND. U.GE.UMIN .AND. V.LE.VMAX .AND. V.GE.VMIN) 1 GO TO 170 165 CONTINUE GO TO 200 170 CALL PWRT (U,V,1HI,1,1,0) 200 RGRID = IGRID c CALL DASHLN (LDASH2) c CALL OPTN(2HIN,MPCOL2) C SHOULD WE BOTHER LIMITING GRID POINTS TRANSFORMED? IF (ICF .NE. 0) GO TO 270 IF (IPROJ.GE.8 .AND. IPROJ.LE.10) GO TO 270 C SET UP TO FIND EXTREMA DLON = RGRID STLON = FLOOR(POLONG/RGRID)*RGRID IF (ISPF.NE.0 .AND. INPF.EQ.0) STLON = STLON+180. RLON = STLON-DLON SPLON = STLON+360. J = 0 PSIGN = 1. C CHECK FOR SOUTH POLE IF (ISPF .NE. 0) PSIGN = -1. C DO WE GRID POLES SPECIALLY? SPLAT(2) = 90.*PSIGN SPLAT(1) = SPLAT(2) C IF BOTH POLES WITHIN FRAME JUMP. IF (INPF.NE.0 .AND. ISPF.NE.0) GO TO 270 C IF EITHER IN FRAME USE AS BASE IF (INPF.NE.0 .OR. ISPF.NE.0) GO TO 230 C NO POLE IS CLOSE TO THE WINDOW J = -1 SPLAT(2) = FLOOR(PHIA/RGRID)*RGRID IF (ABS(SPLAT(2)) .EQ. 90.) SPLAT(2) = 0.0 C SEARCH FOR FIRST POINT WITHIN FRAME. 210 RLON = RLON+DLON DLAT = RGRID RLAT = SPLAT(2)-DLAT 215 RLAT = RLAT+DLAT CALL QCON IF ((U .LE. UMAX) .AND. (U .GE. UMIN) .AND. (V .LE. VMAX) .AND. 1 (V .GE. VMIN)) GO TO 225 IF (ABS(RLAT) .LT. 90.) GO TO 215 IF (DLAT .LT. 0.0) GO TO 220 C REVERSE LATITUDE SEARCH DIRECTION RLAT = SPLAT(2)+DLAT DLAT = -DLAT GO TO 215 C UPDATE LONGITUDE ! QUIT. 220 J = 0 IF (RLON-SPLON) 210,300,300 C SET UP FOR LIMIT SEARCH 225 J = J+1 STLON = RLON RLON = STLON-DLON IF (RLAT .EQ. 0.0) RLAT = SIGN(RLAT,-PSIGN) SPLAT(2) = RLAT SPLAT(1) = SPLAT(2) splat(1) = 90.0 splat(2) = -90.0 stlon = stlon - dlon splon = splon + dlon rlon = stlon 226 if (rlon .le. splon) then rlat2 = 80.0 if (amod(rlon,90.0) .eq. 0.0) rlat2 = 90.0 rlat = -rlat2 ifst = 1 igo = 0 call qvec 227 rlat = rlat + 1.0 if (rlat .le. rlat2) then if (iproj .eq. 1) rlat = rlat2 call qvec if (igo .ne. 0) then if (rlat .lt. splat(1)) splat(1) = rlat if (rlat .gt. splat(2)) splat(2) = rlat endif go to 227 endif rlon = rlon + dlon go to 226 endif if (dlon .ne. 0.0) go to 285 C LONGITUDE LOOP C IGF FLAG TO SIGNAL NO POINTS WITHIN WINDOW. C IPF FLAG SIGNALS WHETHER A POLE LIES WITHIN THE FRAME. C ILF FLAG SIGNALS WHETHER TO PLOT COMPLETE LONGITUDES C (I.E. TO POLE FOR ALL LATITUDES.) 230 RLON = RLON+DLON IF (RLON.GE.SPLON .OR. RLON.LT.STLON) GO TO 285 I1 = IPF I2 = MOD(I1+1,2) TSA = PSIGN DLAT = -PSIGN DX = AMOD(90.,RGRID) IF (DX .EQ. 0.0) DX = RGRID XLAT = 90.-DX IF (ILF.NE.0 .OR. AMOD(RLON,90.).EQ.0.0) XLAT = 90. OLAT = SIGN(AMIN1(ABS(SPLAT(I2+1)),XLAT),SPLAT(I2+1)) IGF = 0 235 IFST = 1 IGO = 0 RLAT = OLAT CALL QVEC C LATITUDE LOOP. 240 RLAT = RLAT+DLAT IGF = MAX0(IGO,IGF) CALL QVEC IF (IGO .NE. 0) GO TO 245 C THIS POINT OUTSIDE THE FRAME IF (RLAT*TSA .LE. SPLAT(I1+1)*TSA) GO TO 250 245 IF (ABS(RLAT) .LT. XLAT) GO TO 240 RLAT = SIGN(AMAX1(ABS(SPLAT(I1+1)),XLAT),SPLAT(I1+1)) C POSSIBLE NEW LATITUDE EXTREME. 250 SPLAT(I1+1) = RLAT C REVERSE LATITUDE SEARCH DIRECTION I1 = I2 I2 = MOD(I1+1,2) TSA = -PSIGN DLAT = PSIGN IF (I1 .NE. 0) GO TO 235 C LATITUDE LOOP FINISHED. IF (ABS(SPLAT(I2+1)) .LT. 90.) GO TO 255 IPF = 1 PSIGN = SIGN(1.,SPLAT(I2+1)) SPLAT(I2+1) = SPLAT(I1+1) SPLAT(I1+1) = 90.*PSIGN 255 IF (IGF .NE. 0) GO TO 230 C LONGITUDE EXTREME REACHED. IF (J .NE. 0) GO TO 260 C CHANGE LONGITUDE DIRECTION. J = 1 SPLON = RLON RLON = STLON DLON = -DLON STLON = SPLON-360. GO TO 230 C SET UP LAST LONGITUDE EXTREME 260 IF (DLON .LT. 0.0) GO TO 265 SPLON = RLON GO TO 285 265 STLON = RLON GO TO 285 C DRAW ALL MERIDIANS. 270 DLON = RGRID STLON = 0.0 SPLON = 360. RLON = 0.0 SPLAT(2) = 90. SPLAT(1) = -90. DX = AMOD(90.,RGRID) IF (DX .EQ. 0.0) DX = RGRID OLAT = 90.-DX 275 RLON = RLON+DLON IGO = 0 IFST = 1 XLAT = OLAT IF (ILF.NE.0 .OR. AMOD(RLON,90.).EQ.0.0) then If (phia.LT.0) then rlat = 0.0 xlat = 90.0 Else if (phia.GE.0) then rlat = 90.0 xlat = 0.0 End if Else rlat = xlat End if CALL QVEC 280 RLAT = RLAT-1. CALL QVEC IF (RLAT .GT. -XLAT) GO TO 280 IF (RLON .LT. SPLON) GO TO 275 C DRAW PARALLELS 285 DLAT = RGRID RLAT = AMIN1(SPLAT(2),SPLAT(1)) OLAT = AMAX1(SPLAT(2),SPLAT(1)) SPLAT(2) = FLOOR(RLAT/RGRID)*RGRID SPLAT(1) = AMIN1(CLING(OLAT/RGRID)*RGRID,90.) RLAT = AMAX1(DLAT-90.,SPLAT(2))-DLAT OLAT = AMIN1(90.-DLAT,SPLAT(1)) c if (dlon .gt. 0.0) then c stlon = stlon - dlon c else c stlon = stlon + dlon c endif DLON = 1. IF (ILF .NE. 0) IPF = 0 290 RLAT = RLAT+DLAT IF (IPF .NE. 0) DLON = 1./COS(DTR*RLAT) IGO = 0 IFST = 1 RLON = STLON CALL QVEC 295 RLON = RLON+DLON CALL QVEC IF (RLON .LE. SPLON) GO TO 295 IF (RLAT .LT. OLAT) GO TO 290 IDOT = IDOTS c CALL DASHLN (LDASH3) c CALL OPTN(2HIN,MPCOL3) C DRAW LIMB LINES 300 IDOTS = IDOT IDOT = 0 GO TO (400,330,305,335,400,340,400,400,400,345, 1 400,400,345),IPROJ C LAMBERT CONFORMAL CONIC [3] 305 DLAT = 1. RLON = PHIO+CON2 OLAT = AMAX1(-90.,SPLAT(2)-DLAT) K = CLING(SPLAT(1)-SPLAT(2)) DO 320 I=1,2 IGO = 0 IFST = 1 RLAT = OLAT CALL QVEC DO 310 J=1,K RLAT = RLAT+DLAT CALL QVEC 310 CONTINUE RLON = PHIO-CON2 320 CONTINUE GO TO 400 C ORTHOGRAPHIC [2] 330 RADIUS = 1. AXIS = 1. GO TO 350 C LAMBERT EQUAL AREA [4] 335 RADIUS = 2. AXIS = 1. GO TO 350 C AZIMUTHAL EQUDISTANT [6] 340 RADIUS = PI AXIS = 1. GO TO 350 C MOLLWEIDE [10,13] 345 RADIUS = 2. AXIS = 0.5 350 U = RADIUS V = 0.0 W = 0.0 IGO = 0 IFST = 1 DO 370 I=1,361 V = AXIS*V IF (U.LE.UMAX .AND. U.GE.UMIN .AND. V.LE.VMAX .AND. V.GE.VMIN) 1 GO TO 355 IGO = 0 GO TO 365 355 IF (IGO .NE. 0) GO TO 360 CALL FRSTPT (U,V) IGO = 1 GO TO 365 360 CALL VECTOR (U,V) 365 V = U*SINLMB+W*COSLMB U = U*COSLMB-W*SINLMB W = V 370 CONTINUE C DRAW BORDER 400 IF (JGR .GT. 0)THEN c CALL DASHLN(LDASH4) c CALL OPTN(2HIN,MPCOL4) c CALL PERIM(1,1,1,1) c call frstpt (umin,vmin) c call vector (umax,vmin) c call vector (umax,vmax) c call vector (umin,vmax) c call vector (umin,vmin) ENDIF IDOT = IDOTS RETURN 905 RETURN END C------------------------------------------------------------------------------- SUBROUTINE QCON C THIS SUBROUTINE TRANSFORMS THE POINT (RLAT,RLON), IN DEGREES, C TO (U,V) ON THE MAP PLANE DEPENDENT UPON THE PROJECTION, IPROJ. COMMON/SUPMP1/PI,TOVPI,DTR,RTD,EPS,OV90,CON1,CON2,PART COMMON/SUPMP3/POLONG,CONE,RLAT,RLON,JGR,ILF,SGN COMMON/SUPMP5/PHIOC,SINO,COSO,SINR,COSR,IPROJ COMMON/SUPMP6/UMIN,UMAX,VMIN,VMAX,UEPS,VEPS COMMON/SUPMP8/U,V,U1,V1,U2,V2 COMMON/SUPMPA/IIER c c common /supmp7/ phio,phia,igrid,idot,ilts c data resl / 2.9765625 / c c DATA OLDU,OLDV / 0.,0./ U = AMOD(RLON+PHIOC,360.)-180. GO TO (50 ,50 ,130,50 ,50 ,50 ,170,50 ,50 ,50 , 1 210,220,230),IPROJ 50 T1 = U*DTR T2 = RLAT*DTR SINPH = SIN(T1) SINLA = SIN(T2) COSPH = COS(T1) COSLA = SQRT(1.-SINLA*SINLA) TCOS = COSLA*COSPH COSA = SINLA*SINO+TCOS*COSO SINA = SQRT(CON1-COSA*COSA) C PATCH TO AVOID DIVIDE BY ZERO OVSINA = 1.E12 IF(SINA.NE.0.0)OVSINA = 1./SINA C END PATCH SINB = COSLA*SINPH*OVSINA COSB = (SINLA*COSO-TCOS*SINO)*OVSINA C PERFORM TRANSFORMATION APPROPRIATE TO THE PROJECTION GO TO (110,120,130,140,150,160,170,180,190,200),IPROJ C STEREOGRAPHIC [1] c 110 R = (1.-COSA)*OVSINA 110 ihemi = 1 if (phia .lt. 0.0) ihemi = -1 call maproj_poster (rlat, rlon, u, v, resl, polong, ihemi, 1) v = 8193 - v go to 305 C ORTHOGRAPHIC [2] 120 R = SINA IF (COSA) 320,320,300 C LAMBERT CONFORMAL CONIC [3] 130 UDIF = ABS(U-OLDU) OLDU = U CHI = 90.-SGN*RLAT IF (CHI .GE. CON2) GO TO 320 R = TAN(0.5*DTR*CHI)**CONE U = U*CONE*DTR V = -R*SGN*COS(U) U = R*SIN(U) GO TO 310 C LAMBERT EQUAL AREA [4] 140 IF (ABS(COSA+1.) .LT. 1.E-6) GO TO 320 R = (1.+COSA)*OVSINA R = 2./SQRT(1.+R*R) GO TO 300 C GNOMONIC [5] 150 IF (COSA .LE. 0.0) GO TO 320 R = SINA/COSA GO TO 300 C AZIMUTHAL EQUIDIDSANT [6] 160 IF (ABS(COSA+1.) .LT. 1.E-6) GO TO 320 R = ACOS(COSA) GO TO 300 C DUMMY -- ERROR [7] 170 IIER = 33 CALL ULIBER (IIER, 1 46H SUPMAP-ATTEMPT TO USE NON-EXISTENT PROJECTION,46) GO TO 320 C CYLINDRICAL EQUIDISTANT, ARBITRARY POLE AND ORIENTATION. 180 IF (ABS(1.-COSA*COSA) .LT. 1.E-4) GO TO 320 U = ATAN2(SINB*COSR+COSB*SINR,SINB*SINR-COSB*COSR)*RTD V = 90.-ACOS(COSA)*RTD GO TO 305 C MERCATOR, ARBITRARY POLE AND ORIENTATION. 190 IF ((1.-COSA*COSA) .LT. 2.E-6) GO TO 320 U = ATAN2(SINB*COSR+COSB*SINR,SINB*SINR-COSB*COSR) V = ALOG((1.+COSA)*OVSINA) GO TO 305 C MOLLWEIDE, ARBITRARY POLE AND ORIENTATION. 200 IF (ABS(1.-COSA*COSA) .LT. 2.E-6) GO TO 320 U = ATAN2(SINB*COSR+COSB*SINR,SINB*SINR-COSB*COSR)*TOVPI UDIF = ABS(U-OLDU) OLDU = U V = COSA C U = U*SQRT(1.-V*V) U = U*SQRT(ABS(1.-V*V)) GO TO 310 C CYLINDRICAL EQUIDISTANT FOR POLAT = ROT = 0. [11] 210 V = RLAT GO TO 305 C MERCATOR [12] 220 U = U*DTR V = ALOG(TAN(0.00872664*(RLAT+90.0001))) GO TO 305 C MOLLWEIDE [13] 230 U = U*OV90 V = SIN(RLAT*DTR) UDIF = ABS(U-OLDU) OLDU = U U = U*SQRT(1.-V*V) GO TO 310 C TERMINAL PHASE [1,2,4,5,6] 300 U = R*(SINB*COSR+COSB*SINR) V = R*(COSB*COSR-SINB*SINR) C CHECK FOR CROSSOVER 305 UDIF = ABS(U-OLDU) OLDU = U 310 VDIF = ABS(V-OLDV) OLDV = V ICROSS = 0 IF (UDIF.GT.UEPS .OR. VDIF.GT.VEPS) ICROSS = 1 RETURN C DISPENSE WITH UNDEFINED POINTS 320 U = 1.E12 ICROSS = 0 IF (ABS(U-OLDU) .GT. UEPS) ICROSS = 1 OLDU = U RETURN END C------------------------------------------------------------------------------- SUBROUTINE QVEC C THIS SUBROUTINE TRANSFORMS AND PLOTS LINE SEGMENTS FOR SUPMAP AND OTHERS C INPUTS (PASSED THROUGH COMMON.) C (RLAT,RLON) NEXT POINT TO BE PLOTTED C IFST A FLAG USED TO SIGNAL THE FIRST POINT OF A LINE SEGMENT C = 0 - START A NEW LINE C = 1 - CONTINUATION OF A LINE C OTHER VARIABLES C (U,V) NEXT POINT TRANSFORMED TO THE VIRTUAL SCREEN BY SUPCONQ C ICROSS A FLAG RETURNED BY SUPCONQ FOR CYLINDRICAL PROJECTIONS C IGO = 0 - LAST POINT NOT PLOTTED C = 1 - LAST POINT WAS PLOTTED. C (U1,V1),(U2,V2) PARAMETERS PASSED TO SUPTRP. COMMON/SUPMP4/IFST,IGO,IGOLD,ICROSS,IOUT,UOLD,VOLD COMMON/SUPMP6/UMIN,UMAX,VMIN,VMAX,UEPS,VEPS COMMON/SUPMP8/U,V,U1,V1,U2,V2 COMMON/SUPMP9/DS,DI,DSRDI SQU(X) = (X)*(X) C TRANSFORM THE POINT CALL QCON C HAVE WE FLIPPED TO OTHER SIDE OF FRAME? IF (ICROSS .NE. 0) IGO = 0 C ARE WE WITHIN THE FRAME? IF (U.GT.UMAX .OR. U.LT.UMIN .OR. V.GT.VMAX .OR. V.LT.VMIN) 1 GO TO 20 IF (IGO .EQ. 0) GO TO 30 C CONTINUE LINE C CHECK PROXIMITY TO PREVIOUS POINT. 5 IF ((SQU(U-UOLD)+SQU(V-VOLD))*DS .LE. DI) RETURN CALL VECPLT 10 UOLD = U VOLD = V IGOLD = IGO RETURN C THIS POINT LIES OUTSIDE THE FRAME 20 IGO = 0 IF (IFST .NE. 0) GO TO 65 IF (IGOLD .EQ. 0) GO TO 10 C IT WAS INSIDE - INTERPOLATE TO EDGE OF FRAME C STATUS OF LAST POINT. IF NOT INSIDE FRAME, GO ON C IF UNINTERPOLATABLE IF (ICROSS .NE. 0) GO TO 70 U1 = UOLD V1 = VOLD U2 = U V2 = V CALL SUPTRP C CHECK PROXIMITY TO PREVIOUS POINT. IF ((SQU(U-UOLD)+SQU(V-VOLD))*DS .LE. DI) GO TO 25 CALL VECPLT 25 UOLD = U2 VOLD = V2 IGOLD = 0 RETURN C THIS POINT IS WITHIN THE FRAME C IS IT THE FIRST POINT OF A LINE? 30 IF (IFST .NE. 0) GO TO 60 IF (IGOLD .EQ. 0) GO TO 50 C THE PREVIOUS POINT WAS INSIDE THE FRAME ON THE OTHER SIDE. C START A NEW LINE 40 CALL FRSTPT (U,V) IGO = 1 GO TO 10 C LAST POINT NOT IN FRAME - THIS ONE IS 50 IF (ICROSS .NE. 0) GO TO 40 C INTERPOLATE BACK TO EDGE U1 = U V1 = V U2 = UOLD V2 = VOLD CALL SUPTRP CALL FRSTPT (U,V) IGO = 1 IGOLD = 1 UOLD = U VOLD = V U = U1 V = V1 GO TO 5 C FIRST POINT ON LINE SEGMENT - CHECK FOR DUPLICATION OF END POINT 60 IF (U.NE.UOLD .OR. V.NE.VOLD) CALL FRSTPT (U,V) IGO = 1 65 IFST = 0 GO TO 10 C IGNORE UNDEFINED POINT 70 IFST = 1 GO TO 10 END C------------------------------------------------------------------------------- SUBROUTINE SUPTRP C THE INTERPOLATION ROUTINE. FINDS (U,V) ON THE EDGE OF THE FRAME NEAREST C (U1,V1). (U1,V1) MUST LIE WITHIN THE FRAME, (U2,V2) WITHOUT. COMMON/SUPMP6/UMIN,UMAX,VMIN,VMAX,UEPS,VEPS COMMON/SUPMP8/U,V,U1,V1,U2,V2 F(V) = (V-V2)*DU/DV+U2 G(U) = (U-U2)*DV/DU+V2 C FIND INDEX TO (U2,V2) C 5 | 4 | 6 C --+---+-- C 2 | 1 | 3 C --+---+-- C 8 | 7 | 9 I = 1 DU = U1-U2 DV = V1-V2 A = U2-UMIN B = U2-UMAX C = V2-VMIN D = V2-VMAX IF (A) 110,140,120 110 I = I+1 GO TO 140 120 IF (B) 140,140,130 130 I = I+2 140 IF (C) 150,200,160 150 I = I+6 GO TO 200 160 IF (D) 200,200,170 170 I = I+3 200 GO TO (900,210,220,230,240,250,260,270,280),I 210 U = UMIN GO TO 300 220 U = UMAX GO TO 300 230 V = VMAX GO TO 350 240 IF (F(VMAX)-UMIN) 210,230,230 250 IF (F(VMAX)-UMAX) 230,230,220 260 V = VMIN GO TO 350 270 IF (F(VMIN)-UMIN) 210,260,260 280 IF (F(VMIN)-UMAX) 260,260,220 C INTERPOLATE 300 V = G(U) RETURN 350 U = F(V) RETURN C ERROR EXIT 900 U = U2 V = V2 RETURN END C------------------------------------------------------------------------------- SUBROUTINE VECPLT C PLOTS THE LINE SEGMENT FROM (UOLD,VOLD) TO (U,V) C INPUTS (PASSED THROUGH COMMON) C (UOLD,VOLD) THE LAST POINT PLOTTED C (U,V) THE NEXT POINT C IDOT CONTROL FLAG [ DOT VS PLOT ] COMMON/SUPMP4/IFST,IGO,IGOLD,ICROSS,IOUT,UOLD,VOLD COMMON/SUPMP7/PHIO,PHIA,IGRID,IDOT,ILTS COMMON/SUPMP8/U,V,U1,V1,U2,V2 COMMON/SUPMP9/DS,DI,DSRDI C DO WE DOT OR PLOT? IF (IDOT .NE. 0) GO TO 10 C PLOT CALL VECTOR (U,V) RETURN C DOT 10 DU = U-UOLD DV = V-VOLD I = (ABS(DU)+ABS(DV))*DSRDI IF (I .LE. 1) GO TO 30 A = 1./FLOAT(I) I = I-1 DU = DU*A DV = DV*A UO = U VO = V U = UOLD V = VOLD DO 20 K=1,I U = U+DU V = V+DV CALL POINT (U,V) 20 CONTINUE U = UO V = VO 30 CALL POINT (U,V) RETURN END C------------------------------------------------------------------------------- SUBROUTINE SUPVEC (XLAT,XLON) C THIS SUBROUTINE ALLOWS THE USER TO DRAW LINES ON THE VIRTUAL SCREEN SET UP BY C SUPMAP, UNENCUMBERED BY DECISIONS AS TO WHETHER IT WILL BE VISIBLE THROUGH C THE WINDOW. C USE SUPFST AND SUPVEC IN EXACTLY THE SAME MANNER AS FRSTPT AND VECTOR. COMMON/SUPMP3/POLONG,CONE,RLAT,RLON,JGR,ILF,SGN RLAT = XLAT RLON = XLON CALL QVEC RETURN END C------------------------------------------------------------------------------- SUBROUTINE SUPFST (XLAT,XLON) COMMON/SUPMP3/POLONG,CONE,RLAT,RLON,JGR,ILF,SGN COMMON/SUPMP4/IFST,IGO,IGOLD,ICROSS,IOUT,UOLD,VOLD IGO = 0 IFST = 1 RLAT = XLAT RLON = XLON CALL QVEC RETURN END C------------------------------------------------------------------------------- SUBROUTINE SUPCON (XLAT,XLON,XU,XV) C THIS SUBROUTINE IS PROVIDED TO RETAIN COMPATIBILITY WITH USER'S PROGRAMS. COMMON/SUPMP3/POLONG,CONE,RLAT,RLON,JGR,ILF,SGN COMMON/SUPMP8/U,V,U1,V1,U2,V2 RLAT = XLAT RLON = XLON CALL QCON XU = U XV = V RETURN END C------------------------------------------------------------------------------- SUBROUTINE SUPMBD C FORCE-LOAD BLOCK DATA COMMON/SUPMP1/PI,TOVPI,DTR,RTD,EPS,OV90,CON1,CON2,PART COMMON/SUPMP4/IFST,IGO,IGOLD,ICROSS,IOUT,UOLD,VOLD COMMON/SUPMP9/DS,DI,DSRDI DATA CON1 / 1.00001/ DATA CON2 / 179.99999/ DATA DI / 16./ DATA DTR / 1.7453292519943E-2/ DATA EPS / 1.E-6/ DATA OV90 / 1.11111111111111E-2/ DATA PI / 3.1415926535898/ DATA RTD / 57.295779513082/ DATA TOVPI / 0.63661977236758/ DATA UOLD / 0.0 / DATA VOLD / 0.0 / DATA PART/0.9/ !SIZE OF PICTURE (90% OF SCREEN) RETURN END C------------------------------------------------------------------------------- SUBROUTINE RESUP(IER) C RECALL SUPMAP USING CONTENTS OF COMMON BLOCKS. COMMON/SUPMP1/PI,TOVPI,DTR,RTD,EPS,OV90,CON1,CON2,PART COMMON/SUPMP3/POLONG,CONE,RLAT,RLON,JGR,ILF,SGN COMMON/SUPMP6/UMIN,UMAX,VMIN,VMAX,UEPS,VEPS COMMON/SUPMP7/PHIO,PHIA,IGRID,IDOT,ILTS COMMON/SUPMPB/X1,Y1,X2,Y2 COMMON/SUPMPC/LPROJ,ROT,JLTS,LGRID,IUS LPRJ=LPROJ PHA=PHIA PLONG=POLONG RT=ROT UMN=UMIN UMX=UMAX VMN=VMIN VMX=VMAX LGRD=LGRID IU=IUS IDT=IDOT PART=AMAX1(X2-X1,Y2-Y1) CALL SUPMAP(LPRJ,PHA,PLONG,RT,UMN,UMX,VMN,VMX,-3,LGRD,IU,IDT,IER) RETURN END c c c subroutine supset (jproj, polat, polon, rot, + pl1, pl2, pl3, pl4, jlts, + jgrid, jout, jdot, jerr) c c call supmap (jproj, polat, polon, rot, + pl1, pl2, pl3, pl4, jlts, + 0, 0, 0, jerr) c c return end