NADCON5-ng  0.0.1
NADCON5 Next Generation
coplotwcv.f
Go to the documentation of this file.
1 c> \ingroup core
2 c> Subroutine to make GMT calls to do a color raster rendering
3 c> of gridded data, with coverage overlaid.
4 c>
5 c> ## Changelog
6 c>
7 c> ### 2016 09 08:
8 c> Had to up the D_FORMAT default to %.3G because tight scalebar ranges
9 c> with the newly allowed "more free" average values were showing
10 c> repeating values when only 2 digits could be shown.
11 c>
12 c> ### 2016 09 07:
13 c> Had to add lines pre/post "makecpt" to change the D_FORMAT. This is because
14 c> I *had* been forcing the "scave" in cpt.f to be ONZD. But that yielded
15 c> bad values sometimes, so I switched it. With that switch, the scave could
16 c> have lots of digits. Well, that means the newly adopted "D_FORMAT" of %.2G
17 c> was insufficient for the CPT table. Who knew that D_FORMAT affected that!
18 c> Anyway, so change "D_FORMAT" pre/post all makecpt calls.
19 c>
20 c> ### 2016 08 30:
21 c> See item #39 in Google ToDo list
22 c> Changed "grdcontour" to have a blank "-R" call so it'll mimic whatever
23 c> decimal places are in the "grdimage" call that came before it.
24 c>
25 c> ### 2016 08 29:
26 c> Updated the initial -R and -B calls to 6 decimal places
27 c>
28 c> ### 2016 08 25:
29 c> - .gmtdefaults4 has been changed so X_ORIGIN is equal to 0.0
30 c> - Center the plot with "-Xc" at grdimage
31 c> - Center the scalebar by setting its Xcoordinate, which runs in "plot frame"
32 c> coordinates (0/0 at lower left) , to be equal to "jm/2"
33 c> - Force the scale bar to be exactly 4 inches wide, always
34 c> - Change the format for "makecpt" from 0.6 to 0.10
35 c>
36 c> ### 2016 07 29:
37 c> - Update to put more data into comment/echo
38 c> - Forced grdcontour with -Ctemp.cpt to align contours with color palette
39 c>
40 c> ### 2016 07 28:
41 c> - For d3 plots, to drop all contours
42 c> - For d3 plots to only use "coverage" part *only* if "ij" is not "1"
43 c> - Same for "09" and "ete" grids
44 c>
45 c> ### 2016 07 21:
46 c> - Set the "JM" code in "grdcontour" to just be "-JM" and let it therefore
47 c> run with whatever JM size is used in "grdimage"
48 c> - Set the "A" code in "grdcontour" to be "-A-" which should turn off
49 c> the labels on all contours
50 c> - Fixed the size of the scale bar
51 c>
52 c> ### 2016 03 01:
53 c> 1. Changed to a continuous color plot
54 c> - Get rid of "-Z" in "makecpt"
55 c> 2. Changed to an 8 color, from 6 color, plot (without changing the RANGE yet...see DRU-12, p. 19)
56 c> - Change varible that comes in from "cptin" to "cptin6"
57 c> - Compute cptin = cptin6 * 0.75d0 immediately
58 c>
59 c> ### 2016 02 29:
60 c> 1. Removed all shading from color plots
61 c> - Get rid of "grdgradient" call
62 c> - Remove from "grdimage" the "-Itempi.grd" part
63 c> - Remove the "rm -f tempi.grd" line
64 c>
65  subroutine coplotwcv(ele,fname,bw,be,bs,bn,jm,b1,b2,maxplots,
66  *olddtm,newdtm,region,elecap,ij,cptlo,cpthi,cptin6,suffixused,
67  *igridsec,fn,cvfname)
68 
69 c - 2016 09 08
70 c Had to up the D_FORMAT default to %.3G because tight scalebar ranges
71 c with the newly allowed "more free" average values were showing
72 c repeating values when only 2 digits could be shown.
73 
74 c - 2016 09 07
75 c Had to add lines pre/post "makecpt" to change the D_FORMAT. This is because
76 c I *had* been forcing the "scave" in cpt.f to be ONZD. But that yielded
77 c bad values sometimes, so I switched it. With that switch, the scave could
78 c have lots of digits. Well, that means the newly adopted "D_FORMAT" of %.2G
79 c was insufficient for the CPT table. Who knew that D_FORMAT affected that!
80 c Anyway, so change "D_FORMAT" pre/post all makecpt calls.
81 
82 c - 2016 08 30: See item #39 in Google ToDo list
83 c Changed "grdcontour" to have a blank "-R" call so it'll mimic whatever
84 c decimal places are in the "grdimage" call that came before it.
85 
86 c - 2016 08 29:
87 c Updated the initial -R and -B calls to 6 decimal places
88 
89 c - 2016 08 25:
90 c .gmtdefaults4 has been changed so X_ORIGIN is equal to 0.0
91 c Center the plot with "-Xc" at grdimage
92 c Center the scalebar by setting its Xcoordinate, which runs in "plot frame"
93 c coordinates (0/0 at lower left) , to be equal to "jm/2"
94 c Force the scale bar to be exactly 4 inches wide, always
95 c Change the format for "makecpt" from 0.6 to 0.10
96 
97 
98 c - 2016 07 29:
99 c Update to put more data into comment/echo
100 c - Forced grdcontour with -Ctemp.cpt to align contours with color palette
101 
102 c - Update 2016 07 28 to:
103 c For d3 plots, to drop all contours
104 c For d3 plots to only use "coverage" part *only* if "ij" is not "1"
105 c Same for "09" and "ete" grids
106 
107 c subroutine coplotwcv(ele,fname,bw,be,bs,bn,jm,b1,b2,maxplots,
108 c *olddtm,newdtm,region,elecap,ij,cptlo,cpthi,cptin,suffixused,
109 c *igridsec,fn,cvfname)
110 
111 c - Update 2016 07 21:
112 c - Set the "JM" code in "grdcontour" to just be "-JM" and let it therefore
113 c run with whatever JM size is used in "grdimage"
114 c - Set the "A" code in "grdcontour" to be "-A-" which should turn off
115 c the labels on all contours
116 c - Fixed the size of the scale bar
117 
118 
119 c - Update 2016 03 01:
120 c - 1) Changed to a continuous color plot
121 c * Get rid of "-Z" in "makecpt"
122 c 2) Changed to an 8 color, from 6 color, plot (without changing the RANGE yet...see DRU-12, p. 19)
123 c * Change varible that comes in from "cptin" to "cptin6"
124 c * Compute cptin = cptin6 * 0.75d0 immediately
125 
126 
127 c - Update 2016 02 29:
128 c - 1) Removed all shading from color plots
129 c * Get rid of "grdgradient" call
130 c * Remove from "grdimage" the "-Itempi.grd" part
131 c * Remove the "rm -f tempi.grd" line
132 c
133 
134  implicit real*8(a-h,o-z)
135 
136 c - Subroutine to make GMT calls to do a color raster rendering
137 c - of gridded data, with coverage overlaid.
138 
139  integer*4 maxplots
140  character*3 ele,elecap
141  character*200 fname,cvfname
142  character*10 olddtm,newdtm,region
143  real*8 bw(maxplots),be(maxplots),bs(maxplots),bn(maxplots)
144  real*4 jm(maxplots)
145  real*4 b1(maxplots),b2(maxplots)
146  character*10 fn(maxplots)
147  real*8 ave,std
148  character*200 suffixused
149  character*20 units
150  character*20 gridnote
151  character*20 gridnote2
152 
153 c ----------------------------------------------------
154 c - FAILSAFES: BEGIN
155 c ----------------------------------------------------
156  if(.not.(fname(1:3).eq.'vmt'.or.fname(1:3).eq.'vst' .or.
157  * fname(1:3).eq.'vmr'.or.fname(1:3).eq.'vsr'))then
158  write(6,2)trim(fname)
159  stop
160  endif
161  2 format('FATAL in coplotwcv. Bad character in spots 1-3: ',a)
162 
163  if(.not.(fname(4:5).eq.'cd' .or. fname(4:5).eq.'dd'))then
164  write(6,3)trim(fname)
165  stop
166  endif
167  3 format('FATAL in coplotwcv. Bad character in spots 4-5: ',a)
168 
169  if(fname(6:8).ne.ele)then
170  write(6,4)trim(fname),ele
171  stop
172  endif
173  4 format('FATAL in coplotwcv. Bad match of fname / ele: ',a,1x,a)
174 
175  if(fname(2:2).eq.'m')then
176  units='meters'
177  elseif(fname(2:2).eq.'s')then
178  units='arcseconds'
179  else
180  write(6,5)trim(fname)
181  stop
182  endif
183  5 format('FATAL in coplotwcv. Bad units in name: ',a)
184 
185  if(igridsec.le.0)then
186  write(6,6)igridsec
187  else
188  write(gridnote,10)igridsec
189  write(gridnote2,11)igridsec
190  endif
191  6 format('FATAL in coplotwcv. Bad igridsec: ',i0)
192  10 format('(',i0,' sec)')
193  11 format(i0,'.')
194 
195 c ----------------------------------------------------
196 c - FAILSAFES: END
197 c ----------------------------------------------------
198 
199 c ------------------------------------------------------------------
200 c - Contour interval and labeling
201 c ------------------------------------------------------------------
202 c - 2016 07 28: Scrapped the whole thing about scaling cptin.
203 c - Just bring in the interval for an 8 color spread, period.
204  cptin = cptin6
205 
206 c - 2016 03 01:
207 c cptin = cptin6 * 0.75d0
208 c - End 2016 03 01
209 
210  conin = cptin
211  conlb = conin / 2.d0
212 
213 c ----------------------------------------------
214 c - GMT COMMANDS: BEGIN
215 c ----------------------------------------------
216 
217 
218 c - Header of commands/echoes:
219 c - 2016 07 29:
220  write(99,991)ele,trim(units),trim(region),trim(fn(ij)),
221  *trim(suffixused),
222  *ele,trim(units),trim(region),trim(fn(ij)),
223  *trim(suffixused)
224 c write(99,991)ele,trim(units),trim(region),trim(fn(ij)),
225 c *ele,trim(units),trim(region),trim(fn(ij))
226 
227 
228 c - GMT call to create the color palette:
229 c - Store the palette as file "temp.cpt", to be
230 c - deleted later.
231  write(99,901)cptlo,cpthi,cptin
232 
233 c - Begin 2016 03 01:
234 c 901 format(
235 c *'makecpt -Crainbow -T',f0.6,'/',f0.6,'/',f0.6,
236 c *' -Z > temp.cpt')
237 c 901 format(
238 c *'makecpt -Crainbow -T',f0.6,'/',f0.6,'/',f0.6,
239 c *' > temp.cpt')
240 c - 2016 08 25
241 c 901 format(
242 c *'makecpt -Crainbow -T',f0.10,'/',f0.10,'/',f0.10,
243 c *' > temp.cpt')
244 c - 2016 09 07
245 c 901 format(
246 c *'gmtset D_FORMAT %.12G',/,
247 c *'makecpt -Crainbow -T',f0.10,'/',f0.10,'/',f0.10,
248 c *' > temp.cpt',/,
249 c *'gmtset D_FORMAT %.2G')
250 c - 2016 09 07
251  901 format(
252  *'gmtset D_FORMAT %.12G',/,
253  *'makecpt -Crainbow -T',f0.10,'/',f0.10,'/',f0.10,
254  *' > temp.cpt',/,
255  *'gmtset D_FORMAT %.3G')
256 
257 
258 c - GMT call to create the gradient file from
259 c - the "grd" formatted version of our input grid.
260 c - Store the gradients as file "tempi.grd" to be
261 c - deleted later.
262 
263 c - Updated 2016 02 29:
264 c write(99,902)fname(1:8),trim(suffixused)
265  902 format(
266  *'grdgradient ',a,'.',a,
267  *'.grd -N4.8 -A90 -M -Gtempi.grd')
268 c - Until here 2016 02 29
269 
270 c - GMT call to create the color render of the
271 c - gridded data. Store as file "plot.ps"
272 c - to be deleted later.
273 
274  write(99,903)fname(1:8),trim(suffixused),
275  *bw(ij),be(ij),bs(ij),bn(ij),
276  *jm(ij),b1(ij),b2(ij),
277  *trim(newdtm),trim(olddtm),elecap,igridsec,
278  *trim(region),trim(fn(ij)),'method noise'
279 
280 
281 c - Commented out for 2016 02 29:
282 c 903 format(
283 c *'grdimage ',a,'.',a,'.grd',
284 c *' -Ei -Itempi.grd',
285 c *' -R',f0.2,'/',f0.2,'/',f0.2,'/',f0.2,
286 c *' -JM',f3.1,'i -B',f0.2,'/',f0.2,
287 c *':."NADCON v5.0 ',a,' minus ',a,' ',a3,
288 c *'(',i0,' sec)',
289 c *1x,a,'-',a,1x,a,
290 c *':" -Ctemp.cpt -K > plot.ps')
291 c - Commented out for 2016 03 01:
292 c 903 format(
293 c *'grdimage ',a,'.',a,'.grd',
294 c *' -Ei ',
295 c *' -R',f0.2,'/',f0.2,'/',f0.2,'/',f0.2,
296 c *' -JM',f3.1,'i -B',f0.2,'/',f0.2,
297 c *':."NADCON v5.0 ',a,' minus ',a,' ',a3,
298 c *'(',i0,' sec)',
299 c *1x,a,'-',a,1x,a,
300 c *':" -Ctemp.cpt -K > plot.ps')
301 c - Below for 2016 08 25:
302 c 903 format(
303 c *'grdimage ',a,'.',a,'.grd',
304 c *' -Ei -Xc',
305 c *' -R',f0.2,'/',f0.2,'/',f0.2,'/',f0.2,
306 c *' -JM',f3.1,'i -B',f0.2,'/',f0.2,
307 c *':."NADCON v5.0 ',a,' minus ',a,' ',a3,
308 c *'(',i0,' sec)',
309 c *1x,a,'-',a,1x,a,
310 c *':" -Ctemp.cpt -K > plot.ps')
311 c - 2016 08 29:
312  903 format(
313  *'grdimage ',a,'.',a,'.grd',
314  *' -Ei -Xc',
315  *' -R',f0.6,'/',f0.6,'/',f0.6,'/',f0.6,
316  *' -JM',f3.1,'i -B',f0.6,'/',f0.6,
317  *':."NADCON v5.0 ',a,' minus ',a,' ',a3,
318  *'(',i0,' sec)',
319  *1x,a,'-',a,1x,a,
320  *':" -Ctemp.cpt -K > plot.ps')
321 
322 c - GMT call to put a color scale on the map. Add it
323 c - to file "plot.ps".
324 
325 c - 2016 07 21:
326 c qqa = jm(ij) - 1.d0
327 c qqb = qqa / 2.d0
328 c write(99,904)qqb,qqa,trim(units)
329 c - 2016 08 25
330  qqa = jm(ij)/2
331  write(99,904)qqa,trim(units)
332 
333  904 format(
334  *'psscale -D',f3.1,'i/-0.4i/4.0i/0.2ih',
335  *' -Ctemp.cpt -I0.4 -B/:',a,
336  *': -O -K >> plot.ps')
337 
338 c write(99,904)trim(units)
339 c 904 format(
340 c *'psscale -D4i/-0.4i/8i/0.2ih -Ctemp.cpt -I0.4 -B/:',a,
341 c *': -O -K >> plot.ps')
342 
343 c - GMT call to put contours on the map.
344 
345 c - 2016 07 28...skip contours if this is a d3, 09 or "ete" plot
346  ll = len(trim(suffixused))
347  if(suffixused(ll-1:ll).eq.'d3')goto 807
348  if(suffixused(ll-1:ll).eq.'09')goto 807
349  if(fname(3:5).eq.'ete')goto 807
350 
351 c - 2016 08 30:
352  write(99,908)fname(1:8),trim(suffixused)
353 c - 2016 07 29:
354 c write(99,908)fname(1:8),trim(suffixused),
355 c *bw(ij),be(ij),bs(ij),bn(ij)
356 c - 2016 07 21:
357 c write(99,908)fname(1:8),trim(suffixused),conin,
358 c *bw(ij),be(ij),bs(ij),bn(ij)
359 c write(99,908)fname(1:8),trim(suffixused),conin,
360 c *bw(ij),be(ij),bs(ij),bn(ij),
361 c *jm(ij),conlb
362 
363 c - 2016 08 30
364  908 format(
365  *'grdcontour ',a,'.',a,'.grd',
366  *' -Ctemp.cpt',
367  *' -R',
368  *' -JM -Wthin',
369  *' -A- -O -K >> plot.ps')
370 c - 2016 07 29:
371 c 908 format(
372 c *'grdcontour ',a,'.',a,'.grd',
373 c *' -Ctemp.cpt',
374 c *' -R',f0.2,'/',f0.2,'/',f0.2,'/',f0.2,
375 c *' -JM -Wthin',
376 c *' -A- -O -K >> plot.ps')
377 c - 2016 07 21:
378 c 908 format(
379 c *'grdcontour ',a,'.',a,'.grd',
380 c *' -C',f0.5,
381 c *' -R',f0.2,'/',f0.2,'/',f0.2,'/',f0.2,
382 c *' -JM -Wthin',
383 c *' -A- -O -K >> plot.ps')
384 c 908 format(
385 c *'grdcontour ',a,'.',a,'.grd',
386 c *' -C',f0.5,
387 c *' -R',f0.2,'/',f0.2,'/',f0.2,'/',f0.2,
388 c *' -JM',f3.1,'i -Wthin',
389 c *' -A',f0.5,' -O -K >> plot.ps')
390 
391 c - 2016 07 28 - Skip to here (skipping contours) for d3, 09 or "ete" plots
392  807 continue
393 
394 c - 2016 07 28 - Skip the coverage if this is a "d3" plot, 09 plot or "ete" plot *and*
395 c - we are at the overall map level of 1 ("entire")
396  if(suffixused(ll-1:ll).eq.'d3' .and. ij.eq.1)goto 808
397  if(suffixused(ll-1:ll).eq.'09' .and. ij.eq.1)goto 808
398  if(fname(3:5).eq.'ete' .and. ij.eq.1)goto 808
399 
400 c - GMT call to overlay the coverage
401  write(99,909)trim(cvfname)
402 
403  909 format(
404  *'psxy ',a,' -R -J -Sc0.02i -Gblack -O -K >> plot.ps')
405 
406  808 continue
407 
408 c - GMT call to create the shoreline. Add it to
409 c - file "plot.ps".
410 c - ALWAYS plot the coast last, as it closes the PS file
411  call plotcoast(region,99)
412 
413 c - GMT call to convert the plot.ps file to plot.jpg.
414  write(99,905)
415  905 format(
416  *'ps2raster plot.ps -Tj -P -A')
417 
418 c - Change name of JPG to its final name:
419  write(99,906)fname(2:8),trim(suffixused),trim(fn(ij))
420  906 format('mv plot.jpg c',a,'.',a,'.',a,'.jpg')
421 
422 c - Remove all the temporary files
423  write(99,907)
424 
425 c - Updated 2016 02 29:
426 c 907 format(
427 c *'rm -f temp.cpt',/,
428 c *'rm -f tempi.grd',/,
429 c *'rm -f plot.ps')
430  907 format(
431  *'rm -f temp.cpt',/,
432  *'rm -f plot.ps')
433 c - Until here 2016 02 29
434 
435 
436  991 format(
437  *'# -----------------------------------------------------',/,
438  *'# color plots with coverage in ',a,1x,a,1x,a,1x,a,
439  * 1x,' Grid:',a/,
440  *'# -----------------------------------------------------',/,
441  *'echo ...color plots with coverage in ',a,1x,a,1x,a,1x,a,
442  * 1x,' Grid:',a)
443 c 991 format(
444 c *'# -----------------------------------------------------',/,
445 c *'# color plots with coverage in ',a,1x,a,1x,a,1x,a,/,
446 c *'# -----------------------------------------------------',/,
447 c *'echo ...color plots with coverage in ',a,1x,a,1x,a,1x,a)
448 
449  return
450  end
subroutine coplotwcv(ele, fname, bw, be, bs, bn, jm, b1, b2, maxplots, olddtm, newdtm, region, elecap, ij, cptlo, cpthi, cptin6, suffixused, igridsec, fn, cvfname)
Subroutine to make GMT calls to do a color raster rendering of gridded data, with coverage overlaid...
Definition: coplotwcv.f:68
subroutine plotcoast(region, ifnum)
Subroutine to write GMT-based commands to create a shoreline Write GMT-based commands to create a sho...
Definition: plotcoast.f:40