subroutine coplot(ele,fname,bw,be,bs,bn,jm,b1,b2,maxplots, *olddtm,newdtm,region,filter,elecap, *ij,cptlo,cpthi,cptin6,suffixused, *igridsec,igridsec2,fn) c - Updated 2018 05 16 for 2 grid spacings c - 2018 03 30: Modified for use in VERTCON 3 with: c a) "filter" now comes in as something to be identified in the title c b) All commands are now GMT 5 compatable c - 2016 09 08 c Had to up the D_FORMAT default to %.3G because tight scalebar ranges c with the newly allowed "more free" average values were showing c repeating values when only 2 digits could be shown. c - 2016 09 07 c Had to add lines pre/post "makecpt" to change the D_FORMAT. This is because c I *had* been forcing the "scave" in cpt.f to be ONZD. But that yielded c bad values sometimes, so I switched it. With that switch, the scave could c have lots of digits. Well, that means the newly adopted "D_FORMAT" of %.2G c was insufficient for the CPT table. Who knew that D_FORMAT affected that! c Anyway, so change "D_FORMAT" pre/post all makecpt calls. c - 2016 08 30: See item #39 in Google ToDo list c Changed "grdcontour" to have a blank "-R" call so it'll mimic whatever c decimal places are in the "grdimage" call that came before it. c - 2016 08 29: c Updated the initial -R and -B calls to 6 decimal places c - 2016 08 25: c .gmtdefaults4 has been changed so X_ORIGIN is equal to 0.0 c Center the plot with "-Xc" at grdimage c Center the scalebar by setting its Xcoordinate, which runs in "plot frame" c coordinates (0/0 at lower left) , to be equal to "jm/2" c Force the scale bar to be exactly 4 inches wide, always c Change the format for "makecpt" from 0.6 to 0.10 c - 2016 07 29: c Update to put more data into comment/echo c Also forced the -Ctemp.cpt option in "grdcontour" to make the contours line up with the color palette c For d3 plots, to drop all contours c For d3 plots to only use "coverage" part *only* if "ij" is not "1" c Same for "09" and "ete" grids c subroutine coplot(ele,fname,bw,be,bs,bn,jm,b1,b2,maxplots, c *olddtm,newdtm,region,elecap,ij,cptlo,cpthi,cptin,suffixused, c *igridsec,fn) c - Update 2016 07 21: c - Set the "JM" code in "grdcontour" to just be "-JM" and let it therefore c run with whatever JM size is used in "grdimage" c - Set the "A" code in "grdcontour" to be "-A-" which should turn off c the labels on all contours c - Fixed size of scale bar c - Update 2016 03 01: c - 1) Changed to a continuous color plot c * Get rid of "-Z" in "makecpt" c 2) Changed to an 8 color, from 6 color, plot (without changing the RANGE yet...see DRU-12, p. 19) c * Change varible that comes in from "cptin" to "cptin6" c * Compute cptin = cptin6 * 0.75d0 immediately c - Update 2016 02 29: c - 1) Removed all shading from color plots c * Get rid of "grdgradient" call c * Remove from "grdimage" the "-Itempi.grd" part c * Remove the "rm -f tempi.grd" line c subroutine coplot(ele,fname,bw,be,bs,bn,jm,b1,b2,maxplots, c *olddtm,newdtm,region,elecap,ij,cptlo,cpthi,cptin,suffix2, c *igridsec,fn) c subroutine coplot(ele,fname,bw,be,bs,bn,jm,b1,b2,maxplots, c *olddtm,newdtm,region,elecap,ij,ave,std,suffix2,igridsec,fn) implicit real*8(a-h,o-z) c - Subroutine to make GMT calls to do a color raster rendering c - of gridded data integer*4 maxplots character*3 ele,elecap character*200 fname character*10 olddtm,newdtm,region character*3 filter real*8 bw(maxplots),be(maxplots),bs(maxplots),bn(maxplots) real*4 jm(maxplots) real*4 b1(maxplots),b2(maxplots) character*10 fn(maxplots) real*8 ave,std character*200 suffixused character*20 units character*23 gridnote logical v2special v2special = .false. if(trim(olddtm).eq.'ngvd29' .and. * trim(newdtm).eq.'navd88' .and. * trim(region).eq.'conus')v2special = .true. c------------------------------------------------------------------- c - New for VERTCON3: c Elements 1-2: No longer have a "vs" c Elements 4-5: Can have "cd" or "dd" still, but new codes for NGVD29/NAVD88/CONUS are: c v2 = Interpolated H88-H29 from VERTCON 2 grid c rz = Diff btwen "cd" and "v2" c Elements 6-8: Type "oht" is added as the only one allowed c------------------------------------------------------------------- c ---------------------------------------------------- c - FAILSAFES: BEGIN c ---------------------------------------------------- if(.not.(fname(1:3).eq.'vmt'.or.fname(1:3).eq.'vst' .or. * fname(1:3).eq.'vmr'.or.fname(1:3).eq.'vsr' .or. * fname(1:3).eq.'vme'.or.fname(1:3).eq.'vse'))then write(6,2)trim(fname) stop endif 2 format('FATAL in coplot. Bad character in spots 1-3: ',a) if(.not.(fname(4:5).eq.'cd' .or. fname(4:5).eq.'dd' .or. * fname(4:5).eq.'rz' .or. fname(4:5).eq.'v2' .or. * fname(4:5).eq.'te' ))then write(6,3)trim(fname) stop endif 3 format('FATAL in coplot. Bad character in spots 4-5: ',a) if(fname(6:8).ne.ele)then write(6,4)trim(fname),ele stop endif 4 format('FATAL in coplot. Bad match of fname / ele: ',a,1x,a) if(fname(2:2).eq.'m')then units='meters' elseif(fname(2:2).eq.'s')then units='arcseconds' else write(6,5)trim(fname) stop endif 5 format('FATAL in coplot. Bad units in name: ',a) if(igridsec.le.0)then write(6,6)igridsec else if(v2special)then if (fname(4:5).eq.'cd')then write(gridnote,11)igridsec,igridsec2 elseif(fname(4:5).eq.'v2')then write(gridnote,10)igridsec2 else write(gridnote,10)igridsec endif else write(gridnote,10)igridsec endif endif 6 format('FATAL in coplot. Bad igridsec: ',i0) 10 format('(',i0,' sec)') 11 format('(rz ',i0,'/v2 ',i0,' sec)') c ---------------------------------------------------- c - FAILSAFES: END c ---------------------------------------------------- c ------------------------------------------------------------------ c - Contour interval and labeling c ------------------------------------------------------------------ c - 2016 07 29: Just use cptin6 as it came in cptin = cptin6 c - Begin 2016 03 01 c cptin = cptin6 * 0.75d0 c - End 2016 03 01 conin = cptin conlb = conin / 2.d0 c ---------------------------------------------- c - GMT COMMANDS: BEGIN c ---------------------------------------------- c -------- HEADER: BEGIN ------------------------------------------- c - Header of commands/echoes: c - 2016 07 29: write(99,991)ele,trim(units),trim(region),trim(fn(ij)), *trim(suffixused), *ele,trim(units),trim(region),trim(fn(ij)), *trim(suffixused) c -------- HEADER: END ------------------------------------------- c -------- COLOR PALETTE: BEGIN ------------------------------------- c - GMT call to create the color palette: c - Store the palette as file "temp.cpt", to be c - deleted later. write(99,901)cptlo,cpthi,cptin c - 2018 03 30: For GMT5: 901 format( *'gmtset FORMAT_FLOAT_OUT %.12G',/, *'makecpt -Crainbow -T',f0.10,'/',f0.10,'/',f0.10, *' > temp.cpt',/, *'gmtset FORMAT_FLOAT_OUT %.3G') c - 2016 09 08 c 901 format( c *'gmtset D_FORMAT %.12G',/, c *'makecpt -Crainbow -T',f0.10,'/',f0.10,'/',f0.10, c *' > temp.cpt',/, c *'gmtset D_FORMAT %.3G') c - GMT call to create the gradient file from c - the "grd" formatted version of our input grid. c - Store the gradients as file "tempi.grd" to be c - deleted later. c -------- COLOR PALETTE: END ------------------------------------- c -------- GRADIENTS for SHADING: BEGIN ------------------------------------- c - Although "grdgradient" was dropped long ago, if we use it again, I've updated format 902 for GMT5: c - Updated 2016 02 29: c write(99,902)fname(1:8),trim(suffixused) 902 format( *'grdgradient ',a,'.',a, *'.grd -N4.8 -A90 -fg -Gtempi.grd') c 902 format( c *'grdgradient ',a,'.',a, c *'.grd -N4.8 -A90 -M -Gtempi.grd') c - Until here 2016 02 29 c -------- GRADIENTS for SHADING: END ------------------------------------- c -------- GRDIMAGE: BEGIN ------------------------------------- c - 2018 03 30: VERTCON3: Added "filter" to title in "grdimage" call c - 2018 03 30: GMT5: As far as I can tell, the "grdimage" all is same from GMT4 to GMT5 c - GMT call to create the color render of the c - gridded data. Store as file "plot.ps" c - to be deleted later. write(99,903)fname(1:8),trim(suffixused), *bw(ij),be(ij),bs(ij),bn(ij), *jm(ij),b1(ij),b2(ij), *trim(newdtm),trim(olddtm),elecap,trim(gridnote), *trim(region),trim(fn(ij)),fname(2:5),filter c write(99,903)fname(1:8),trim(suffixused), c *bw(ij),be(ij),bs(ij),bn(ij), c *jm(ij),b1(ij),b2(ij), c *trim(newdtm),trim(olddtm),elecap,igridsec, c *trim(region),trim(fn(ij)),fname(2:5),filter c write(99,903)fname(1:8),trim(suffixused), c *bw(ij),be(ij),bs(ij),bn(ij), c *jm(ij),b1(ij),b2(ij), c *trim(newdtm),trim(olddtm),elecap,igridsec, c *trim(region),trim(fn(ij)),fname(2:5) c - Commented this out for 2016 03 01: c 903 format( c *'grdimage ',a,'.',a,'.grd', c *' -Ei -Itempi.grd', c *' -R',f0.2,'/',f0.2,'/',f0.2,'/',f0.2, c *' -JM',f3.1,'i -B',f0.2,'/',f0.2, c *':."NADCON v5.0 ',a,' minus ',a,' ',a3, c *'(',i0,' sec)', c *1x,a,'-',a,1x,a, c *':" -Ctemp.cpt -K > plot.ps') c 903 format( c *'grdimage ',a,'.',a,'.grd', c *' -Ei ', c *' -R',f0.2,'/',f0.2,'/',f0.2,'/',f0.2, c *' -JM',f3.1,'i -B',f0.2,'/',f0.2, c *':."NADCON v5.0 ',a,' minus ',a,' ',a3, c *'(',i0,' sec)', c *1x,a,'-',a,1x,a, c *':" -Ctemp.cpt -K > plot.ps') c - Below for 2016 08 25: c 903 format( c *'grdimage ',a,'.',a,'.grd', c *' -Ei -Xc', c *' -R',f0.2,'/',f0.2,'/',f0.2,'/',f0.2, c *' -JM',f3.1,'i -B',f0.2,'/',f0.2, c *':."NADCON v5.0 ',a,' minus ',a,' ',a3, c *'(',i0,' sec)', c *1x,a,'-',a,1x,a, c *':" -Ctemp.cpt -K > plot.ps') c - 2016 08 29: c 903 format( c *'grdimage ',a,'.',a,'.grd', c *' -Ei -Xc', c *' -R',f0.6,'/',f0.6,'/',f0.6,'/',f0.6, c *' -JM',f3.1,'i -B',f0.6,'/',f0.6, c *':."NADCON v5.0 ',a,' minus ',a,' ',a3, c *'(',i0,' sec)', c *1x,a,'-',a,1x,a, c *':" -Ctemp.cpt -K > plot.ps') c 903 format( c *'grdimage ',a,'.',a,'.grd', c *' -Ei -Xc', c *' -R',f0.6,'/',f0.6,'/',f0.6,'/',f0.6, c *' -JM',f3.1,'i -B',f0.6,'/',f0.6, c *':."VERTCON v3.0 ',a,' minus ',a,' ',a3, c *'(',i0,' sec)', c *1x,a,'-',a,1x,a,'-f',a3, c *':" -Ctemp.cpt -K > plot.ps') 903 format( *'grdimage ',a,'.',a,'.grd', *' -Ei -Xc', *' -R',f0.6,'/',f0.6,'/',f0.6,'/',f0.6, *' -JM',f3.1,'i -B',f0.6,'/',f0.6, *':."VERTCON v3.0 ',a,' minus ',a,' ',a3, *a, *1x,a,'-',a,1x,a,'-f',a3, *':" -Ctemp.cpt -K > plot.ps') c -------- GRDIMAGE: END ------------------------------------- c -------- SCALE BAR: BEGIN ----------------------------------- c - 2018 03 30: GMT5 - psscale call changed c - GMT call to put a color scale on the map. Add it c - to file "plot.ps". c - This version, without the "-Ba" option, will label the scale bar c - exactly as per the CPT, but the problem is that sometimes c - the "0" in the CPT is really like "10e-16", which then throws c - 16 digits on everything on the scale bar and you can't read it. c - This is a problem inside of "makecpt" and I couldn't solve it. c - As such, I have the next version of 904 which uses "-Ba" to c - automatically pick labels in psscale. c - 2016 07 21 c write(99,904)trim(units) c qqa = jm(ij) - 1.d0 c qqb = qqa / 2.d0 c write(99,904)qqb,qqa,trim(units) c - 2016 08 25 qqa = jm(ij)/2 write(99,904)qqa,trim(units) 904 format( *'psscale -D',f3.1,'i/-0.4i/4.0i/0.2ih', *' -Ctemp.cpt -I0.4 -Ba -By+l',a, *': -O -K >> plot.ps') c 904 format( c *'psscale -D',f3.1,'i/-0.4i/4.0i/0.2ih', c *' -Ctemp.cpt -I0.4 -B/:',a, c *': -O -K >> plot.ps') c 904 format( c *'psscale -D4i/-0.4i/8i/0.2ih -Ctemp.cpt -I0.4 -B/:',a, c *': -O -K >> plot.ps') c -------- SCALE BAR: END ----------------------------------- c -------- CONTOURS: BEGIN ------------------------------------------- c - 2018 03 30: GMT5 - As far as I can tell, the "grdcontour" call is same as GMT4 c - 2018 03 30: VERTCON3: Skip contours entirely. Too local/noisy. c - GMT call to put contours on the map. c - 2016 07 29 c - 2016 07 21 goto 807 c 2016 07 29: skip contours if this is a d3, 09 or "ete" plot ll = len(trim(suffixused)) if(suffixused(ll-1:ll).eq.'d3')goto 807 if(suffixused(ll-1:ll).eq.'09')goto 807 if(fname(3:5).eq.'ete')goto 807 write(99,908)fname(1:8),trim(suffixused) c write(99,908)fname(1:8),trim(suffixused), c *bw(ij),be(ij),bs(ij),bn(ij) c write(99,908)fname(1:8),trim(suffixused),conin, c *bw(ij),be(ij),bs(ij),bn(ij) c write(99,908)fname(1:8),trim(suffixused),conin, c *bw(ij),be(ij),bs(ij),bn(ij), c *jm(ij),conlb c - 2016 08 30: 908 format( *'grdcontour ',a,'.',a,'.grd', *' -Ctemp.cpt', *' -R', *' -JM -Wthin', *' -A- -O -K >> plot.ps') c - 2016 07 29 c 908 format( c *'grdcontour ',a,'.',a,'.grd', c *' -Ctemp.cpt', c *' -R',f0.2,'/',f0.2,'/',f0.2,'/',f0.2, c *' -JM -Wthin', c *' -A- -O -K >> plot.ps') c - 2016 07 21 c 908 format( c *'grdcontour ',a,'.',a,'.grd', c *' -C',f0.5, c *' -R',f0.2,'/',f0.2,'/',f0.2,'/',f0.2, c *' -JM -Wthin', c *' -A- -O -K >> plot.ps') c 908 format( c *'grdcontour ',a,'.',a,'.grd', c *' -C',f0.5, c *' -R',f0.2,'/',f0.2,'/',f0.2,'/',f0.2, c *' -JM',f3.1,'i -Wthin', c *' -A',f0.5,' -O -K >> plot.ps') c -------- CONTOURS: END ------------------------------------------- c - 2016 07 28 - Skip to here (skipping contours) for d3, 09 or "ete" plots 807 continue c -------- SHORELINE: BEGIN ------------------------------------------- c - 2018 03 30: VERTCON3 - Updated "plotcoast" subroutine for new regions c - GMT call to create the shoreline. Add it to c - file "plot.ps". c - ALWAYS plot the coast last, as it closes the PS file call plotcoast(region,99) c -------- SHORELINE: END ------------------------------------------- c -------- CONVERT to JPG: BEGIN ------------------------------------------- c - 2018 03 30: GMT5 - "ps2raster" becomes "psconvert" c - GMT call to convert the plot.ps file to plot.jpg. c - 2018 03 30: VERTCON3 - Updated "plotcoast" subroutine for new regions write(99,905) 905 format( *'psconvert plot.ps -Tj -P -A') c 905 format( c *'ps2raster plot.ps -Tj -P -A') c -------- CONVERT to JPG: END ------------------------------------------- c - Change name of JPG to its final name: write(99,906)fname(2:8),trim(suffixused),trim(fn(ij)) 906 format('mv plot.jpg c',a,'.',a,'.',a,'.jpg') c - Remove all the temporary files write(99,907) c - Update 2016 02 29: c 907 format( c *'rm -f temp.cpt',/, c *'rm -f tempi.grd',/, c *'rm -f plot.ps') 907 format( *'rm -f temp.cpt',/, *'rm -f plot.ps') c - Until here 2016 02 29 c - 2016 07 29: 991 format( *'# -----------------------------------------------------',/, *'# color plots in ',a,1x,a,1x,a,1x,a,1x,' Grid:',a /, *'# -----------------------------------------------------',/, *'echo ...color plots in ',a,1x,a,1x,a,1x,a,1x,' Grid:',a) c 991 format( c *'# -----------------------------------------------------',/, c *'# color plots in ',a,1x,a,1x,a,1x,a,/, c *'# -----------------------------------------------------',/, c *'echo ...color plots in ',a,1x,a,1x,a,1x,a) return end