subroutine getmapbounds(mapflag,maxplots,region,nplots, *olddtm,newdtm, *bw,be,bs,bn,jm,b1,b2,fn,lrv,rv0x,rv0y,rl0y) c - Cannibalized from NADCON 5 version for VERTCON 3 c - 2016 08 29: Taking in olddtm and newdtm now, and adding c - code to use that to filter out "Saint" regions in Alaska c - when plotting transformations not supported in those regions. c - 2016 08 26: Used actual mercator projection math to compute c - the exact reference vector and label locations 1/2 inch and c - 3/4 inch respectively below the S/W corner of the plot. c - Updated 2016 07 21: Two new columns added to "map.parameters", which c - have the location of the reference vector. Return a logical "lrv" c - as true if there is an optional special location for the reverence vector. c - Return as false if not. If true, return lon/lat coords of ref vector c - origin in rv0x/rv0y. If false, return zeros in those fields. c - c - Also, compute "jm" on the fly, ignoring what is in the table. All plots c - will now be forced PORTRAIT and forced no wider than 6" and no taller c - than 8", while maintaining proper X/Y ratios in a Mercator projection. c - That means, make the biggest plot possible, with the right ratio, that c - is neither wider than 6" nor taller than 8" and then, whatever the width c - of that largest plot is -- return that width in the "jm" field. c c - Subroutine to collect up the MAP boundaries c - for use in creating NADCON 5 c c - This CAN BE different than the GRID boundaries c - as such: c c - GRID boundaries will be just four values (n/s/w/e) for any region c c - MAP boundaries will allow multiple maps to be made and may or may c - not align with the GRID boundaries. Used to allow for more c - "close up" maps and such, without the need to screw up the c - MAP boundaries. c - Input: c region, maxplots, mapflag c olddtm, newdtm c - Output: c number of plots, their bounds, their GMT scales and c label deltas c - Version for NADCON 5 c - Built upon the original version used in GEOCON v2.0 c - Do not use with GEOCON v2.0 c - Broken down into sub-subroutines to make it easier c - to swap out when I make different choices. character*10 region,olddtm,newdtm real*8 bw(maxplots),be(maxplots),bs(maxplots),bn(maxplots) real*4 jm(maxplots) real*4 b1(maxplots),b2(maxplots) character*10 fn(maxplots) character*1 mapflag c - 2016 07 21: logical lrv(maxplots) real*8 rv0x(maxplots),rv0y(maxplots),rl0y(maxplots) character*3 c3x,c3y character*2 c2x,c2y real*8 pi,d2r character*200 card pi = 2.d0*dasin(1.d0) d2r = pi/180.d0 ifile = 3 open(ifile, *file='/home/dru/VERTCON/VERTCON3/Data/map.parameters', *status='old',form='formatted') c - Get numerical version of mapflag read(mapflag,*)iflag c - Get the header out of the way read(ifile,'(a)')card read(ifile,'(a)')card c - Loop over all map parameters, finding the ones relevant to c - our region and mapflag. i = 0 1 read(ifile,'(a)',end=2)card if(card( 1:10).ne.region)goto 1 read(card(14:14),*)iflag0 if(iflag0.gt.iflag)goto 1 i = i + 1 fn(i) = trim(card( 16: 25)) read(card( 27: 30),*)iwd read(card( 32: 33),*)iwm read(card( 35: 38),*)ied read(card( 40: 41),*)iem read(card( 43: 45),*)isd read(card( 47: 48),*)ism read(card( 50: 52),*)ind read(card( 54: 55),*)inm read(card( 57: 60),*)xjm read(card( 62: 64),*)ib1d read(card( 66: 67),*)ib1m read(card( 69: 71),*)ib2d read(card( 73: 74),*)ib2m c - 2016 07 21: read(card( 76: 78),'(a3)')c3x read(card( 80: 81),'(a2)')c2x read(card( 83: 85),'(a3)')c3y read(card( 87: 88),'(a2)')c2y bw(i) = dble(iwd) + dble(iwm)/60.d0 be(i) = dble(ied) + dble(iem)/60.d0 if(ind.lt.0)inm = -inm if(isd.lt.0)ism = -ism bn(i) = dble(ind) + dble(inm)/60.d0 bs(i) = dble(isd) + dble(ism)/60.d0 jm(i) = xjm b1(i) = dble(ib1d) + dble(ib1m)/60.d0 b2(i) = dble(ib2d) + dble(ib2m)/60.d0 c - 2016 07 21 (Optional reference vector locations): if(c3x.eq.'---')then lrv(i) = .false. rv0x(i) = 0 rv0y(i) = 0 else lrv(i) = .true. read(c3x,*)rv0xd read(c2x,*)rv0xm read(c3y,*)rv0yd read(c2y,*)rv0ym rv0x(i) = dble(rv0xd) + dble(rv0xm)/60.d0 rv0y(i) = dble(rv0yd) + dble(rv0ym)/60.d0 endif c - 2016 07 21 (Forcing the map to keep its proper Mercator-projected X/Y ratio, c - while filling the maximum space that does not exceed 6" wide nor 8" high) dx = (be(i) - bw(i)) * d2r q1 = dtan(bn(i)*d2r) q2 = 1.d0 / dcos(bn(i)*d2r) ! Secant yn = log(q1 + q2) q1 = dtan(bs(i)*d2r) q2 = 1.d0 / dcos(bs(i)*d2r) ! Secant ys = log(q1 + q2) dy = yn - ys ratio = dx/dy c - Max height and width c xmxht = 8 xmxht = 6 xmxwd = 6 if(xmxht*ratio .gt. xmxwd)then jm(i) = xmxwd c - And height = xmxwd / ratio... else jm(i) = xmxht*ratio c - And height = xmxht endif c - 2016 08 26 c See DRU-12, p. 56-57 c - dy/dphi is linear for latitudes south of 65. All of our plots c - have at most a 65 degree south border. So I can use that for c - a pretty good approximation to compute the latitude of the c - reference vector, in degrees, to send back as "rv0y". c First, get the ratio of dx/dy in the lower left part of the c - plot. Use 1% of the E/W span and apply that both ways. ddum = 0.01d0*(be(i) - bw(i)) ydum = bs(i) + ddum xdum = bw(i) + ddum xw = bw(i) * d2r xe = xdum * d2r dx = (xe - xw) q1 = dtan(ydum*d2r) q2 = 1.d0 / dcos(ydum*d2r) ! Secant yn = log(q1 + q2) q1 = dtan(bs(i)*d2r) q2 = 1.d0 / dcos(bs(i)*d2r) ! Secant ys = log(q1 + q2) c - dy is in radians dy = yn - ys c - Now I have a new "dx" and "dy" representing the lower left 1% of the plot, c - in radians/radians. Compute the new ratio for this part. ratio = dx/dy c - Since radians to inches does NOT change in the E/W direction, and since c - I know exactly how wide the plot is, I can compute a radians/inch conversion c - in E/W, for 1% of the width of the plot. r2iew = (0.01d0 * jm(i)) / dx c - The ratio of dx/dy, is going to be identical to the ratio of r2iew/r2ins: r2ins = r2iew / ratio c - Convert r2ins into degrees-to-inches d2ins = r2ins * d2r c - Now, the south edge is known in degrees. I want the reference vector c - to be 1/2 inch below that. I want the reference vector label to c - be an additional 1/4 inches down. rv0y(i) = bs(i) - (0.50d0)/d2ins rl0y(i) = bs(i) - (0.65d0)/d2ins c - And set the reference vector's left edge to equal the left edge of the plot rv0x(i) = bw(i) c - Flags for debugging write(6,8) trim(fn(i)) 8 format('***************************',/, * 'Inside gridstats.f....',/, * ' Reference Vector Computations for ', * 'Sub-Region: ',a) write(6,*) ' yn,ys,dy = ',yn,ys,dy write(6,*) ' xe,xw,dx = ',xe,xw,dx write(6,*) ' E/W width (jm) = ',jm(i) write(6,*) ' SW ratio = ',ratio write(6,*) ' SW bounds,Lat = ',ydum,bs(i) write(6,*) ' SW bounds,Lon = ',xdum,bw(i) write(6,*) ' r2i E/W = ',r2iew write(6,*) ' r2i N/S = ',r2ins write(6,*) ' d2i N/S = ',d2ins write(6,*) ' Southwest corner = ',bs(i),bw(i) write(6,*) ' Ref Vector = ',rv0y(i),rv0x(i) write(6,*) ' Ref Label = ',rl0y(i),rv0x(i) write(6,9) 9 format('***************************') goto 1 2 nplots = i close(ifile) if(nplots.eq.0)then write(6,100)trim(region) stop 10001 endif 100 format(6x,'FATAL Subroutine getmapbounds: Unknown region: ',a) return end