c - Program "makeplotfiles02" c - 2018 05 16: Updated to allow for TWO grid spacings: c Under "v2special" they will be the spacings of the RZ and CD grids, independent of each other c Otherwise, they will be IDENTICAL c c c - 2018 05 04 - Undid any previous decisions to falsely store "rz" data under variable names with "cd" in them... c - See DRU-12, p. 133 c c - Note, there's a bit of schizophrenia here... c c - This program was built off of the NADCON one, which made VECTOR plots. But with VERTCON3, I've decided c - to plot COLORED DOTS instead of B/W vectors. As such, the code which I called "lazy coding", where I spun through the work c - file JUST to compute statistics of the vectors so I could determine a "reference vector" length, c - needs to stay, but instead of determining a "reference vector length", I instead generate a "color palette" c - based on those stats. The question I ask is: since this program ALSO interrogates the GRIDS created FROM c - the point data to determine a color palette, couldn't I just use the same color palette? The answer I provide c - to myself is "maybe". I think there is danger in doing so, because gridded data may have some big deviations which c - occur away from the point data, so the color palettes could be different to each make a good map. On the other c - hand, there's a certain symmetry in identical color palettes, so that the point and grid data can be directly c - compared to one another. c - As of today (5/9/2018) I am going to try "identical" color palettes as follows: c - 1) For the "CD" point data and the grid of "CD" data (all cases) c - 2) For the "RZ" point data and the grid of "RZ" data (v2special only) c - 3) We won't do a "V2" point data plot, as that seems kind of stupid and unnecessary. However, I would c enforce the "V2" grid to be colored exactly as the "CD" grid so they can be direclty compared. cccccccccccccccccccccccccccccccccccccccccccccccccccccccc c - 2018 03 29: Cannibalized version of this code taken from NADCON 5 to use in VERTCON 3 c - Note, that for VERTCON, using just orthometric heights, I'm switching from vectors c - to colored dots. c c - Program presumes, as input c - the "old datum", "new datum" , "region" and "filter" and "gridsec" c - and "mapflag" arguments fed into "doit3.bat" c - Examples: c olddatum = 'ngvd29' c newdatum = 'navd88' c region = 'conus' c filter = '0' c gridsec = '900' c - The logical argument "v2special" means we are in NGVD29/NAVD88/CONUS. It means that c - a special remove/compute/restore process was used specfically for that triad which c - builds around the existing VERTCON 2 grid. In such a case, we have extra files c - rather than *just* the Coordinate Differnce ("cd") files. We also have "v2" (from c - the VERTCON 2 grid) and "rz" (Residuals between cd and v2 values). c c - Creates a batch file called "gmtbat03.(olddtm).(newdtm).(region).(filter).(igridsec) c - That batch file will create JPGs of: c c - 1)B/W plots of thinned coverage of points that went into the T=0.4 transformation grid c If NOT v2special: c -> These are just thinned "cd" datasets that were pulled at T=0.4 to make the transformation grid c If v2special: c -> Because "coverage" is identical for "rz" or "cd" data, we just stick with "cd" data as above c c - 2)B/W plots of dropped coverage of points that did not go into the T=0.4 transformation grid c If NOT v2special: c -> These are just thinned "cd" datasets that were pulled at T=0.4 to make the transformation grid c If v2special: c -> Because "coverage" is identical for "rz" or "cd" data, we just stick with "cd" data as above c c - 3)Color plots of thinned vectors (as colored dots) that went into the T=0.4 transformation grid c If NOT v2special: c -> These are just thinned "cd" datasets that were pulled at T=0.4 to make the transformation grid c If v2special: c -> These are the thinned "rz" datasets that were pulled at T=0.4 to *contribute* to the transformation grid c c - 4)Color plots of dropped vectors (as colored dots) that did not go into the T=0.4 transformation grid c If NOT v2special: c -> These are just dropped "cd" datasets that didn't go into the transformation grid c If v2special: c -> These are the dropped "rz" datasets that didn't contribute to transformation grid c c - 5)Color Plots of the doht grid at T=0.4 c If NOT v2special: c -> Just plot "cd" T=0.4 grid (actual gridded "cd" values pulled at T=0.4) c If v2special: c -> Make 3 plots: "cd" , "rz" and "v2" grid, BUT note this: c * The "rz" T=0.4 grid is actual residuals gridded at T=0.4 c * The "cd" T=0.4 grid is the sum of the "rz" T=0.4 grid (above) c and the Vertcon 2 original grid. It's NOT actually a gridding c of "cd" values at T=0.4! c * The "v2" T=0.4 grid is just a falsely named copy of the original VERTCON 2 grid c c - 6)Color Plots of the "method noise" grids (the "d3" grids, see DRU-11, p. 150) with thinned coverage overlaid c If NOT v2special: c -> Just plot "d3" grid from gridding "cd" values at 0.0 and 1.0 c If v2special: c -> Just plot "d3" grid from gridding "rz" values at 0.0 and 1.0. There is no "method noise" c associated with "cd" nor "v2" data. c --------------------------------------------------------------------- c - Below are comments from NADCON5, not terribly relevant to now, but c - left here for historic reference. c --------------------------------------------------------------------- c - 2016 08 26 c - Added new code to do reference vectors consistently c - See DRU-12, p. 56-57 c - Also fixed a typo in reference vector length for "vmtcdeht" 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 07 29: c Scrapped code about personalized reference vectors. Just put all c reference vectors outside/below plot c - Also moved "gridstats" and "vecstats" out into the /Subs directory c - to be used by other programs (like "makeplotfiles03.f") c - 2016 07 28: c Changed code to build the color palette of the "d3" grids c around the median, and not ave or std. c See DRU-12, p. 48 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 - Updated 01/21/2016 to get the CPT values fixed in d3 grids, so that c - (cpthi - cptlo) is exactly divisible by cptin at (2 x csm) c - Updated 10/27/2015 to work with the new naming scheme (see DRU-11, p. 150) c - Updated 10/05/2015 to work with the new naming scheme (see DRU-11, p. 139) c c - Part of the NADCON5 process c c - Built upon the skeleton of "makeplotem.f" for GEOCON v2.0 c - But built specifically for NADCON v5.0. So different c - in file names and expanded plot creation that it was c - given the new name "makeplotfiles02.f" to align with c - another NADCON5 program "makeplotfiles01.f" c --------------------------------------------------------------------- c - End of historic NADCON5 comments c --------------------------------------------------------------------- implicit real*8(a-h,o-z) parameter(maxplots=60) character*35 fname0 character*34 dirnam character*10 olddtm,newdtm,od,nd character*2 state,stdum character*3 ele,elelat,elelon,eleeht,elehor,ele0 character*5 agridsec character*5 agridsec2 integer*4 igridsec integer*4 igridsec2 character*200 suffix1,suffix2,suffix3,suffix2d3 character*200 suffix2t04 character*200 wfname,gmtfile character*200 gfncvtcdoht character*200 gfncvdcdoht character*200 bfnvmtcdoht character*200 bfnvmtrzoht character*200 bfnvmtv2oht c - 10/27/2015: "sad" = "Scaled AbsoluteValue Differential", aka "d3" grids c - See DRU-11, p. 150 character*200 sadbfnvmtcdoht character*200 gfnvmtcdoht character*200 gfnvmdcdoht character*200 sadbfnvmtrzoht character*200 gfnvmtrzoht character*200 gfnvmdrzoht character*200 sadbfnvmtv2oht character*200 gfnvmtv2oht character*200 gfnvmdv2oht 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 rejoht real*8 xlath,xlonh,xohth real*8 dohtm,v2dohtm,residual character*10 olddtm,newdtm,region c - Units to put on scale on GMT color plots character*17 scunlat,scunlon,scuneht,scunhor character*1 mapflag character*3 filter c - Plot/Scale stuff real*8 lorvopc real*8 lorvogohtm,lorvogohtmv2,lorvogohtmrz character*1 h1f1,h1f2,h1f3,h1f4,h2f1,h2f2,h2f3,h2f4 c - v2special (29/88/conus, remove/compute/restore stuff): logical v2special character*200 v2fnameg,dumfng character*200 v2fnameb,dumfnb c character*4 dumgridsec c character*4 dumgridsec2 c ------------------------------------------------------------------ c - BEGIN PROGRAM c ------------------------------------------------------------------ write(6,1001) 1001 format('BEGIN program makeplotfiles02.f') c-------------------------------------------- c - Some necessary constants c c - "lorvopc" = Length of Reference Vector on Paper in **CM** c - *** CRITICAL: The units of centimeters align with c - the use, in .gmtdefaults4, of: c - MEASURE_UNIT = cm c - If you change that value, then the units of lorvopc c - will be whatever you change it too, which may prove c - displeasing. c - c - "csm" is the Color Sigma Multiplier...the c - number of Standard Deviations away from c - the average for the color palette to c - span on color plots. c-------------------------------------------- lorvopc = 1.d0 csm = 3.d0 pi = 2.d0*dasin(1.d0) d2r = pi/180.d0 re = 6371000.d0 MultiplierLorvog = 2 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 - 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) suffix2t04=trim(suffix2)//'.04' suffix2d3=trim(suffix2)//'.d3' suffix3=trim(suffix2)//'.'//trim(mapflag) c goto 888 c ------------------------------------------------------------------ c - LAZY CODING: Spin through work file, collecting stats c - which we will use to scale the map and vectors, etc. c ------------------------------------------------------------------ c ------------------------------------------------------------------ c - Open the work file c ------------------------------------------------------------------ wfname='work.'//trim(suffix1) open(1,file='Work/'//wfname,status='old',form='formatted') write(6,1004)trim(wfname) 1004 format(6x,'makeplotfiles02.f: Opening work file ',a) ndoht = 0 avedohtm = 0.d0 rmsdohtm = 0.d0 stddohtm = 0.d0 ndohtv2 = 0 avedohtmv2 = 0.d0 rmsdohtmv2 = 0.d0 stddohtmv2 = 0.d0 ndohtrz = 0 avedohtmrz = 0.d0 rmsdohtmrz = 0.d0 stddohtmrz = 0.d0 891 if(.not.v2special)then read(1,104,end=892)pid,state,rejoht, * xlath,xlonh,xohth,dohtm,olddtm,newdtm, * h1f1,h1f2,h1f3,h1f4,h2f1,h2f2,h2f3,h2f4 else read(1,2104,end=892)pid,state,rejoht, * xlath,xlonh,xohth,dohtm,olddtm,newdtm, * h1f1,h1f2,h1f3,h1f4,h2f1,h2f2,h2f3,h2f4, * v2dohtm,residual endif if(rejoht.eq.' ')then c - Because they can be +/-, we and this is JUST for scaling the map, c - we want the average MAGNITUDE...use ABS. This means that c - I don't need to accept my own latter advice of returning c - here and using RMS. avedohtm = avedohtm + dabs(dohtm) rmsdohtm = rmsdohtm + dabs(dohtm**2) ndoht = ndoht + 1 if(v2special)then avedohtmv2 = avedohtmv2 + dabs(v2dohtm) rmsdohtmv2 = rmsdohtmv2 + dabs(v2dohtm**2) ndohtv2 = ndohtv2 + 1 avedohtmrz = avedohtmrz + dabs(residual) rmsdohtmrz = rmsdohtmrz + dabs(residual**2) ndohtrz = ndohtrz + 1 endif endif goto 891 892 continue c - Compute final stats if(ndoht.gt.0)then avedohtm = avedohtm / dble(ndoht) rmsdohtm = dsqrt(rmsdohtm / dble(ndoht)) fac = dble(ndoht) / dble(ndoht-1) stddohtm = dsqrt(fac*(rmsdohtm**2-avedohtm**2)) if(v2special)then avedohtmv2 = avedohtmv2 / dble(ndohtv2) rmsdohtmv2 = dsqrt(rmsdohtmv2 / dble(ndohtv2)) fac = dble(ndohtv2) / dble(ndohtv2-1) stddohtmv2 = dsqrt(fac*(rmsdohtmv2**2-avedohtmv2**2)) avedohtmrz = avedohtmrz / dble(ndohtv2) rmsdohtmrz = dsqrt(rmsdohtmrz / dble(ndohtrz)) fac = dble(ndohtrz) / dble(ndohtrz-1) stddohtmrz = dsqrt(fac*(rmsdohtmrz**2-avedohtmrz**2)) endif else avedohtm = 0.d0 stddohtm = 0.d0 if(v2special)then avedohtmv2 = 0.d0 stddohtmv2 = 0.d0 avedohtmrz = 0.d0 stddohtmrz = 0.d0 endif endif if(.not.v2special)then write(6,893)ndoht,avedohtm,stddohtm else write(6,2893)ndoht,avedohtm,stddohtm, * avedohtmv2,stddohtmv2, * avedohtmrz,stddohtmrz endif 893 format(6x,'makeplotfiles02.f: Vector Stats: ',/, *10x,'Number of Good Ort. Ht. Vectors : ',i10,/, *10x,'Average length (meters) : ',f10.3,/, *10x,'STD length (meters) : ',f10.3) 2893 format(6x,'makeplotfiles02.f: Vector Stats: ',/, *10x,'Number of Good Ort. Ht. Vectors : ',i10,/, *10x,'Average length (meters) : ',f10.3,/, *10x,'STD length (meters) : ',f10.3,/, *10x,'Average length, Vertcon2 (meters) : ',f10.3,/, *10x,'STD length, Vertcon2 (meters) : ',f10.3,/, *10x,'Average length, Residual (meters) : ',f10.3,/, *10x,'STD length, Residual (meters) : ',f10.3) 104 format(a6,1x,a2,a1,1x,f14.10,1x,f14.10,1x,f9.4,1x, *f9.4,1x,a10,1x,a10,1x,8(a1,1x)) 2104 format(a6,1x,a2,a1,1x,f14.10,1x,f14.10,1x,f9.4,1x, * f9.4,1x,a10,1x,a10,1x,8(a1,1x),1x,f9.4,1x,f9.4) c ------------------------------------------------------------------ c - END OF LAZY CODING c ------------------------------------------------------------------ 888 continue 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) write(6,1006)trim(region) 1006 format(6x,'makeplotfiles02.f: Calling getmapbounds for region ',a) c - Note, we do not NEED "lorv" (length of reference vector) computations c - at ALL, but this code is here as an artifact. It gets sent to c - "coplotvc.f" subroutine but never used) c --------------------------------------------------------------- c - Compute and report various things about our coverage c - and vector plots. c --------------------------------------------------------------- c - Reference vector set to length of the average c - absolute value of horizontal or ell ht shift. c - Note "gm2pc" = "ground meters to paper centimeters" c - "gs2pc" = "ground arcseconds to paper centimeters" c - "lorvog*s" = "length of reference vector on ground (element) in arcseconds" c - "lorvog*m" = "length of reference vector on ground (element) in meters" if(ndoht.ne.0)then lorvogohtm = onzd2(MultiplierLorvog*avedohtm) gm2pcoht = lorvopc / lorvogohtm c - VERTCON3, NGVD29/NAVD88/CONUS specials: if(v2special)then lorvogohtmv2 = onzd2(MultiplierLorvog*avedohtmv2) gm2pcohtv2 = lorvopc / lorvogohtmv2 lorvogohtmrz = onzd2(MultiplierLorvog*avedohtmrz) gm2pcohtrz = lorvopc / lorvogohtmrz endif endif c - Report: write(6,894)nplots 894 format(6x,'makeplotfiles02.f: Info about plots:',/, *8x,'Number of sub-region 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(ndoht .ne.0)then write(6,899)lorvogohtm,gm2pcoht if(v2special)then write(6,2899)lorvogohtmv2,gm2pcohtv2 write(6,3899)lorvogohtmrz,gm2pcohtrz endif else write(6,900) endif 895 continue 896 format( *8x,'Plot # ',i2,/, *10x,'S/N/W/E/N-S/E-W = ',6f7.1) 899 format( *10x,'Ort. Height plots: ',/, *12x,'Reference Vector ( m) = ',f10.2,/, *12x,'Ground M to Paper cm = ',f20.10) 2899 format( *10x,'Ort. Height plots (vertcon2 grid): ',/, *12x,'Reference Vector ( m) = ',f10.2,/, *12x,'Ground M to Paper cm = ',f20.10) 3899 format( *10x,'Ort. Height plots (residuals): ',/, *12x,'Reference Vector ( m) = ',f10.2,/, *12x,'Ground M to Paper cm = ',f20.10) 900 format( *10x,'Ort. Height plots: ',/, *12x,'No orthometric height 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 - The output batch file full of all the GMT commands c - needed to create the plots of interest c ------------------------------------------------------------------ gmtfile = 'gmtbat03.'//trim(suffix3) open(99,file=gmtfile,status='new',form='formatted') write(6,1011)trim(gmtfile) 1011 format(6x,'makeplotfiles02.f: Creating GMT batch file ',a) write(99,1030)trim(gmtfile) 1030 format('echo BEGIN batch file ',a,/, *'gmtset MAP_GRID_PEN_PRIMARY 0.25p,-',/, *'gmtset MAP_FRAME_TYPE fancy',/, *'gmtset FONT_TITLE 12p,Helvetica',/, *'gmtset MAP_TITLE_OFFSET 0.5c') 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 - c - OR...we could just use the statistics from the POINT data c - which we already collected... c ------------------------------------------------------------------ gfncvtcdoht = 'cvtcdoht.'//trim(suffix2) gfncvdcdoht = 'cvdcdoht.'//trim(suffix2) bfnvmtcdoht = 'vmtcdoht.'//trim(suffix2t04)//'.b' sadbfnvmtcdoht = 'vmtcdoht.'//trim(suffix2d3)//'.b' gfnvmtcdoht = 'vmtcdoht.'//trim(suffix2) gfnvmdcdoht = 'vmdcdoht.'//trim(suffix2) bfnvmtrzoht = 'vmtrzoht.'//trim(suffix2t04)//'.b' sadbfnvmtrzoht = 'vmtrzoht.'//trim(suffix2d3)//'.b' gfnvmtrzoht = 'vmtrzoht.'//trim(suffix2) gfnvmdrzoht = 'vmdrzoht.'//trim(suffix2) bfnvmtv2oht = 'vmtv2oht.'//trim(suffix2t04)//'.b' sadbfnvmtv2oht = 'vmtv2oht.'//trim(suffix2d3)//'.b' gfnvmtv2oht = 'vmtv2oht.'//trim(suffix2) gfnvmdv2oht = 'vmdv2oht.'//trim(suffix2) if(v2special)then call vecstats(gfnvmtrzoht,nthinoht) else call vecstats(gfnvmtcdoht,nthinoht) endif c------------------------------------- c - Stats of the TRANSFORMATION grids ("CD" if not v2special, "RZ" if v2special) c------------------------------------- if(nthinoht.ne.0)then c - First off, whether I'm in "v2special" or not, I need a color palette for the "CD" data c write(6,1012)trim(bfnvmtcdoht) c call gridstats(bfnvmtcdoht,ave,std,xmd) c aveohtm = ave c stdohtm = std write(6,*)' Calling cpt for CD values' c - Using grid stats: c call cpt(aveohtm,stdohtm,csm,cptloohtm,cpthiohtm,cptinohtm) c - Using point stats: c call cpt(aveohtm,stdohtm,csm,cptloohtm,cpthiohtm,cptinohtm) call cpt(avedohtm,stddohtm,csm,cptloohtm, * cpthiohtm,cptinohtm) c - Now, for v2special, I do this: c a) Set the color palette for the "rz" data c b) Later I'll use the "cd" color palette on "v2" data if(v2special)then c write(6,1012)trim(bfnvmtrzoht) c call gridstats(bfnvmtrzoht,ave,std,xmd) c aveohtmrz = ave c stdohtmrz = std write(6,*)' Calling cpt for gridded RZ values' c - Using grid stats: c call cpt(aveohtmrz,stdohtmrz,csm,cptloohtmrz, c * cpthiohtmrz,cptinohtmrz) c - Using point stats: call cpt(avedohtmrz,stddohtmrz,csm,cptloohtmrz, * cpthiohtmrz,cptinohtmrz) endif endif c------------------------------------- c - Stats of the METHOD NOISE ("d3") grids c - In this case, we definitely want the gridded data to c - rule the color palette. And we want that grid c - to be "cd" if no v2special, but "rz" if it is v2special c------------------------------------- if(nthinoht.ne.0)then if(v2special)then write(6,1012)trim(sadbfnvmtrzoht) call gridstats(sadbfnvmtrzoht,ave,std,xmd) aveohtmsadrz = ave stdohtmsadrz = std xmdohtmsadrz = xmd write(6,1042)trim(sadbfnvmtrzoht) 1042 format(8x,'makeplotfiles02.f: Calling cpt2 on file: ',a) call cpt2(xmdohtmsadrz,2.d0, * cptloohtmsadrz,cpthiohtmsadrz,cptinohtmsadrz) else write(6,1012)trim(sadbfnvmtcdoht) call gridstats(sadbfnvmtcdoht,ave,std,xmd) aveohtmsad = ave stdohtmsad = std xmdohtmsad = xmd write(6,1042)trim(sadbfnvmtcdoht) call cpt2(xmdohtmsad,2.d0, * cptloohtmsad,cpthiohtmsad,cptinohtmsad) endif endif 1012 format(6x,'makeplotfiles02.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 Thinned data c - Orthometric Height c - VERTCON3: Added "filter" to call. Also, coverage is ALWAYS c - *just* done on "cd" data. if(nthinoht.ne.0) * call bwplotcv('oht',gfncvtcdoht,bw,be,bs,bn,jm, * b1,b2,maxplots,olddtm,newdtm,region,filter, * 'OHT',ij, * igridsec,fn) c - Make B/W Coverage Plots of Dropped data c - Orthometric Height if(nthinoht.ne.0) * call bwplotcv('oht',gfncvdcdoht,bw,be,bs,bn,jm, * b1,b2,maxplots,olddtm,newdtm,region,filter, * 'OHT',ij, * igridsec,fn) c ------------------ c - VECTORS PLOTS (Not really, for VERTCON 3 -- we c - have switched to colored dots) c ------------------ 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-------------------------------- c - Plots of THINNED data c-------------------------------- c - Orthometric height if(nthinoht.ne.0)then c - Always do the Hnew-Hold ("cd") plot... c - But for "v2special" we do not "thin" or "drop" the "v2" vector c - files nor the true "cd" vector files, only the "rz" files. As such c - we do not PLOT thin/drop POINTS (vectors) for the "v2" nor "cd" files, c - but we will do the color GRIDS of them (later on in this program). 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 - If we are in "v2special" mode, then we plot RESIDUALS (rz) at this point, not COORD DIFFS (cd) c - (We'll get there eventually though) if(v2special)then c call cpt(avedohtmrz,stddohtmrz,csm,cptlorz,cpthirz,cptinrz) call coplotvc('oht',gfnvmtrzoht,bw,be,bs,bn, * jm,b1,b2,maxplots, * olddtm,newdtm,region,filter,'OHT',ij,xvlon, * xvlat,xllon,xllat, * lorvogohtmrz,lorvopc,igridsec,igridsec2,fn, * cptloohtmrz,cpthiohtmrz,cptinohtmrz) else c call cpt(avedohtm,stddohtm,csm,cptlocd,cpthicd,cptincd) call coplotvc('oht',gfnvmtcdoht,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 endif c-------------------------------- c - Plots of DROPPED data c-------------------------------- c - Orthometric Height if(nthinoht.ne.0)then c - Always do the Hnew-Hold ("cd") plot... ccccccccc c - 2018 05 04: FALSE c - In "v2special", however, we have put "rz" values into the c - variable called "cd" for ease of coding... ccccccccc c - Furthermore, we do not "thin" or "drop" the "v2" vector files nor the true "cd" vector files. As such c - we do not PLOT thin/drop for the "v2" nor "cd" files. c - If we are in "v2special" mode, then we plot RESIDUALS (rz), not COORD DIFFS (cd) if(v2special)then c call cpt(averzdohtm,stdrzdohtm,csm,cptlorz,cpthirz,cptinrz) call coplotvc('oht',gfnvmdrzoht,bw,be,bs,bn, * jm,b1,b2,maxplots, * olddtm,newdtm,region,filter,'OHT',ij,xvlon, * xvlat,xllon,xllat, * lorvogohtmrz,lorvopc,igridsec,igridsec2,fn, * cptloohtmrz,cpthiohtmrz,cptinohtmrz) else cccc c write(6,1041)trim(bfnvmtrzoht) c call cpt(avedohtm,stddohtm,csm,cptlocd,cpthicd,cptincd) 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") call coplotvc('oht',gfnvmdcdoht,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 endif c -------------------- c - Color gridded data (Transformation grid) c -------------------- c - Make color plots of gridded (T=0.4), thinned data, with no c - points. c - Orthometric Height c - VERTCON3: Added "filter" to the call. Also, updated subroutine c - to be GMT5 if(nthinoht.ne.0)then c - No matter what, we will plot a "cd" grid: call coplot('oht',bfnvmtcdoht,bw,be,bs,bn,jm,b1,b2,maxplots, * olddtm,newdtm,region,filter,'OHT', * ij,cptloohtm,cpthiohtm,cptinohtm, * suffix2t04,igridsec,igridsec2,fn) cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c - 2018 05 04 : FALSE now about faking "rz" data into "cd" variable names c c - Special code for 29/88/CONUS (v2special): For this particular c - combo, so far we've been thinning, dropping c - and gridding "rz" data data , though we've been carrying c - it forward in the "cd" file names, like "fvnvmtCDoht", even c - though it's "rz" data. Now here, at the end of the c - program, when we create the color plots of the gridded c - in "rz" data, we also want to create PARALLEL plots c - of the "v2" data and the true "cd" data. c - c - The "v2" data is pre-stored elsewhere: c /home/dru/VERTCON/VERTCON2/vertcon2.0.v3bounds.sec####.meters.grd c where #### = integer from 0060 (1 min) through 3600 (60 min) c - However, due to the rigid way I wrote "coplot", requiring specific c - file names, I can't just plot the above named files. Rather, c - I copy them into an alias. See DRU-12, p. 130. c c - c - The "cd" grid is created inside of "gmtbat02" calls (as generated c - by "mymedian5", so it will exist by the time we get here. c - Note: The "cd" grid (at T=0.4) will be one of two things, but c - it WILL exist either way: c If NOT v2special: c -> Will be the grid from pulling thinned "cd" data at 0.4 c If v2special: c -> Will be sum of "v2" data and grid from pulling thinned "rz" data at 0.4 c c Either way, it will be named: "vmtcdoht.ngvd29.navd88.conus...04.b" c c - Why "04"? Because the "rz" data were pulled at 0.4 tension, and this new c - "cd" grid is just a sum of that and the old VERTCON2 grid, so we maintain c - it, even though it was NOT generated by pulling "cd" data at T=0.4 c - c - So, let's replicate the above "coplot" calls, only using the file names c - for "v2" grid and "cd" grid: v2fnameg,dumfng cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c write(dumgridsec,'(i4.4)')igridsec if(v2special)then c - V2 data: c - Use the color palette from the "cd" data, so they can be compared c - easier! c call cpt(avev2dohtm,stdv2dohtm,csm, c * cptlov2dohtm,cpthiv2dohtm,cptinv2dohtm) c call coplot('oht',v2fnameb,bw,be,bs,bn,jm,b1,b2,maxplots, c * olddtm,newdtm,region,filter,'OHT', c * ij,cptlov2dohtm,cpthiv2dohtm,cptinv2dohtm, c * suffix2t04,igridsec,fn) call coplot('oht',bfnvmtv2oht,bw,be,bs,bn,jm,b1,b2,maxplots, * olddtm,newdtm,region,filter,'OHT', * ij,cptloohtm,cpthiohtm,cptinohtm, * suffix2t04,igridsec,igridsec2,fn) c - Finally, the "rz" data c - Need to get the color palette for it: c call cpt(avev2dohtm,stdv2dohtm,csm, c * cptlov2dohtm,cpthiv2dohtm,cptinv2dohtm) c call coplot('oht',dumfnb,bw,be,bs,bn,jm,b1,b2,maxplots, c * olddtm,newdtm,region,filter,'OHT', c * ij,cptlov2dohtm,cpthiv2dohtm,cptinv2dohtm, c * suffix2t04,igridsec,fn) call coplot('oht',bfnvmtrzoht,bw,be,bs,bn,jm,b1,b2,maxplots, * olddtm,newdtm,region,filter,'OHT', * ij,cptloohtmrz,cpthiohtmrz,cptinohtmrz, * suffix2t04,igridsec,igridsec2,fn) endif endif c -------------------- c - Color gridded data ("method noise", AKA "d3" grid) c -------------------- c - Make color plots of "method noise" *with* data overlain. See c - DRU-11, p.150 c - Orthometric Height c - VERTCON3: Added "filter" to the call. Also, updated subroutine c - to be GMT5 if(nthinoht.ne.0)then if(v2special)then call coplotwcv('oht',sadbfnvmtrzoht,bw,be,bs,bn,jm,b1,b2, * maxplots,olddtm,newdtm,region,filter, * 'OHT',ij,cptloohtmsadrz, * cpthiohtmsadrz,cptinohtmsadrz,suffix2d3,igridsec,igridsec2, * fn, * gfncvtcdoht) else call coplotwcv('oht',sadbfnvmtcdoht,bw,be,bs,bn,jm,b1,b2, * maxplots,olddtm,newdtm,region,filter, * 'OHT',ij,cptloohtmsad, * cpthiohtmsad,cptinohtmsad,suffix2d3,igridsec,igridsec2, * fn, * gfncvtcdoht) endif 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 makeplotfiles02.f') end c c ------------------------------------------------------ c c include 'Subs/bwplotvc.f' include 'Subs/coplotvc.f' include 'Subs/bwplotcv.f' include 'Subs/coplot.f' include 'Subs/coplotwcv.f' include 'Subs/plotcoast.f' include 'Subs/getmapbounds.f' include 'Subs/getmag.f' include 'Subs/cpt.f' c include '/home/dru/Subs/onzd2.f' ! don't need this cuz it's in "cpt" already include '/home/dru/NumRec/Modified/select2.for' include 'Subs/cpt2.f' include 'Subs/gridstats.f' include 'Subs/vecstats.f'