c - Program "makeplotfiles03" c - Cannibalized version from NADCON 5 for VERTCON 3 c - NADCON 5 comments below c------------------------------------------------------------- c - 2016 10 19 c - FIX for the HARN/FBN transformation. As it stood, the "gridstats" was c - returning "0.0" as the median of the post-masked "ete" grid (which c - is true, but unfortunate.) I've changed the call to "gridstats" c - to send it the "PREMASKED" version of the "ete" grid, so the median c - won't be zero. c - 2016 08 26 c - Added new code to do reference vectors consistently c - See DRU-12, p. 56-57 c - Fixed an error where "lorvogehtm" was declared real*8 but past 72 column, c - so defaulting to integer*4 and coming out as "0.000" on plots c - Also changing the call to "getmapbounds" to give it "olddatum" and "newdatum" c - to aide in filtering out things like the Saint regions in Alaska c - for unsupported transformations. c - 2016 08 02: Changed the color palette for "09" grids from 2xMedian to 3xMedian c - 2016 07 29: Dropped code about personalized reference vectors and just c - let them be below/outside map c - 2016 08 01: Moved "gridstats" to subroutines c - Also, completely removed the in-program computations of c the color palette, and instead relied on "cpt" and "cpt2" c as per "makeplotfiles02" c - 2016 07 21: c - Added code to allow for optional placement of reference vectors, coming from c "map.parameters" as read in subroutine "getmapbounds" c c - Part of the NADCON5 process c c - Program presumes, as input c - the "old datum", "new datum" , "region", "gridsec" and "mapflag" c - arguments fed into "doit4.bat" c - Examples: c olddatum = 'ussd' c newdatum = 'nad27' c region = 'conus' c gridsec= '900' c mapflag = '0' c c - Updated 2016 01 21 to fix the CPT for the "09" (data noise) c - grids so that (cpthi - cptin) is exactly divisible by cptin c c - Updated 2015 11 09 to adopt new naming scheme (yes, again) c - (See DRU-11, p. 150), as well as creating plots of the c - total error grid. c - Latest version: 2015 10 07 which had new naming scheme and c - dual-computations of arcseconds and meters for lat/lon c - See DRU-11, pl. 139. c - Creates a batch file called "gmtbat06.(olddtm).(newdtm).(region).(igridsec).(mapflag)" c - That batch file will create JPGs of: c - 1)Color Plots of the rddlat/rddlon/rddeht grids c - 2)B/W plots of coverage of RMS'd differential vectors that went into the grid c - 3)B/W plots of RMS'd differential vectors that went into the grid implicit real*8(a-h,o-z) parameter(maxplots=60) character*35 fname0 character*34 dirnam character*10 od,nd character*2 state,stdum character*3 ele,elelat,elelon,eleoht,ele0 character*5 agridsec character*5 agridsec2 integer*4 igridsec integer*4 igridsec2 character*200 suffix1,suffix2,suffix3 character*200 suffix2t09,suffix2d3 character*200 wfname,gmtfile c - ".b" of RMS'd double difference data character*200 bfnvsrddlat,bfnvsrddlon,bfnvmrddoht character*200 bfnvmrddlat,bfnvmrddlon c - "grd" of RMS'd double difference data character*200 rfnvsrddlat,rfnvsrddlon,rfnvmrddoht character*200 rfnvmrddlat,rfnvmrddlon c - "coverage" of RMS'd data (needs to be created in "myrms.f") character*200 gfncvrddlat,gfncvrddlon,gfncvrddoht c - RMS'd double difference vectors character*200 gfnvsrddlat,gfnvsrddlon,gfnvmrddoht character*200 gfnvmrddlat,gfnvmrddlon c - "Method Noise" grids: character*200 sadbfnvmtcdlat1000,sadbfnvmtcdlon1000 character*200 sadbfnvmtcdoht1000 character*200 sadbfnvstcdlat1000,sadbfnvstcdlon1000 c -" Total Error" grids: character*200 bfnvsetelat,bfnvsetelon character*200 bfnvmetelat,bfnvmetelon character*200 bfnvmeteoht c - 2016 10 19: c - For the HARN/FBN situation *only* c - "Total Error" grids, pre-mask: character*200 bfnvsetelat0,bfnvsetelon0 character*200 bfnvmetelat0,bfnvmetelon0 character*200 bfnvmeteeht0 c - 2015 10 09: character*200 gfnvsrddhor,gfnvmrddhor c - Stats file character*200 sfn c - GMT stuff real*8 bw(maxplots),be(maxplots),bs(maxplots),bn(maxplots) real*4 jm(maxplots) real*4 b1(maxplots),b2(maxplots) character*10 fn(maxplots) c - 2016 07 21: logical lrv(maxplots) real*8 rv0x(maxplots),rv0y(maxplots),rl0y(maxplots) c - Stuff that is in the work file records: character*6 pid character*2 state character*1 rejlat,rejlon,rejeht real*8 xlath,xlonh,xohth real*8 dlatsec,dlonsec,dohtm real*8 dlatm,dlonm character*10 olddtm,newdtm,region character*3 filter c - Units to put on scale on GMT color plots character*17 scunlat,scunlon,scunoht character*1 mapflag real*8 lorvopc,lorvoglats,lorvoglatm,lorvoglons,lorvoglonm real*8 lorvogohtm c - 2015 10 09 real*8 lorvoghors,lorvoghorm logical v2special c ------------------------------------------------------------------ c - BEGIN PROGRAM c ------------------------------------------------------------------ write(6,1001) 1001 format('BEGIN program makeplotfiles03.f') c ------------------------------------------------------------------ c - User-supplied input c ------------------------------------------------------------------ read(5,'(a)')olddtm read(5,'(a)')newdtm read(5,'(a)')region read(5,'(a)')filter read(5,'(a)')agridsec read(5,'(a)')agridsec2 read(5,'(a)')mapflag v2special = .false. if(trim(region).eq.'conus' .and. * trim(olddtm).eq.'ngvd29' .and. * trim(newdtm).eq.'navd88') v2special = .true. c ------------------------------------------------------------------ c - Some necessary constants c ------------------------------------------------------------------ lorvopc = 1.d0 csm = 3.d0 pi = 2.d0*dasin(1.d0) d2r = pi/180.d0 MultiplierLorvog = 2 c ------------------------------------------------------------------ c - Generate the suffixes used in all our files c ------------------------------------------------------------------ read(agridsec,*)igridsec read(agridsec2,*)igridsec2 suffix1=trim(olddtm)//'.'//trim(newdtm)//'.'//trim(region) *//'.'//trim(filter) c suffix2=trim(suffix1)//'.'//trim(agridsec) suffix2=trim(suffix1)//'.'//trim(agridsec)//'.'//trim(agridsec2) suffix3=trim(suffix2)//'.'//trim(mapflag) suffix2d3=trim(suffix2)//'.d3' suffix2t09=trim(suffix2)//'.09' c ------------------------------------------------------------------ c - Got clever here...I stored # of double differences and the RMS of those c - vectors into "dvstats...." back in the running of "checkgrid.f". c - Now let's re-acquire this data and use it to scale double differenced vectors c - in this program. Granted, this won't be EXACTLY the same c - as what is in the data (the RMS'd DD vectors are RMS'd in cells/bins, c - but what is in "dvstats" is just an area-wide RMS value) but it's c - good enough to scale things in the ballpark. sfn = 'dvstats.'//trim(suffix2) open(1,file=sfn,status='old',form='formatted') write(6,1004) 1004 format(6x,'makeplotfiles03.f: Acquiring vectors ', *'stats from the dvstats file') c read(1,*)nlat,rlats c read(1,*)nlon,rlons read(1,*)noht,rohtm c read(1,*)nlat0,rlatm c read(1,*)nlon0,rlonm c - 2015 10 09 c read(1,*)nhor ,rhors c read(1,*)nhor0,rhorm c if(nlat0.ne.nlat)stop 10001 c if(nlon0.ne.nlon)stop 10002 c if(nhor0.ne.nhor)stop 10003 c if(nhor.ne.nlat)stop 10004 c if(nhor.ne.nlon)stop 10005 write(6,893)noht,rohtm 893 format(6x,'makeplotfiles03.f: Vector Stats: ',/, *10x,'Number of Double Differenced Vectors in OHT: ',i10,/, *10x,'RMS length (meters) : ',f10.3) c write(6,893)nlat,rlats,rlatm, c *nlon,rlons,rlonm, c *nhor,rhors,rhorm, c *neht,rehtm c 893 format(6x,'makeplotfiles03.f: Vector Stats: ',/, c *10x,'Number of Double Differenced Vectors in LAT: ',i10,/, c *10x,'RMS length (arcseconds) : ',f10.3,/, c *10x,'RMS length (meters) : ',f10.3,/, c *10x,'Number of Double Differenced Vectors in LON: ',i10,/, c *10x,'RMS length (arcseconds) : ',f10.3,/, c *10x,'RMS length (meters) : ',f10.3,/, c *10x,'Number of Double Differenced Vectors in HOR: ',i10,/, c *10x,'RMS length (arcseconds) : ',f10.3,/, c *10x,'RMS length (meters) : ',f10.3,/, c *10x,'Number of Double Differenced Vectors in EHT: ',i10,/, c *10x,'RMS length (meters) : ',f10.3) c ------------------------------------------------------------------ c -------------------------------------------------------------- c - GMT Batch file: Open the file c -------------------------------------------------------------- gmtfile = 'gmtbat06.'//trim(suffix3) open(99,file=gmtfile,status='new',form='formatted') write(6,1011)trim(gmtfile) 1011 format(6x,'makeplotfiles03.f: Creating GMT batch file ',a) write(99,1030)trim(gmtfile) 1030 format('echo BEGIN batch file ',a) c - The "method noise" grids c - UNNEEDED! c sadbfnvmtcdlat1000 = 'vmtcdlat.'//trim(suffix2d3)//'.b' c sadbfnvstcdlat1000 = 'vstcdlat.'//trim(suffix2d3)//'.b' c sadbfnvmtcdlon1000 = 'vmtcdlon.'//trim(suffix2d3)//'.b' c sadbfnvstcdlon1000 = 'vstcdlon.'//trim(suffix2d3)//'.b' c sadbfnvmtcdoht1000 = 'vmtcdoht.'//trim(suffix2d3)//'.b' c - The "data noise" grids bfnvsrddlat = 'vsrddlat.'//trim(suffix2t09)//'.b' bfnvsrddlon = 'vsrddlon.'//trim(suffix2t09)//'.b' bfnvmrddoht = 'vmrddoht.'//trim(suffix2t09)//'.b' bfnvmrddlat = 'vmrddlat.'//trim(suffix2t09)//'.b' bfnvmrddlon = 'vmrddlon.'//trim(suffix2t09)//'.b' c - The "total noise" grids bfnvsetelat ='vsetelat.'//trim(suffix2)//'.b' bfnvsetelon ='vsetelon.'//trim(suffix2)//'.b' bfnvmetelat ='vmetelat.'//trim(suffix2)//'.b' bfnvmetelon ='vmetelon.'//trim(suffix2)//'.b' bfnvmeteoht ='vmeteoht.'//trim(suffix2)//'.b' c -------------------------------------------------------------- c - GMT Batch file: Get color palatte for the "total error" c - grids. Keep this separate from the color palette for c - the "data noise" grids. c -------------------------------------------------------------- c - VERTCON3: Skip to the OHT part... goto 880 c - 2016 08 01: All below changed to rely on "cpt2" for "ete" grids if(nlat.ne.0)then c - 2016 10 19 c - Lazy coding here: You can ONLY have "olddtm" = "nad83_harn" if c - you also have "nad83_fbn" as "newdtm" and "conus" c - as "region", so I won't bother checking all 3: if(trim(olddtm).ne.'nad83_harn')then write(6,1012)trim(bfnvmetelat) call gridstats(bfnvmetelat,ave,std,xmd) else bfnvmetelat0 = trim(bfnvmetelat)//'.premask' write(6,1012)trim(bfnvmetelat0) call gridstats(bfnvmetelat0,ave,std,xmd) endif avelatmE = ave stdlatmE = std xmdlatmE = xmd call cpt2(xmdlatmE,2.d0, * cptlolatmE,cpthilatmE,cptinlatmE) c cptinlatmE=onzd2(std) c scaledave =onzd2(ave) c cptlolatmE = 0.d0 c cpthilatmE = scaledave + csm*cptinlatmE c - 2016 10 19 c - Lazy coding here: You can ONLY have "olddtm" = "nad83_harn" if c - you also have "nad83_fbn" as "newdtm" and "conus" c - as "region", so I won't bother checking all 3: if(trim(olddtm).ne.'nad83_harn')then write(6,1012)trim(bfnvsetelat) call gridstats(bfnvsetelat,ave,std,xmd) else bfnvsetelat0 = trim(bfnvsetelat)//'.premask' write(6,1012)trim(bfnvsetelat0) call gridstats(bfnvsetelat0,ave,std,xmd) endif avelatsE = ave stdlatsE = std xmdlatsE = xmd call cpt2(xmdlatsE,2.d0, * cptlolatsE,cpthilatsE,cptinlatsE) c cptinlatsE=onzd2(std) c scaledave =onzd2(ave) c cptlolatsE = 0.d0 c cpthilatsE = scaledave + csm*cptinlatsE endif if(nlon.ne.0)then c - 2016 10 19 c - Lazy coding here: You can ONLY have "olddtm" = "nad83_harn" if c - you also have "nad83_fbn" as "newdtm" and "conus" c - as "region", so I won't bother checking all 3: if(trim(olddtm).ne.'nad83_harn')then write(6,1012)trim(bfnvmetelon) call gridstats(bfnvmetelon,ave,std,xmd) else bfnvmetelon0 = trim(bfnvmetelon)//'.premask' write(6,1012)trim(bfnvmetelon0) call gridstats(bfnvmetelon0,ave,std,xmd) endif avelonmE = ave stdlonmE = std xmdlonmE = xmd call cpt2(xmdlonmE,2.d0, * cptlolonmE,cpthilonmE,cptinlonmE) c cptinlonmE=onzd2(std) c scaledave =onzd2(ave) c cptlolonmE = 0.d0 c cpthilonmE = scaledave + csm*cptinlonmE c - 2016 10 19 c - Lazy coding here: You can ONLY have "olddtm" = "nad83_harn" if c - you also have "nad83_fbn" as "newdtm" and "conus" c - as "region", so I won't bother checking all 3: if(trim(olddtm).ne.'nad83_harn')then write(6,1012)trim(bfnvsetelon) call gridstats(bfnvsetelon,ave,std,xmd) else bfnvsetelon0 = trim(bfnvsetelon)//'.premask' write(6,1012)trim(bfnvsetelon0) call gridstats(bfnvsetelon0,ave,std,xmd) endif avelonsE = ave stdlonsE = std xmdlonsE = xmd call cpt2(xmdlonsE,2.d0, * cptlolonsE,cpthilonsE,cptinlonsE) c cptinlonsE=onzd2(std) c scaledave =onzd2(ave) c cptlolonsE = 0.d0 c cpthilonsE = scaledave + csm*cptinlonsE endif c - VERTCON3: Til here 880 continue if(noht.ne.0)then c - 2016 10 19 c - Lazy coding here: You can ONLY have "olddtm" = "nad83_harn" if c - you also have "nad83_fbn" as "newdtm" and "conus" c - as "region", so I won't bother checking all 3: c - VERTCON3: No need for the HARN/FBN code c if(trim(olddtm).ne.'nad83_harn')then write(6,1012)trim(bfnvmeteoht) call gridstats(bfnvmeteoht,ave,std,xmd) c else c bfnvmeteeht0 = trim(bfnvmeteeht)//'.premask' c write(6,1012)trim(bfnvmeteeht0) c call gridstats(bfnvmeteeht0,ave,std,xmd) c endif aveohtmE = ave stdohtmE = std xmdohtmE = xmd call cpt2(xmdohtmE,2.d0, * cptloohtmE,cpthiohtmE,cptinohtmE) c cptinehtmE=onzd2(std) c scaledave =onzd2(ave) c cptloehtmE = 0.d0 c cpthiehtmE = scaledave + csm*cptinehtmE endif c -------------------------------------------------------------- c - GMT Batch file: Determine Number of areas to plot (=nplots) c - and Boundaries of plots and other stuff. c - Report the information out. c -------------------------------------------------------------- c - 2016 08 29: call getmapbounds(mapflag,maxplots,region,nplots, *olddtm,newdtm, *bw,be,bs,bn,jm,b1,b2,fn,lrv,rv0x,rv0y,rl0y) c - 2016 08 26, DRU-12, p.56-57 c call getmapbounds(mapflag,maxplots,region,nplots, c *bw,be,bs,bn,jm,b1,b2,fn,lrv,rv0x,rv0y,rl0y) c - See DRU-11, p. 126 for choices on grid boundaries for NADCON v5.0 c - 2016 07 21 c call getmapbounds(mapflag,maxplots,region,nplots, c *bw,be,bs,bn,jm,b1,b2,fn,lrv,rv0x,rv0y) c call getmapbounds(mapflag,maxplots,region,nplots, c *bw,be,bs,bn,jm,b1,b2,fn) write(6,1006)trim(region) 1006 format(6x,'makeplotfiles03.f: Calling getmapbounds for region ',a) c --------------------------------------------------------------- c - Compute and report various things about our coverage c - and vector plots. c --------------------------------------------------------------- c - How many centimeters long will our reference arrows be? c - Base this on the actual data if(nlat.ne.0)then c lorvoglats = onzd2(MultiplierLorvog*rlats) c gs2pclat = lorvopc / lorvoglats c lorvoglatm = onzd2(MultiplierLorvog*rlatm) c gm2pclat = lorvopc / lorvoglatm endif if(nlon.ne.0)then c lorvoglons = onzd2(MultiplierLorvog*rlons) c gs2pclon = lorvopc / lorvoglons c lorvoglonm = onzd2(MultiplierLorvog*rlonm) c gm2pclon = lorvopc / lorvoglonm endif if(noht.ne.0)then lorvogohtm = onzd2(MultiplierLorvog*rohtm) gm2pcoht = lorvopc / lorvogohtm endif if(nlon.ne.0 .and. nlat.ne.0)then c lorvoghors = onzd2(MultiplierLorvog*rhors) c gs2pchor = lorvopc / lorvoghors c lorvoghorm = onzd2(MultiplierLorvog*rhorm) c gm2pchor = lorvopc / lorvoghorm endif c - Report: write(6,894)nplots 894 format(6x,'makeplotfiles03.f: Info about plots:',/, *8x,'Number of sub-area plot sets to cover this region: ',i2) do 895 i=1,nplots dns = bn(i) - bs(i) dew = be(i) - bw(i) write(6,896)i,bs(i),bn(i),bw(i),be(i),dns,dew if(nlat.ne.0)then c write(6,897)'lat',lorvoglatm,lorvoglats,gm2pclat,gs2pclat else c write(6,898)'lat' endif if(nlon.ne.0)then c write(6,897)'lon',lorvoglonm,lorvoglons,gm2pclon,gs2pclon else c write(6,898)'lon' endif if(noht.ne.0)then write(6,899)'oht',lorvogohtm,gm2pcoht else write(6,898)'oht' endif if(nlat.ne.0 .and. nlon.ne.0)then c write(6,897)'hor',lorvoghorm,lorvoghors,gm2pchor,gs2pchor else c write(6,898)'hor' endif 895 continue 896 format( *8x,'Plot # ',i2,/, *10x,'S/N/W/E/N-S/E-W = ',6f7.1) 897 format( *10x,a3,' plots',13x,':',/, *12x,'Reference Vector ( m) = ',f10.2,/, *12x,'Reference Vector ( s) = ',f10.6,/, *12x,'Ground M to Paper cm = ',f20.10,/, *12x,'Ground S to Paper cm = ',f20.10) 899 format( *10x,a3,' plots',13x,':',/, *12x,'Reference Vector ( m) = ',f10.2,/, *12x,'Ground M to Paper cm = ',f20.10) 898 format( *10x,a3,' plots: No data available for plotting') c ------------------------------------------------------------------ c - pvlon and pvlat are the percentage of the total lon/lat c - span of the plot, from which the reference vector c - begins, starting with the Lower Left corner c ------------------------------------------------------------------ pvlon = 10.d0 pvlat = 10.d0 pvlat = (pvlat / 100.d0) pvlon = (pvlon / 100.d0) c -------------------------------------------------------------- c - GMT Batch file: Get color palatte for the "data noise", aka "09" c - grids. c -------------------------------------------------------------- c ------------------------------------------------------------------ c - Because our color palette will depend upon the statistics of c - our grids, we need to access those grids, and store a few c - statistics. c ------------------------------------------------------------------ rfnvsrddlat = 'vsrddlat.'//trim(suffix2t09)//'.grd' rfnvsrddlon = 'vsrddlon.'//trim(suffix2t09)//'.grd' rfnvmrddoht = 'vmrddoht.'//trim(suffix2t09)//'.grd' rfnvmrddlat = 'vmrddlat.'//trim(suffix2t09)//'.grd' rfnvmrddlon = 'vmrddlon.'//trim(suffix2t09)//'.grd' gfnvsrddlat = 'vsrddlat.'//trim(suffix2) gfnvsrddlon = 'vsrddlon.'//trim(suffix2) gfnvmrddoht = 'vmrddoht.'//trim(suffix2) gfnvmrddlat = 'vmrddlat.'//trim(suffix2) gfnvmrddlon = 'vmrddlon.'//trim(suffix2) gfnvsrddhor = 'vsrddhor.'//trim(suffix2) gfnvmrddhor = 'vmrddhor.'//trim(suffix2) gfncvrddlat = 'cvrddlat.'//trim(suffix2) gfncvrddlon = 'cvrddlon.'//trim(suffix2) gfncvrddoht = 'cvrddoht.'//trim(suffix2) c - NOTE ON COLOR PALLETTES, BELOW: c - Because all of the differential vectors are RMS'd values at c - this point, they will have a minimum value of zero. As such, c - set the "cptlo" values to zero. c - 2016 08 01: Changed to use "cpt2" for "09" (RMS) grids c - 2016 08 02: Changed the multiplier in the 09 grids from 2 to 3 c - VERTCON3: Skip lat/lon goto 881 if(nlat.ne.0)then write(6,1012)trim(bfnvsrddlat) call gridstats(bfnvsrddlat,ave,std,xmd) avelats = ave stdlats = std xmdlats = xmd call cpt2(xmdlats,3.d0, * cptlolats,cpthilats,cptinlats) write(6,1012)trim(bfnvmrddlat) call gridstats(bfnvmrddlat,ave,std,xmd) avelatm = ave stdlatm = std xmdlatm = xmd call cpt2(xmdlatm,3.d0, * cptlolatm,cpthilatm,cptinlatm) endif if(nlon.ne.0)then write(6,1012)trim(bfnvsrddlon) call gridstats(bfnvsrddlon,ave,std,xmd) avelons = ave stdlons = std xmdlons = xmd call cpt2(xmdlons,3.d0, * cptlolons,cpthilons,cptinlons) write(6,1012)trim(bfnvmrddlon) call gridstats(bfnvmrddlon,ave,std,xmd) avelonm = ave stdlonm = std xmdlonm = xmd call cpt2(xmdlonm,3.d0, * cptlolonm,cpthilonm,cptinlonm) endif c - VERTCON 3: Til here 881 continue if(noht.ne.0)then write(6,1012)trim(bfnvmrddoht) call gridstats(bfnvmrddoht,ave,std,xmd) aveohtm = ave stdohtm = std xmdohtm = xmd call cpt2(xmdohtm,3.d0, * cptloohtm,cpthiohtm,cptinohtm) endif 1012 format(6x,'makeplotfiles03.f: Grabbing stats of grid: ',a) c ------------------------------------------------------------------ c ------------------------------------------------------------------ c ------------------------------------------------------------------ c - MAIN LOOP: Loop over number of sub-areas (nplots) I've chosen c - for this region (see "getmapbounds.f" subroutine). c ------------------------------------------------------------------ c ------------------------------------------------------------------ c ------------------------------------------------------------------ do 2001 ij=1,nplots write(99,990)trim(region),trim(fn(ij)), * trim(region),trim(fn(ij)) c -------------------- c - B/W COVERAGE PLOTS c -------------------- c - Make B/W Coverage Plots of RMS'd data c - Latitude c if(nlat.ne.0) c * call bwplotcv('lat',gfncvrddlat,bw,be,bs,bn,jm, c * b1,b2,maxplots,olddtm,newdtm,region,'LAT',ij, c * igridsec,fn) c - Longitude c if(nlon.ne.0) c * call bwplotcv('lon',gfncvrddlon,bw,be,bs,bn,jm, c * b1,b2,maxplots,olddtm,newdtm,region,'LON',ij, c * igridsec,fn) c - VERTCON3: c - Orthometric Height c - Note everything is on "igridsec2" now if(noht.ne.0) * call bwplotcv('oht',gfncvrddoht,bw,be,bs,bn,jm, * b1,b2,maxplots,olddtm,newdtm,region,filter, * 'OHT',ij, * igridsec2,fn) c ------------------ c - B/W VECTOR PLOTS c ------------------ c - VERTCON3: B/W Vector plots were replaced with color dot plots for easier reading. c - VERTCON3: Added "filter" to call, and have new subroutine "coplotvc" now... c (even though "vc" means "vector"...just hand-wave that colored dots replace "vectors") c - Where is the start of the reference vector: c - 2016 07 29: Decided to scrap personalized locations c - and just put all reference vectors outside/below c - the plots c - 2016 08 26 xvlon = rv0x(ij) xvlat = rv0y(ij) xllon = xvlon xllat = rl0y(ij) c - 2016 07 21: c if(lrv(ij))then c xvlon = rv0x(ij) c xvlat = rv0y(ij) c else cccc xvlon = bw(ij) + (pvlon*(be(ij)-bw(ij))) c xvlat = bs(ij) + (pvlat*(bn(ij)-bs(ij))) cccc xvlat = bs(ij) - (pvlat*(bn(ij)-bs(ij))) c endif c - Where is the start of the label for the ref vector: cccc xllon = xvlon c - To prevent the label from going off the page on c - plots that have very little N/S span, put in a c - failsafe c if(lrv(ij))then c xllat = xvlat - 0.1d0 c else cccc if(bn(ij)-bs(ij).lt.2.0)then c xllat = bs(ij) + ((pvlat*0.75d0)*(bn(ij)-bs(ij))) cccc xllat = bs(ij) - ((pvlat*1.25d0)*(bn(ij)-bs(ij))) cccc else c xllat = bs(ij) + (pvlat*(bn(ij)-bs(ij))) - 0.1d0 cccc xllat = bs(ij) - (pvlat*(bn(ij)-bs(ij))) - 0.1d0 cccc endif c endif c - Make B/W Vector Plots of RMS'd data c - Latitude if(nlat.ne.0)then c call bwplotvc('lat',gfnvsrddlat,bw,be,bs,bn,jm,b1,b2,maxplots, c * olddtm,newdtm,region,'LAT',ij,xvlon,xvlat,xllon,xllat, c * lorvoglats,lorvopc,igridsec,fn) c call bwplotvc('lat',gfnvmrddlat,bw,be,bs,bn,jm,b1,b2,maxplots, c * olddtm,newdtm,region,'LAT',ij,xvlon,xvlat,xllon,xllat, c * lorvoglatm,lorvopc,igridsec,fn) endif c - Longitude if(nlon.ne.0)then c call bwplotvc('lon',gfnvsrddlon,bw,be,bs,bn,jm,b1,b2,maxplots, c * olddtm,newdtm,region,'LON',ij,xvlon,xvlat,xllon,xllat, c * lorvoglons,lorvopc,igridsec,fn) c call bwplotvc('lon',gfnvmrddlon,bw,be,bs,bn,jm,b1,b2,maxplots, c * olddtm,newdtm,region,'LON',ij,xvlon,xvlat,xllon,xllat, c * lorvoglonm,lorvopc,igridsec,fn) endif c - Ellipsoid height if(noht.ne.0)then c call bwplotvc('eht',gfnvmrddeht,bw,be,bs,bn,jm,b1,b2,maxplots, c * olddtm,newdtm,region,'EHT',ij,xvlon,xvlat,xllon,xllat, c * lorvogehtm,lorvopc,igridsec,fn) c - VERTCON3: call coplotvc('oht',gfnvmrddoht,bw,be,bs,bn, * jm,b1,b2,maxplots, * olddtm,newdtm,region,filter,'OHT',ij,xvlon, * xvlat,xllon,xllat, * lorvogohtm,lorvopc,igridsec,igridsec2,fn, * cptloohtm,cpthiohtm,cptinohtm) endif c - Horizontal if(nlon.ne.0 .and. nlat.ne.0)then c call bwplotvc('hor',gfnvsrddhor,bw,be,bs,bn,jm,b1,b2,maxplots, c * olddtm,newdtm,region,'HOR',ij,xvlon,xvlat,xllon,xllat, c * lorvoghors,lorvopc,igridsec,fn) c call bwplotvc('hor',gfnvmrddhor,bw,be,bs,bn,jm,b1,b2,maxplots, c * olddtm,newdtm,region,'HOR',ij,xvlon,xvlat,xllon,xllat, c * lorvoghorm,lorvopc,igridsec,fn) endif c ------------------------------------------------ c - Color gridded data of RMS grids ("data noise") c ------------------------------------------------ c - Make color plots of gridded, RMS'd data, with no c - points. write(99,3101) 3101 format('echo Color Plots of RMS data...') c - Latitude if(nlat.ne.0)then c call coplot('lat',bfnvsrddlat,bw,be,bs,bn,jm,b1,b2,maxplots, c * olddtm,newdtm,region,'LAT',ij,cptlolats,cpthilats,cptinlats, c * suffix2t09,igridsec,fn) c call coplot('lat',bfnvmrddlat,bw,be,bs,bn,jm,b1,b2,maxplots, c * olddtm,newdtm,region,'LAT',ij,cptlolatm,cpthilatm,cptinlatm, c * suffix2t09,igridsec,fn) endif c - Longitude if(nlon.ne.0)then c call coplot('lon',bfnvsrddlon,bw,be,bs,bn,jm,b1,b2,maxplots, c * olddtm,newdtm,region,'LON',ij,cptlolons,cpthilons,cptinlons, c * suffix2t09,igridsec,fn) c call coplot('lon',bfnvmrddlon,bw,be,bs,bn,jm,b1,b2,maxplots, c * olddtm,newdtm,region,'LON',ij,cptlolonm,cpthilonm,cptinlonm, c * suffix2t09,igridsec,fn) endif c - Ellipsoid Height c - VERTCON3: updated call to coplot if(noht.ne.0)then call coplot('oht',bfnvmrddoht,bw,be,bs,bn,jm,b1,b2,maxplots, * olddtm,newdtm,region,filter,'OHT', * ij,cptloohtm,cpthiohtm,cptinohtm, * suffix2t09,igridsec,igridsec2,fn) endif c ------------------------------------------------ c - Color gridded data of "Total Error" grids c ------------------------------------------------ write(99,3102) 3102 format('echo Color Plots of Total Error grids...') c - Latitude if(nlat.ne.0)then c call coplot('lat',bfnvsetelat,bw,be,bs,bn,jm,b1,b2,maxplots, c * olddtm,newdtm,region,'LAT',ij,cptlolatsE,cpthilatsE, c * cptinlatsE, c * suffix2,igridsec,fn) c call coplot('lat',bfnvmetelat,bw,be,bs,bn,jm,b1,b2,maxplots, c * olddtm,newdtm,region,'LAT',ij,cptlolatmE,cpthilatmE, c * cptinlatmE, c * suffix2,igridsec,fn) endif c - Longitude if(nlon.ne.0)then c call coplot('lon',bfnvsetelon,bw,be,bs,bn,jm,b1,b2,maxplots, c * olddtm,newdtm,region,'LON',ij,cptlolonsE,cpthilonsE, c * cptinlonsE, c * suffix2,igridsec,fn) c call coplot('lon',bfnvmetelon,bw,be,bs,bn,jm,b1,b2,maxplots, c * olddtm,newdtm,region,'LON',ij,cptlolonmE,cpthilonmE, c * cptinlonmE, c * suffix2,igridsec,fn) endif c - Orthometric Height c - VERTCON3: Updated call to coplot if(noht.ne.0)then call coplot('oht',bfnvmeteoht,bw,be,bs,bn,jm,b1,b2,maxplots, * olddtm,newdtm,region,filter,'OHT', * ij,cptloohtmE,cpthiohtmE, * cptinohtmE, * suffix2,igridsec,igridsec2,fn) endif 2001 continue 990 format( *'# ------------------------------',/, *'# Plots for region: ',a,', sub-region: ',a,/, *'# ------------------------------',/, *'echo Creating plots for region: ',a,', sub-region: ',a) write(99,1031)trim(gmtfile) 1031 format('echo END batch file ',a) close(99) write(6,9999) 9999 format('END program makeplotfiles03.f') end c c ---------------------------------------------------- c include 'Subs/getmapbounds.f' include 'Subs/getmag.f' include 'Subs/coplot.f' c include 'Subs/bwplotvc.f' include 'Subs/coplotvc.f' include 'Subs/bwplotcv.f' include 'Subs/plotcoast.f' include '/home/dru/Subs/onzd2.f' include 'Subs/gridstats.f' include 'Subs/cpt2.f' include '/home/dru/NumRec/Modified/select2.for'