NADCON5 Core Library

Programs and Functions which are responsible for computing the grid transformations used to build NADCON5. More...

Programs, Functions, and Subroutines
program	addem
	Part of the NADCON5 NADCON5 Core Library , adds one grid to another. More...
program	b2xyz
	Part of the NADCON5 NADCON5 Core Library , converts *.b grid to xyz More...
program	convlv
	Part of the NADCON5 NADCON5 Core Library , Convolves two grids. More...
program	decimate
	Part of the NADCON5 NADCON5 Core Library , Extract a reduced (1 of n) dataset. More...
program	gabs
	Part of the NADCON5 NADCON5 Core Library , Convert values in a *.b grid to absolute value. More...
program	gscale
	Part of the NADCON5 NADCON5 Core Library , Scales a grid by a factor. More...
program	gsqr
	Part of the NADCON5 NADCON5 Core Library , Squares values in a *.b grid. More...
program	gsqrt
	Part of the NADCON5 NADCON5 Core Library , Square Root of values in a *.b grid. More...
program	regrd2
	Part of the NADCON5 NADCON5 Core Library , regrid data. More...
program	subtrc
	Part of the NADCON5 NADCON5 Core Library , Subtract one grid from another. More...
program	xyz2b
	Part of the NADCON5 NADCON5 Core Library , Converts GMT *.grd to a *.b NADCON style grid file. More...
subroutine	bicubic (z, glamn, glomn, dla, dlo, nla, nlo, maxla, maxlo, xla, xlo, val)
	Subroutine to perform a 2-D cubic ("bicubic") interpolation. More...
subroutine	bilin (data, glamn, glomn, dla, dlo, nla, nlo, maxla, maxlo, xla, xlo, val)
	Subroutine to perform bilinear interpolation. More...
subroutine	biquad (z, glamn, glomn, dla, dlo, nla, nlo, maxla, maxlo, xla, xlo, val)
	Subroutine to perform a 2-D quadratic ("biquadratic") interpolation. More...
subroutine	bwplotcv (ele, fname, bw, be, bs, bn, jm, b1, b2, maxplots, olddtm, newdtm, region, elecap, ij, igridsec, fn)
	Subroutine to make GMT calls to do a B/W coverage plot. More...
subroutine	bwplotvc (ele, fname, bw, be, bs, bn, jm, b1, b2, maxplots, olddtm, newdtm, region, elecap, ij, xvlon, xvlat, xllon, xllat, lorvog, lorvopc, igridsec, fn)
	Subroutine to make GMT calls to do a B/W vector plot. More...
subroutine	coplot (ele, fname, bw, be, bs, bn, jm, b1, b2, maxplots, olddtm, newdtm, region, elecap, ij, cptlo, cpthi, cptin6, suffixused, igridsec, fn)
	Subroutine to make GMT calls to do Color Raster Rendering of Gridded Data. More...
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. More...
subroutine	cpt (ave, std, csm, xlo, xhi, xin)
	This subroutine generates the color pallette variables for a GMT color plot. More...
subroutine	cpt2 (med, csm, xlo, xhi, xin)
	This subroutine generates the color pallette variables for a GMT color plot. More...
real function	cubterp (x, f0, f1, f2, f3)
	This function fits a cubic function through four points. More...
subroutine	getgridbounds (region, xn, xs, xw, xe)
	Subroutine to collect up the GRID boundaries for use in creating NADCON 5. More...
subroutine	getmag (x, ix)
	Subroutine to return the magnitude of a double precision value. More...
subroutine	getmapbounds (mapflag, maxplots, region, nplots, olddtm, newdtm, bw, be, bs, bn, jm, b1, b2, fn, lrv, rv0x, rv0y, rl0y)
	Subroutine to collect up the MAP boundaries for use in creating NADCON 5. More...
subroutine	gridstats (fname, ave, std, med)
	Subroutine to print grid statistics to stdout. More...
subroutine	indexxd (n, nd, arr, indx)
	Subroutine to perform ?? indexing on floating point data (double precision) More...
subroutine	indexxi (n, nd, arr, indx)
	Subroutine to perform ?? indexing on integer data. More...
integer *2 function	iselect2 (k, n, arr, nmax)
	Function to select an element of a partially filled, but packed multi dimensional array, integer*2 More...
real *4 function	onzd (x)
	Function to round a digit to one significant figure (one non zero digit), single precision. More...
real *8 function	onzd2 (x)
	Function to round a digit to one significant figure (one non zero digit), double precision. More...
subroutine	plotcoast (region, ifnum)
	Subroutine to write GMT-based commands to create a shoreline Write GMT-based commands to create a shoreline based on region. More...
real function	qterp (x, f0, f1, f2)
	This function fits a quadratic function through 3 points. More...
real *8 function	select2 (k, n, arr, nmax)
	Function to select an element of a partially filled, but packed multi dimensional array, double precision. More...
real *4 function	select2 (k, n, arr, nmax)
	Function to select an element of a partially filled, but packed multi dimensional array, single precision. More...
subroutine	vecstats (fname, n)
	Subroutine to tell us how many thinned vectors were used to make a grid. More...

Detailed Description

Programs and Functions which are responsible for computing the grid transformations used to build NADCON5.

The programs here can be considered "Library" Programs, or functions, which do the raw work of computing various transformations and conversions from *.b NADCON binary grid to a GMT style *.grd

Function Documentation

program addem

(

)

Part of the NADCON5 NADCON5 Core Library , adds one grid to another.

Program arguments

Arguments are newline terminated and read from standard input

When run from the command line, the program prints a prompt for each argument

They are enumerated here

Parameters

infileA	First Input File Name
infileB	Second Input File Name
outfile	Output File Name of A+B

Program Inputs:

• lin1 Input File A
• lin2 Input File B
• lout Output File to Write A+B

Definition at line 20 of file addem.f.

program b2xyz

(

)

Part of the NADCON5 NADCON5 Core Library , converts *.b grid to xyz

Program to convert standard "*.b" grid formatted data to a binary xyz (lon, lat, value) list, which can then be used by GMT for various things (like running the GMT routine "xyz2grd", to get a "*.grd" file, which is useful for plotting, etc)

Program arguments

Arguments are newline terminated and read from standard input

When run from the command line, the program prints a prompt for each argument

They are enumerated here

Parameters

infile

Input File Name

Program Inputs:

• Input File defined by infile

Program Outputs:

• temp.xyz

Definition at line 27 of file b2xyz.f.

subroutine bicubic	(	real*4, dimension(maxla,maxlo)	z,
			glamn,
			glomn,
			dla,
			dlo,
			nla,
			nlo,
			maxla,
			maxlo,
			xla,
			xlo,
			val
	)

Subroutine to perform a 2-D cubic ("bicubic") interpolation.

Performs interpolation at location "xla,xlo" off of grid "z", whose header information is the standard ".b" header information with additional inputs

Parameters

[in]	z	Input Grid
[in]	glamn	minimum latitude (real*8 decimal degrees)
[in]	glomn	minimum longitude (real*8 decimal degrees)
[in]	dla	latitude spacing (real*8 decimal degrees)
[in]	dlo	longitude spacing (real*8 decimal degrees)
[in]	nla	number of lat rows (integer*4)
[in]	nlo	number of lon cols (integer*4)
[in]	maxla	actual dimensioned size of "z" in rows
[in]	maxlo	actual dimensioned size of "z" in cols
[in]	xla	lat of pt for interpolation (real*8 dec. deg)
[in]	xlo	lon of pt for interpolation (real*8 dec. deg)
[out]	val	interpolated value (real*8)

Method:

Fit a 4x4 window over the random point. Unless the point is less than one grid spacing from the edge of the grid, it will fall in the inner 2x2 cell, and the 4x4 cell will be centered upon that.

Thus, if our point of interest is the asterisk, the window will look like this:

       .  .  .  .
       .  .  .  .
       .  .* .  . 
       .  .  .  .

Definition at line 37 of file bicubic.f.

Referenced by checkgrid().

c - Data is assumed real*4
      implicit real*8 (a-h,o-z)
      real*4 z(maxla,maxlo)

c - Internal use
      real*4 x,y,fx0,fx1,fx2,fx3,cubterp

c - Find which row should be LEAST when fitting
c - a 4x4 window around xla.  
      difla = (xla - glamn)
      ratla = difla / dla
      jla = int(ratla)
c - Fix any edge overlaps
      if(jla.lt.1)jla = 1
      if(jla.gt.nla-3)jla = nla-3

c - Find which column should be LEAST when fitting
c - a 4x4 window around xlo.
      diflo = (xlo - glomn)
      ratlo = diflo / dlo
      jlo = int(ratlo)
c - Fix any edge overlaps
      if(jlo.lt.1)jlo = 1
      if(jlo.gt.nlo-3)jlo = nlo-3

c - In the range of 0(westernmost) to
c - 3(easternmost) col, where is our
c - random lon value?  That is, x must
c - be real and fall between 0 and 3.
c - (Note, unless it is near an edge,
c - it will always be between 1 and 2)

      x=(xlo-dlo*(jlo-1)-glomn)/dlo

      if(x.lt.0.0)then
        write(6,100)x
        stop
      endif
  100 format(
     *'FATAL in bicubic:  x<0 : ',f20.10)

      if(x.gt.3.0)then
        write(6,101)x
        stop
      endif
  101 format(
     *'FATAL in bicubic:  x>3 : ',f20.10)

c - In the range of 0(southernmost) to
c - 3(northernmost) row, where is our
c - random lat value?  That is, x must
c - be real and fall between 0 and 3.
c - (Note, unless it is near an edge,
c - it will always be between 1 and 2)

      y=(xla-dla*(jla-1)-glamn)/dla

      if(y.lt.0.0)then
        write(6,102)y
        stop
      endif
  102 format(
     *'FATAL in bicubic:  y<0 : ',f20.10)

      if(y.gt.3.0)then
        write(6,103)y
        stop
      endif
  103 format(
     *'FATAL in bicubic:  y>3 : ',f20.10)

c - Now do the interpolation.  First, build a 1-D cubic function
c - east-west in the southermost row and interpolate to longitude
c - "xlo" (at "x" for 0<x<3).  Then do it in the 2nd row,
c - then again in the 3rd, and finally the 4th and most northern row.  
c - The last step is to fit a 1-D cubic function north-south at the
c - four previous interpolations, but this time to interpolate to 
c - latitude "xla" (at "y" for 0<y<3).  Obviously we 
c - could reverse the order, doing 4 north-south cubics
c - followed by one east-east and we'd get the same answer.


      fx0=cubterp(x,z(jla  ,jlo  ),z(jla  ,jlo+1),
     *              z(jla  ,jlo+2),z(jla  ,jlo+3))

      fx1=cubterp(x,z(jla+1,jlo  ),z(jla+1,jlo+1),
     *              z(jla+1,jlo+2),z(jla+1,jlo+3))

      fx2=cubterp(x,z(jla+2,jlo  ),z(jla+2,jlo+1),
     *              z(jla+2,jlo+2),z(jla+2,jlo+3))

      fx3=cubterp(x,z(jla+3,jlo  ),z(jla+3,jlo+1),
     *              z(jla+3,jlo+2),z(jla+3,jlo+3))

      val=dble(cubterp(y,fx0,fx1,fx2,fx3))

      
c     write(6,1001)glamn,xla,jla,y
c     write(6,1002)glomn,xlo,jlo,x
c     write(6,1003)z(jla+3,jlo  ),z(jla+3,jlo+1),
c    *             z(jla+3,jlo+2),z(jla+3,jlo+3)
c     write(6,1003)z(jla+2,jlo  ),z(jla+2,jlo+1),
c    *             z(jla+2,jlo+2),z(jla+2,jlo+3)
c     write(6,1003)z(jla+1,jlo  ),z(jla+1,jlo+1),
c    *             z(jla+1,jlo+2),z(jla+1,jlo+3)
c     write(6,1003)z(jla  ,jlo  ),z(jla  ,jlo+1),
c    *             z(jla  ,jlo+2),z(jla  ,jlo+3)
c     write(6,1004)fx0,fx1,fx2,fx3,val
c     pause

 1001 format(f14.10,1x,f14.10,1x,i8,1x,f8.5)
 1002 format(f14.10,1x,f14.10,1x,i8,1x,f8.5)
 1003 format(4(f15.8,1x))
 1004 format(
     *'fx0 = ',f15.8,/,
     *'fx1 = ',f15.8,/,
     *'fx2 = ',f15.8,/,
     *'fx3 = ',f15.8,/,
     *'val = ',f15.8)


      
      return

Here is the caller graph for this function:

subroutine bilin	(	real*4, dimension(maxla,maxlo)	data,
			glamn,
			glomn,
			dla,
			dlo,
			nla,
			nlo,
			maxla,
			maxlo,
			xla,
			xlo,
			val
	)

Subroutine to perform bilinear interpolation.

Performs a bilinear interpolation at location xla,xlo off of grid data, whose header information is the standard .b header information

Parameters

[in]	data	Input Data assumed to be real*4
[in]	glamn	minimum latitude (real*8 decimal degrees) *.b
[in]	glomn	minimum longitude (real*8 decimal degrees) *.b
[in]	dla	latitude spacing (real*8 decimal degrees) *.b
[in]	dlo	longitude spacing (real*8 decimal degrees) *.b
[in]	nla	number of lat rows (integer*4) *.b
[in]	nlo	number of lon cols (integer*4) *.b
[in]	maxla	actual dimensioned size of "data" in rows *.b
[in]	maxlo	actual dimensioned size of "data" in cols *.b
[in]	xla	lat of pt for interpolation (real*8 dec. def)
[in]	xlo	lon of pt for interpolation (real*8 dec. def)
[out]	val	interpolated value (real*8)

Definition at line 23 of file bilin.f.

Referenced by checkgrid().

c - Subroutine to perform a bilinear interpolation
c - at location "xla,xlo" off of grid "data", whose
c - header information is the standard ".b" header
c - information of
c -    glamn = minimum latitude   (real*8 decimal degrees)
c -    glomn = minimum longitude  (real*8 decimal degrees)
c -    dla   = latitude spacing   (real*8 decimal degrees)
c -    dlo   = longitude spacing  (real*8 decimal degrees)
c -    nla   = number of lat rows (integer*4)
c -    nlo   = number of lon cols (integer*4)
c -    maxla = actual dimensioned size of "data" in rows
c -    maxlo = actual dimensioned size of "data" in cols

c - data is assumed real*4
      implicit real*8 (a-h,o-z)
      real*4 data(maxla,maxlo)
      logical onedlat,onedlon

c - HACK 
c     do 8787 i=1,nla
c       do 8788 j=1,nlo
c         write(6,8789)i,j,data(i,j)
c8788   continue
c8787 continue
c8789 format(i5,1x,i5,1x,f20.10)

c - Find the row of xla
      difla = (xla - glamn) 
      ratla = difla / dla
      ila = int(ratla)+1
      gla0 = glamn + (ila-1)*dla

c - Find the col of xlo
      diflo = (xlo - glomn) 
      ratlo = diflo / dlo
      ilo = int(ratlo)+1
      glo0 = glomn + (ilo-1)*dlo


c - Do the interpolation.  
      t=(xlo-glo0)/dlo
      u=(xla-gla0)/dla
     
      val = (1.d0-t)*(1.d0-u)*data(ila  ,ilo  )  
     *    + (     t)*(1.d0-u)*data(ila  ,ilo+1)
     *    + (     t)*(     u)*data(ila+1,ilo+1)
     *    + (1.d0-t)*(     u)*data(ila+1,ilo  )

c - HACK
      return  
      write(6,198)xla,xlo
      write(6,200)ila,ilo,gla0,glo0
      write(6,199)xla-gla0,xlo-glo0
      write(6,201)ila,ilo,gla0,glo0,
     *data(ila,ilo)
      write(6,201)ila,ilo+1,gla0,glo0+dlo,
     *data(ila,ilo+1)
      write(6,201)ila+1,ilo,gla0+dla,glo0,
     *data(ila+1,ilo)
      write(6,201)ila+1,ilo+1,gla0+dla,glo0+dlo,
     *data(ila+1,ilo+1)
      write(6,202)'SW',(1.d0-t)*(1.d0-u)
      write(6,202)'SE',(     t)*(1.d0-u)
      write(6,202)'NW',(     t)*(     u)
      write(6,202)'NE',(1.d0-t)*(     u)
  202 format(
     *'Weight to ',a,' corner = ',f20.15)

  198 format(
     *6x,'bilin.f: lat/lon = ',f10.6,1x,f10.6)
  199 format(
     *6x,'bilin.f: diff from SW corner of cell = ',f10.6,1x,f10.6)
  200 format(
     *6x,'bilin.f: SW corner: ',i6,1x,i6,1x,f10.6,1x,f10.6)
  201 format(
     *6x,2(i6,1x),2(f10.6,1x),f20.15)
      

c     write(6,801)xla,xlo,ila,ilo,gla0,glo0,t,u,val
  801 format(f15.8,1x,f15.8,1x,i6,1x,i6,1x,
     *f15.8,1x,f15.8,1x,f8.6,1x,f8.6,1x,f20.10)

c     write(6,802)data(ila  ,ilo  ),
c    *data(ila  ,ilo+1),data(ila+1,ilo+1),
c    *data(ila+1,ilo  )
  802 format(4(f15.8,1x))
      
      return

Here is the caller graph for this function:

subroutine biquad	(	real*4, dimension(maxla,maxlo)	z,
			glamn,
			glomn,
			dla,
			dlo,
			nla,
			nlo,
			maxla,
			maxlo,
			xla,
			xlo,
			val
	)

Subroutine to perform a 2-D quadratic ("biquadratic") interpolation.

Performs a biquadratic interpolation at location xla,xlo off of grid z, whose header information is the standard ".b" header information

Parameters

[in]	z	Input Grid
[in]	glamn	minimum latitude (real*8 decimal degrees)
[in]	glomn	minimum longitude (real*8 decimal degrees)
[in]	dla	latitude spacing (real*8 decimal degrees)
[in]	dlo	longitude spacing (real*8 decimal degrees)
[in]	nla	number of lat rows (integer*4)
[in]	nlo	number of lon cols (integer*4)
[in]	maxla	actual dimensioned size of "z" in rows
[in]	maxlo	actual dimensioned size of "z" in cols
[in]	xla	lat of pt for interpolation (real*8 dec. deg)
[in]	xlo	lon of pt for interpolation (real*8 dec. deg)
[out]	val	interpolated value (real*8)

Method:

Fit a 3x3 window over the random point. The closest 2x2 points will surround the point. But based on which quadrant of that 2x2 cell in which the point falls, the 3x3 window could extend NW, NE, SW or SE from the 2x2 cell.

Definition at line 29 of file biquad.f.

Referenced by checkgrid().

c - Subroutine to perform a 2-D quadratic ("biquadratic")
c - interpolation at location "xla,xlo" off of grid "z", whose
c - header information is the standard ".b" header
c - information of
c -    glamn = minimum latitude   (real*8 decimal degrees)
c -    glomn = minimum longitude  (real*8 decimal degrees)
c -    dla   = latitude spacing   (real*8 decimal degrees)
c -    dlo   = longitude spacing  (real*8 decimal degrees)
c -    nla   = number of lat rows (integer*4)
c -    nlo   = number of lon cols (integer*4)
c -    maxla = actual dimensioned size of "z" in rows
c -    maxlo = actual dimensioned size of "z" in cols
c - Further input: 
c -    xla   = lat of pt for interpolation (real*8 dec. deg)
c -    xlo   = lon of pt for interpolation (real*8 dec. deg)
c - Output:
c -    val   = interpolated value (real*8)
c - 
c - Method:
c -   Fit a 3x3 window over the random point.  The closest
c -   2x2 points will surround the point.  But based on which
c -   quadrant of that 2x2 cell in which the point falls, the
c -   3x3 window could extend NW, NE, SW or SE from the 2x2 cell.

c - Data is assumed real*4
      implicit real*8 (a-h,o-z)
      real*4 z(maxla,maxlo)

c - Internal use
      real*4 x,y,fx0,fx1,fx2,qterp

c - Find which row should be LEAST when fitting
c - a 3x3 window around xla.
      difla = (xla - glamn)
      ratla = difla / (dla/2.d0)
      ila = int(ratla)+1
      if(mod(ila,2).ne.0)then
        jla = (ila+1)/2 - 1
      else
        jla = (ila  )/2
      endif
c - Fix any edge overlaps
      if(jla.lt.1)jla = 1
      if(jla.gt.nla-2)jla = nla-2

c - Find which column should be LEAST when fitting
c - a 3x3 window around xlo.
      diflo = (xlo - glomn)
      ratlo = diflo / (dlo/2.d0)
      ilo = int(ratlo)+1
      if(mod(ilo,2).ne.0)then
        jlo = (ilo+1)/2 - 1
      else
        jlo = (ilo  )/2
      endif
c - Fix any edge overlaps
      if(jlo.lt.1)jlo = 1
      if(jlo.gt.nlo-2)jlo = nlo-2

c - In the range of 0(westernmost) to
c - 2(easternmost) col, where is our
c - random lon value?  That is, x must
c - be real and fall between 0 and 2.

      x=(xlo-dlo*(jlo-1)-glomn)/dlo

      if(x.lt.0.0)then
        val = 1.d30
        write(6,100)x,val
c       stop
        return
      endif
  100 format(
     *'FATAL in biquad:  x<0 : ',f20.10,/,
     *' --> Returning with val = ',f40.1)

      if(x.gt.2.0)then
        val = 1.d30
        write(6,101)x,val
c       stop
        return
      endif
  101 format(
     *'FATAL in biquad:  x>2 : ',f20.10,/,
     *' --> Returning with val = ',f40.1)

c - In the range of 0(southernmost) to
c - 2(northernmost) row, where is our
c - random lat value?  That is, x must
c - be real and fall between 0 and 2.

      y=(xla-dla*(jla-1)-glamn)/dla

      if(y.lt.0.0)then
        val = 1.d30
        write(6,102)y,val
c       stop
        return
      endif
  102 format(
     *'FATAL in biquad:  y<0 : ',f20.10,/,
     *' --> Returning with val = ',f40.1)

      if(y.gt.2.0)then
        val = 1.d30
        write(6,103)y,val
c       stop
        return
      endif
  103 format(
     *'FATAL in biquad:  y>2 : ',f20.10,/,
     *' --> Returning with val = ',f40.1)

c - Now do the interpolation.  First, build a paraboloa
c - east-west the southermost row and interpolate to longitude
c - "xlo" (at "x" for 0<x<2).  Then do it in the middle
c - row, then finally the northern row.  The last step
c - is to fit a parabola north-south at the three previous
c - interpolations, but this time to interpolate to 
c - latitude "xla" (at "y" for 0<y<2).  Obviously we 
c - could reverse the order, doing 3 north-south parabolas
c - followed by one east-east and we'd get the same answer.

      fx0=qterp(x,z(jla  ,jlo),z(jla  ,jlo+1),z(jla  ,jlo+2))
      fx1=qterp(x,z(jla+1,jlo),z(jla+1,jlo+1),z(jla+1,jlo+2))
      fx2=qterp(x,z(jla+2,jlo),z(jla+2,jlo+1),z(jla+2,jlo+2))
      val=dble(qterp(y,fx0,fx1,fx2))

c     goto 888

c     write(6,*) ' '
c     write(6,*) ' ---------------------------------'
c     write(6,*) ' Inside biquad.f'
c     write(6,*) ' '

c     write(6,*) ' Southern Row of 3x3 window: '
c     write(6,*) ' jla  ,jlo  ,z(jla  ,jlo  ) = ',
c    *jla  ,jlo  ,z(jla  ,jlo  )
c     write(6,*) ' jla  ,jlo+1,z(jla  ,jlo+1) = ',
c    *jla  ,jlo+1,z(jla  ,jlo+1)
c     write(6,*) ' jla  ,jlo+2,z(jla  ,jlo+2) = ',
c    *jla  ,jlo+2,z(jla  ,jlo+2)
c     write(6,*) ' fx0 = ',fx0

c     write(6,*) ' Middle   Row of 3x3 window: '
c     write(6,*) ' jla+1,jlo  ,z(jla+1,jlo  ) = ',
c    *jla+1,jlo  ,z(jla+1,jlo  )
c     write(6,*) ' jla+1,jlo+1,z(jla+1,jlo+1) = ',
c    *jla+1,jlo+1,z(jla+1,jlo+1)
c     write(6,*) ' jla+1,jlo+2,z(jla+1,jlo+2) = ',
c    *jla+1,jlo+2,z(jla+1,jlo+2)
c     write(6,*) ' fx1 = ',fx1

c     write(6,*) ' Northern Row of 3x3 window: '
c     write(6,*) ' jla+2,jlo  ,z(jla+2,jlo  ) = ',
c    *jla+2,jlo  ,z(jla+2,jlo  )
c     write(6,*) ' jla+2,jlo+1,z(jla+2,jlo+1) = ',
c    *jla+2,jlo+1,z(jla+2,jlo+1)
c     write(6,*) ' jla+2,jlo+2,z(jla+2,jlo+2) = ',
c    *jla+2,jlo+2,z(jla+2,jlo+2)
c     write(6,*) ' fx2 = ',fx2

c     write(6,*) ' Final interpolated value = ',val

c 888 continue


      return

Here is the caller graph for this function:

subroutine bwplotcv	(	character*3	ele,
		character*200	fname,
		real*8, dimension(maxplots)	bw,
		real*8, dimension(maxplots)	be,
		real*8, dimension(maxplots)	bs,
		real*8, dimension(maxplots)	bn,
		real*4, dimension(maxplots)	jm,
		real*4, dimension(maxplots)	b1,
		real*4, dimension(maxplots)	b2,
		integer*4	maxplots,
		character*10	olddtm,
		character*10	newdtm,
		character*10	region,
		character*3	elecap,
			ij,
			igridsec,
		character*10, dimension(maxplots)	fn
	)

Subroutine to make GMT calls to do a B/W coverage plot.

Changelog

2016 08 29:

Updated the -R and -B initial calls to 6 decimal places

2016 08 25:

.gmtdefaults4 has been changed so X_ORIGIN is equal to 0.0 Center the plot with "-Xc" at first "psxy" call Remove all "-JM**i+" references, and just use the actual width (jm) that came out of the "getmapbounds" routine and was sent here.

2016 07 21:

Modified use of JM command based on new forced sizes.

2015 02 15:

Updated to allow this subroutine to work earlier (in makeplotfiles01()), before igridsec was defined. See DRU-11, p. 139

Definition at line 26 of file bwplotcv.f.

References plotcoast().

Referenced by makeplotfiles01(), makeplotfiles02(), and makeplotfiles03().

c - 2016 08 29:
c   Updated the -R and -B initial 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 first "psxy" call
c     Remove all "-JM**i+" references, and just use the actual
c       width (jm) that came out of the "getmapbounds" routine and
c       was sent here.

c - 2016 07 21:
c - Modified use of JM command based on new forced sizes.  

c - Subroutine to make GMT calls to do a B/W coverage plot

c - Updated 10/2/2015 to allow this subroutine to work
c - earlier (in makeplotfiles01.f), before "igridsec" was defined.
c - See DRU-11, p. 139

      integer*4 maxplots
      character*3 ele,elecap
      character*200 fname
      character*10 olddtm,newdtm,region
      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*5 extra
      character*20 gridnote

c ----------------------------------------------------
c - All coverage type files begin with "cv". 
c - The 3rd element tells me all, thinned, dropped or RMS 
c - The 4th/5th elements tell me if these are coordinate
c - differneces or double differences.  Elements 6,7,8
c - tell me what data we're plotting.
c ----------------------------------------------------


c ----------------------------------------------------
c - FAILSAFES: BEGIN
c ----------------------------------------------------
      if(fname(3:3).eq.'t')then
        extra='-thin'
      elseif(fname(3:3).eq.'d')then
        extra='-drop'
      elseif(fname(3:3).eq.'a')then
        extra='-all '
      elseif(fname(3:3).eq.'r')then
        extra='-RMSd'
      else
        write(6,1)trim(fname)
        stop
      endif
    1 format('FATAL in bwplotcv.  Bad character in spot 3: ',a)

      if(fname(1:2).ne.'cv')then
        write(6,2)trim(fname)
        stop
      endif
    2 format('FATAL in bwplotcv.  Bad character in spots 1-2: ',a)
      
      if(.not.(fname(4:5).eq.'cd' .or. fname(4:5).eq.'dd'))then
        write(6,3)trim(fname)
        stop
      endif
    3 format('FATAL in bwplotcv.  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 bwplotcv.  Bad match of fname / ele: ',a,1x,a)

c - Just in case we forgot to set "igridsec" to be -1
c - when we came here from makeplotfiles01:
      if(fname(1:3).eq.'cvacd')then
        igridsec = -1
      endif

      if(igridsec.le.0)then
        gridnote = ''
      else
        write(gridnote,10)igridsec
      endif
   10 format('(',i0,' sec)')



c ----------------------------------------------------
c - FAILSAFES: END
c ----------------------------------------------------



c ----------------------------------------------
c - GMT COMMANDS: BEGIN
c ----------------------------------------------

c - Header of commands/echoes:
      write(99,991)ele,extra,trim(region),trim(fn(ij)),
     *ele,extra,trim(region),trim(fn(ij))

c - Write out the actual dots, and the title, with PSXY command
      write(99,904)trim(fname),bw(ij),be(ij),bs(ij),bn(ij),
     *jm(ij),b1(ij),b2(ij),trim(newdtm),trim(olddtm),elecap,
     *extra,trim(gridnote),trim(region),trim(fn(ij))

c - ALWAYS plot the coast last, as it closes the PS file
      call plotcoast(region,99)

c - Convert PS to JPG
      write(99,905)

c - Rename the JPG to my naming scheme 
      write(99,910)trim(fname),trim(fn(ij))
c ----------------------------------------------
c - GMT COMMANDS: END
c ----------------------------------------------



  991 format(
     *'# -----------------------------------------------------',/,
     *'# coverage in ',a,a,1x,a,1x,a,/,
     *'# -----------------------------------------------------',/,
     *'echo ...coverage in ',a,a,1x,a,1x,a)


c - Plot the actual coverage points
c - 2016 08 29
  904 format('psxy ',a,' -Xc -R',f0.6,'/',f0.6,'/',sp,f0.6,'/',f0.6,
     *  ss,' -JM',f3.1,'i -B',f0.6,'/',f0.6,':."',
     *  'NADCON v5.0 ',a,' minus ',a,' ',a3,a5,a,
     *  1x,a,'-',a,
     *  '": -Sc0.02i ',
     *  '-Gblack -K > plot.ps')



c - 2016 08 25:
c 904 format('psxy ',a,' -Xc -R',f0.2,'/',f0.2,'/',sp,f0.2,'/',f0.2,
c    *  ss,' -JM',f3.1,'i -B',f0.2,'/',f0.2,':."',
c    *  'NADCON v5.0 ',a,' minus ',a,' ',a3,a5,a,
c    *  1x,a,'-',a,
c    *  '": -Sc0.02i ',
c    *  '-Gblack -K > plot.ps')
c - 2016 07 21
c 904 format('psxy ',a,' -R',f0.2,'/',f0.2,'/',sp,f0.2,'/',f0.2,
c    *  ss,' -JM',f3.1,'i -B',f0.2,'/',f0.2,':."',
c    *  'NADCON v5.0 ',a,' minus ',a,' ',a3,a5,a,
c    *  1x,a,'-',a,
c    *  '": -Sc0.02i ',
c    *  '-Gblack -K > plot.ps')
c 904 format('psxy ',a,' -R',f0.2,'/',f0.2,'/',sp,f0.2,'/',f0.2,
c    *  ss,' -JM',f3.1,'i+ -B',f0.2,'/',f0.2,':."',
c    *  'NADCON v5.0 ',a,' minus ',a,' ',a3,a5,a,
c    *  1x,a,'-',a,
c    *  '": -Sc0.02i ',
c    *  '-Gblack -K > plot.ps')

c - 905 = Convert PS to JPG
  905 format('ps2raster plot.ps -Tj -P -A ')


c - Renaming 
  910 format('mv -f plot.jpg ',a,'.',a,'.jpg',/,
     *'rm -f plot.ps')
c 910 format('mv -f plot.jpg cv',a1,a3,'.',a,'.',a,'.',a,
c    *'.',i0,'.',a,'.jpg',/,'rm -f plot.ps')
 1910 format('mv -f plot.jpg cv',a3,'.',a,'.',a,'.',a,
     *'.',a,'.jpg',/,'rm -f plot.ps')

      return

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine bwplotvc	(	character*3	ele,
		character*200	fname,
		real*8, dimension(maxplots)	bw,
		real*8, dimension(maxplots)	be,
		real*8, dimension(maxplots)	bs,
		real*8, dimension(maxplots)	bn,
		real*4, dimension(maxplots)	jm,
		real*4, dimension(maxplots)	b1,
		real*4, dimension(maxplots)	b2,
		integer*4	maxplots,
		character*10	olddtm,
		character*10	newdtm,
		character*10	region,
		character*3	elecap,
			ij,
			xvlon,
			xvlat,
			xllon,
			xllat,
		real*8	lorvog,
		real*8	lorvopc,
			igridsec,
		character*10, dimension(maxplots)	fn
	)

Subroutine to make GMT calls to do a B/W vector plot.

Changelog

2016 08 29:

Expanded the refence vector calls to be 6 decimal places, as well as the initial -R call for S/N/W/E and also the -B part of that call

2016 08 25:

.gmtdefaults4 has been changed so X_ORIGIN is equal to 0.0 Center the plot with "-Xc" at first "psxy" call Remove all "-JM**i+" references, and just use the actual width (jm) that came out of the "getmapbounds" routine and was sent here.

2016 07 29:

Updated the reference vector call to have the "-N" option, so it'll plot outside the map

2016 07 21:

Modified use of JM command based on new forced sizes.

2015 02 15:

Updated to allow this subroutine to work earlier (in makeplotfiles01()), before igridsec was defined. See DRU-11, p. 139

Definition at line 33 of file bwplotvc.f.

References plotcoast().

Referenced by checkgrid(), makeplotfiles01(), makeplotfiles02(), and makeplotfiles03().

c - 2016 08 29:
c     Expanded the refence vector calls to be 6 decimal places, as well
c     as the initial -R call for S/N/W/E
c     and also the -B part of that call

c - 2016 08 25:
c     .gmtdefaults4 has been changed so X_ORIGIN is equal to 0.0
c     Center the plot with "-Xc" at first "psxy" call
c     Remove all "-JM**i+" references, and just use the actual
c       width (jm) that came out of the "getmapbounds" routine and
c       was sent here.


c - 2016 07 29:  Updated the reference vector call to 
c - have the "-N" option, so it'll plot outside the map

c - 2016 07 21:  Updated JM usage for new forced sizes

c     subroutine bwplotvc(ele,fname,bw,be,bs,bn,jm,b1,b2,maxplots,
c    *olddtm,newdtm,region,elecap,ij,xvlon,xvlat,xllon,xllat,ncm,
c    *q1,igridsec,fn)
c - Subroutine to make GMT calls to do a B/W vector plot

c - Updated 10/2/2015 to allow this subroutine to work
c - earlier (in makeplotfiles01.f), before "igridsec" was defined.
c - See DRU-11, p. 139

      implicit real*8(a-h,o-z)
      integer*4 maxplots
      character*3 ele,elecap
      character*200 fname
      character*10 olddtm,newdtm,region
      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*5 extra
      character*20 units
      character*20 gridnote
      character*20 gridnote2

      real*8 lorvopc,lorvog

c ----------------------------------------------------
c - All vector  type files begin with either "vm"
c - (meters) or "vs" (arcseconds).
c - The 3rd element tells me all, thinned, dropped or RMS
c - The 4th/5th elements tell me if these are coordinate
c - differnces or double differences.  Elements 6,7,8
c - tell me what data we're plotting.
c ----------------------------------------------------


c ----------------------------------------------------
c - FAILSAFES: BEGIN
c ----------------------------------------------------
      if(fname(3:3).eq.'t')then
        extra='-thin'
      elseif(fname(3:3).eq.'d')then
        extra='-drop'
      elseif(fname(3:3).eq.'a')then
        extra='-all '
      elseif(fname(3:3).eq.'r')then
        extra='-RMSd'
      else
        write(6,1)trim(fname)
        stop
      endif
    1 format('FATAL in bwplotvc.  Bad character in spot 3: ',a)

      if(.not.(fname(1:2).eq.'vm'.or.fname(1:2).eq.'vs'))then
        write(6,2)trim(fname)
        stop
      endif
    2 format('FATAL in bwplotvc.  Bad character in spots 1-2: ',a)

      if(.not.(fname(4:5).eq.'cd' .or. fname(4:5).eq.'dd' .or.
     *         fname(4:5).eq.'gi'))then
        write(6,3)trim(fname)
        stop
      endif
    3 format('FATAL in bwplotvc.  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 bwplotvc.  Bad match of fname / ele: ',a,1x,a)

c - Just in case we forgot to set "igridsec" to be -1
c - when we came here from makeplotfiles01:
      if(fname(1:3).eq.'vmacd' .or. fname(1:3).eq.'vsacd')then
        igridsec = -1
      endif

      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 bwplotvc.  Bad units in name: ',a)

      if(igridsec.le.0)then
        gridnote = ''
        gridnote2 = ''
      else
        write(gridnote,10)igridsec
        write(gridnote2,11)igridsec
      endif
   10 format('(',i0,' sec)')
   11 format(i0,'.')

c ----------------------------------------------------
c - FAILSAFES: END
c ----------------------------------------------------


c ----------------------------------------------
c - GMT COMMANDS: BEGIN
c ----------------------------------------------

c - Header of commands/echoes:
      write(99,991)ele,extra,trim(units),trim(region),trim(fn(ij)),
     *ele,extra,trim(units),trim(region),trim(fn(ij))

c - Plot the actual vectors, and the title, with PSXY command
      write(99,903)trim(fname),bw(ij),be(ij),bs(ij),bn(ij),
     *jm(ij),b1(ij),b2(ij),trim(newdtm),trim(olddtm),elecap,
     *extra,trim(gridnote),trim(region),trim(fn(ij)),fname(2:5)



c - Plot the reference vector
      write(99,908)xvlon,xvlat,90.d0,lorvopc
c     write(99,908)xvlon,xvlat,90.d0,q1

c - Label the reference vector
c - HACK
c     write(6,807)trim(fname),trim(units)
c 807 format('File: ',a,/,'Units: ',a)
      if    (trim(units).eq.'arcseconds')then
        write(99,1909)xllon,xllat,lorvog,trim(units)
      elseif(trim(units).eq.'meters')then
        write(99, 909)xllon,xllat,lorvog,trim(units)
      endif

c - ALWAYS plot the coast last, as it closes the PS file
      call plotcoast(region,99)

c - Convert PS to JPG
      write(99,905)

c - Rename the JPG to my naming scheme
c     write(99,906)trim(fname),trim(gridnote2),trim(fn(ij))
      write(99,906)trim(fname),trim(fn(ij))

c ----------------------------------------------
c - GMT COMMANDS: END
c ----------------------------------------------

  991 format(
     *'# -----------------------------------------------------',/,
     *'# vectors in ',a,a,1x,a,1x,a,1x,a,/,
     *'# -----------------------------------------------------',/,
     *'echo ...vectors in ',a,a,1x,a,1x,a,1x,a)

c - Plot the actual vectors
c - 2016 08 29
  903 format('psxy ',a,' -Xc -R',f0.6,'/',f0.6,'/',sp,f0.6,'/',f0.6,
     *  ss,' -JM',f3.1,'i -B',f0.6,'/',f0.6,':."',
     *  'NADCON v5.0 ',a,' minus ',a,' ',a3,a5,a,
     *  1x,a,'-',a,1x,a,
     *'": -SV0.0001i/0.02i/0.02i ',
     *  '-Gblack -K > plot.ps')
c - 2016 08 25:
c 903 format('psxy ',a,' -Xc -R',f0.2,'/',f0.2,'/',sp,f0.2,'/',f0.2,
c    *  ss,' -JM',f3.1,'i -B',f0.2,'/',f0.2,':."',
c    *  'NADCON v5.0 ',a,' minus ',a,' ',a3,a5,a,
c    *  1x,a,'-',a,1x,a,
c    *'": -SV0.0001i/0.02i/0.02i ',
c    *  '-Gblack -K > plot.ps')

c 903 format('psxy ',a,' -R',f0.2,'/',f0.2,'/',sp,f0.2,'/',f0.2,
c    *  ss,' -JM',f3.1,'i -B',f0.2,'/',f0.2,':."',
c    *  'NADCON v5.0 ',a,' minus ',a,' ',a3,a5,a,
c    *  1x,a,'-',a,1x,a,
c    *'": -SV0.0001i/0.02i/0.02i ',
c    *  '-Gblack -K > plot.ps')
c 903 format('psxy ',a,' -R',f0.2,'/',f0.2,'/',sp,f0.2,'/',f0.2,
c    *  ss,' -JM',f3.1,'i+ -B',f0.2,'/',f0.2,':."',
c    *  'NADCON v5.0 ',a,' minus ',a,' ',a3,a5,a,
c    *  1x,a,'-',a,1x,a,
c    *'": -SV0.0001i/0.02i/0.02i ',
c    *  '-Gblack -K > plot.ps')

c - 908 = Plot a REFERENCE vector
c - 2016 08 29:  More decimal places
  908 format('psxy -SV0.0001i/0.02i/0.02i -N -R -O -K -JM',
     *  ' -Gred >> plot.ps << !',/,
     *  f10.6,1x,f10.6,1x,f5.1,1x,f9.1,/,
     *  '!')
c - 2016 07 29:  Make reference vector appear outside of map
c 908 format('psxy -SV0.0001i/0.02i/0.02i -N -R -O -K -JM',
c    *  ' -Gred >> plot.ps << !',/,
c    *  f6.2,1x,f6.2,1x,f5.1,1x,f9.1,/,
c    *  '!')
c 908 format('psxy -SV0.0001i/0.02i/0.02i -R -O -K -JM',
c    *  ' -Gred >> plot.ps << !',/,
c    *  f6.2,1x,f6.2,1x,f5.1,1x,f9.1,/,
c    *  '!')

c - 909 = Label the REFERENCE vector
c - 2016 08 29:  More decimal places
  909 format('pstext -N -O -K -R -JM -Gred >> plot.ps << !',/,
     *  f10.6,1x,f10.6,1x,'12 0 1 TL ',f10.3,1x,a,/,
     *  '!')
c - 2016 07 29:  Use -N to plot it outside
c 909 format('pstext -N -O -K -R -JM -Gred >> plot.ps << !',/,
c    *  f6.2,1x,f6.2,1x,'12 0 1 TL ',f10.3,1x,a,/,
c    *  '!')
c - 2016 08 29:  More decimal places
 1909 format('pstext -N -O -K -R -JM -Gred >> plot.ps << !',/,
     *  f10.6,1x,f10.6,1x,'12 0 1 TL ',f10.6,1x,a,/,
     *  '!')
c1909 format('pstext -N -O -K -R -JM -Gred >> plot.ps << !',/,
c    *  f6.2,1x,f6.2,1x,'12 0 1 TL ',f10.6,1x,a,/,
c    *  '!')
c 909 format('pstext -O -K -R -JM -Gred >> plot.ps << !',/,
c    *  f6.2,1x,f6.2,1x,'12 0 1 TL ',f10.3,1x,a,/,
c    *  '!')
c1909 format('pstext -O -K -R -JM -Gred >> plot.ps << !',/,
c    *  f6.2,1x,f6.2,1x,'12 0 1 TL ',f10.6,1x,a,/,
c    *  '!')

c 909 format('pstext -O -K -R -JM -Gred >> plot.ps << !',/,
c    *  f6.2,1x,f6.2,1x,'12 0 1 TL ',i6,' cm',/,
c    *  '!')

c - 905 = Convert PS to JPG
  905 format('ps2raster plot.ps -Tj -P -A ')

c - Rename to our naming scheme
c 906 format('mv -f plot.jpg ',a,
c    *'.',a,a,'.jpg',/,'rm -f plot.ps')
  906 format('mv -f plot.jpg ',a,
     *'.',a,'.jpg',/,'rm -f plot.ps')

      return

Here is the call graph for this function:

Here is the caller graph for this function:

program convlv

(

)

Part of the NADCON5 NADCON5 Core Library , Convolves two grids.

Convolves one grid against another

\[ c(i,j) = a(i,j) * b(i,j) \]

Program arguments

Arguments are newline terminated and read from standard input

When run from the command line, the program prints a prompt for each argument

They are enumerated here

Parameters

infileA	First Input File Name
infileB	Second Input File Name
outfile	Output File Name of A*B

Program Inputs:

• lin1 Input File A
• lin2 Input File B
• lout Output File to Write A+B

Definition at line 25 of file convlv.f.

subroutine coplot	(	character*3	ele,
		character*200	fname,
		real*8, dimension(maxplots)	bw,
		real*8, dimension(maxplots)	be,
		real*8, dimension(maxplots)	bs,
		real*8, dimension(maxplots)	bn,
		real*4, dimension(maxplots)	jm,
		real*4, dimension(maxplots)	b1,
		real*4, dimension(maxplots)	b2,
		integer*4	maxplots,
		character*10	olddtm,
		character*10	newdtm,
		character*10	region,
		character*3	elecap,
			ij,
			cptlo,
			cpthi,
			cptin6,
		character*200	suffixused,
			igridsec,
		character*10, dimension(maxplots)	fn
	)

Subroutine to make GMT calls to do Color Raster Rendering of Gridded Data.

Changelog

2016 09 08:

Had to up the D_FORMAT default to %.3G because tight scalebar ranges with the newly allowed "more free" average values were showing repeating values when only 2 digits could be shown.

2016 09 07:

Had to add lines pre/post "makecpt" to change the D_FORMAT. This is because I had been forcing the "scave" in cpt.f to be ONZD. But that yielded bad values sometimes, so I switched it. With that switch, the scave could have lots of digits. Well, that means the newly adopted "D_FORMAT" of %.2G was insufficient for the CPT table. Who knew that D_FORMAT affected that! Anyway, so change "D_FORMAT" pre/post all makecpt calls.

2016 08 30:

See item #39 in Google ToDo list Changed "grdcontour" to have a blank "-R" call so it'll mimic whatever decimal places are in the "grdimage" call that came before it.

2016 08 29:

Updated the initial -R and -B calls to 6 decimal places

2016 08 25:

-.gmtdefaults4 has been changed so X_ORIGIN is equal to 0.0

• Center the plot with "-Xc" at grdimage
• Center the scalebar by setting its Xcoordinate, which runs in "plot frame" coordinates (0/0 at lower left) , to be equal to "jm/2"
• Force the scale bar to be exactly 4 inches wide, always
• Change the format for "makecpt" from 0.6 to 0.10

2016 07 29:

Update to put more data into comment/echo

Also,

• forced the -Ctemp.cpt option in "grdcontour" to make the contours line up with the color palette
• For d3 plots, to drop all contours
• For d3 plots to only use "coverage" part only if "ij" is not "1"
• Same for "09" and "ete" grids

2016 07 21:

• Set the "JM" code in "grdcontour" to just be "-JM" and let it therefore run with whatever JM size is used in "grdimage"
• Set the "A" code in "grdcontour" to be "-A-" which should turn off the labels on all contours
• Fixed size of scale bar

2016 03 01:

1. Changed to a continuous color plot
   • Get rid of "-Z" in "makecpt"
2. Changed to an 8 color, from 6 color, plot (without changing the RANGE yet...see DRU-12, p. 19)
   • Change varible that comes in from "cptin" to "cptin6"
   • Compute cptin = cptin6 * 0.75d0 immediately

Update 2016 02 29:

1. Removed all shading from color plots
   • Get rid of "grdgradient" call
   • Remove from "grdimage" the "-Itempi.grd" part
   • Remove the "rm -f tempi.grd" line

Definition at line 70 of file coplot.f.

References plotcoast().

Referenced by makeplotfiles02(), and makeplotfiles03().

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
      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*20 gridnote
      character*20 gridnote2

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.'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
        write(gridnote,10)igridsec
        write(gridnote2,11)igridsec
      endif
    6 format('FATAL in coplot.  Bad igridsec: ',i0)
   10 format('(',i0,' sec)')
   11 format(i0,'.')

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 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     write(99,991)ele,trim(units),trim(region),trim(fn(ij)),
c    *ele,trim(units),trim(region),trim(fn(ij))

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 - Begin 2016 03 01:
c 901 format(
c    *'makecpt -Crainbow -T',f0.6,'/',f0.6,'/',f0.6,
c    *' -Z > temp.cpt')
c 901 format(
c    *'makecpt -Crainbow -T',f0.6,'/',f0.6,'/',f0.6,
c    *' > temp.cpt')
c - 2016 08 25
c 901 format(
c    *'makecpt -Crainbow -T',f0.10,'/',f0.10,'/',f0.10,
c    *' > temp.cpt')
c - 2016 09 07
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 %.2G')
c - 2016 09 08
  901 format(
     *'gmtset D_FORMAT %.12G',/,
     *'makecpt -Crainbow -T',f0.10,'/',f0.10,'/',f0.10,
     *' > temp.cpt',/,
     *'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 - Updated 2016 02 29:
c     write(99,902)fname(1:8),trim(suffixused)
  902 format(
     *'grdgradient ',a,'.',a,
     *'.grd -N4.8 -A90 -M -Gtempi.grd')
c - Until here 2016 02 29

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,igridsec,
     *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:
  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,
     *':."NADCON v5.0 ',a,' minus ',a,' ',a3,
     *'(',i0,' sec)',
     *1x,a,'-',a,1x,a,
     *':" -Ctemp.cpt -K > plot.ps')

c - GMT call to put a color scale on the map.  Add it
c - to file "plot.ps".

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 -B/:',a,
     *': -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 - GMT call to put contours on the map.
c - 2016 07 29
c - 2016 07 21

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 - 2016 07 28 - Skip to here (skipping contours) for d3, 09 or "ete"  plots
  807 continue

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 - GMT call to convert the plot.ps file to plot.jpg.
      write(99,905)
  905 format(
     *'ps2raster plot.ps -Tj -P -A')

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

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine coplotwcv	(	character*3	ele,
		character*200	fname,
		real*8, dimension(maxplots)	bw,
		real*8, dimension(maxplots)	be,
		real*8, dimension(maxplots)	bs,
		real*8, dimension(maxplots)	bn,
		real*4, dimension(maxplots)	jm,
		real*4, dimension(maxplots)	b1,
		real*4, dimension(maxplots)	b2,
		integer*4	maxplots,
		character*10	olddtm,
		character*10	newdtm,
		character*10	region,
		character*3	elecap,
			ij,
			cptlo,
			cpthi,
			cptin6,
		character*200	suffixused,
			igridsec,
		character*10, dimension(maxplots)	fn,
		character*200	cvfname
	)

Subroutine to make GMT calls to do a color raster rendering of gridded data, with coverage overlaid.

Changelog

2016 09 08:

Had to up the D_FORMAT default to %.3G because tight scalebar ranges with the newly allowed "more free" average values were showing repeating values when only 2 digits could be shown.

2016 09 07:

Had to add lines pre/post "makecpt" to change the D_FORMAT. This is because I had been forcing the "scave" in cpt.f to be ONZD. But that yielded bad values sometimes, so I switched it. With that switch, the scave could have lots of digits. Well, that means the newly adopted "D_FORMAT" of %.2G was insufficient for the CPT table. Who knew that D_FORMAT affected that! Anyway, so change "D_FORMAT" pre/post all makecpt calls.

2016 08 30:

See item #39 in Google ToDo list Changed "grdcontour" to have a blank "-R" call so it'll mimic whatever decimal places are in the "grdimage" call that came before it.

2016 08 29:

Updated the initial -R and -B calls to 6 decimal places

2016 08 25:

• .gmtdefaults4 has been changed so X_ORIGIN is equal to 0.0
• Center the plot with "-Xc" at grdimage
• Center the scalebar by setting its Xcoordinate, which runs in "plot frame" coordinates (0/0 at lower left) , to be equal to "jm/2"
• Force the scale bar to be exactly 4 inches wide, always
• Change the format for "makecpt" from 0.6 to 0.10

2016 07 29:

• Update to put more data into comment/echo
• Forced grdcontour with -Ctemp.cpt to align contours with color palette

2016 07 28:

• For d3 plots, to drop all contours
• For d3 plots to only use "coverage" part only if "ij" is not "1"
• Same for "09" and "ete" grids

2016 07 21:

• Set the "JM" code in "grdcontour" to just be "-JM" and let it therefore run with whatever JM size is used in "grdimage"
• Set the "A" code in "grdcontour" to be "-A-" which should turn off the labels on all contours
• Fixed the size of the scale bar

2016 03 01:

1. Changed to a continuous color plot
   • Get rid of "-Z" in "makecpt"
2. Changed to an 8 color, from 6 color, plot (without changing the RANGE yet...see DRU-12, p. 19)
   • Change varible that comes in from "cptin" to "cptin6"
   • Compute cptin = cptin6 * 0.75d0 immediately

2016 02 29:

1. Removed all shading from color plots
   • Get rid of "grdgradient" call
   • Remove from "grdimage" the "-Itempi.grd" part
   • Remove the "rm -f tempi.grd" line

Definition at line 68 of file coplotwcv.f.

References plotcoast().

Referenced by makeplotfiles02().

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 -   Forced grdcontour with -Ctemp.cpt to align contours with color palette

c - Update 2016 07 28 to:
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 coplotwcv(ele,fname,bw,be,bs,bn,jm,b1,b2,maxplots,
c    *olddtm,newdtm,region,elecap,ij,cptlo,cpthi,cptin,suffixused,
c    *igridsec,fn,cvfname)

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 the size of the 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  

      implicit real*8(a-h,o-z)

c - Subroutine to make GMT calls to do a color raster rendering
c - of gridded data, with coverage overlaid.

      integer*4 maxplots
      character*3 ele,elecap
      character*200 fname,cvfname
      character*10 olddtm,newdtm,region
      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*20 gridnote
      character*20 gridnote2

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'))then
        write(6,2)trim(fname)
        stop
      endif
    2 format('FATAL in coplotwcv.  Bad character in spots 1-3: ',a)

      if(.not.(fname(4:5).eq.'cd' .or. fname(4:5).eq.'dd'))then
        write(6,3)trim(fname)
        stop
      endif
    3 format('FATAL in coplotwcv.  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 coplotwcv.  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 coplotwcv.  Bad units in name: ',a)

      if(igridsec.le.0)then
        write(6,6)igridsec
      else
        write(gridnote,10)igridsec
        write(gridnote2,11)igridsec
      endif
    6 format('FATAL in coplotwcv.  Bad igridsec: ',i0)
   10 format('(',i0,' sec)')
   11 format(i0,'.')

c ----------------------------------------------------
c - FAILSAFES: END
c ----------------------------------------------------

c ------------------------------------------------------------------
c - Contour interval and labeling
c ------------------------------------------------------------------
c - 2016 07 28:  Scrapped the whole thing about scaling cptin.
c - Just bring in the interval for an 8 color spread, period.
      cptin = cptin6

c - 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 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     write(99,991)ele,trim(units),trim(region),trim(fn(ij)),
c    *ele,trim(units),trim(region),trim(fn(ij))


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 - Begin 2016 03 01:
c 901 format(
c    *'makecpt -Crainbow -T',f0.6,'/',f0.6,'/',f0.6,
c    *' -Z > temp.cpt')
c 901 format(
c    *'makecpt -Crainbow -T',f0.6,'/',f0.6,'/',f0.6,
c    *' > temp.cpt')
c - 2016 08 25
c 901 format(
c    *'makecpt -Crainbow -T',f0.10,'/',f0.10,'/',f0.10,
c    *' > temp.cpt')
c - 2016 09 07
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 %.2G')
c - 2016 09 07
  901 format(
     *'gmtset D_FORMAT %.12G',/,
     *'makecpt -Crainbow -T',f0.10,'/',f0.10,'/',f0.10,
     *' > temp.cpt',/,
     *'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 - Updated 2016 02 29:
c     write(99,902)fname(1:8),trim(suffixused)
  902 format(
     *'grdgradient ',a,'.',a,
     *'.grd -N4.8 -A90 -M -Gtempi.grd')
c - Until here 2016 02 29

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,igridsec,
     *trim(region),trim(fn(ij)),'method noise'


c - Commented out for 2016 02 29:
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 - Commented out for 2016 03 01:
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:
  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,
     *':."NADCON v5.0 ',a,' minus ',a,' ',a3,
     *'(',i0,' sec)',
     *1x,a,'-',a,1x,a,
     *':" -Ctemp.cpt -K > plot.ps')

c - GMT call to put a color scale on the map.  Add it
c - to file "plot.ps".

c - 2016 07 21:
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 -B/:',a,
     *': -O -K >> plot.ps')

c     write(99,904)trim(units)
c 904 format(
c    *'psscale -D4i/-0.4i/8i/0.2ih -Ctemp.cpt -I0.4 -B/:',a,
c    *': -O -K >> plot.ps')

c - GMT call to put contours on the map.

c - 2016 07 28...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
 
c - 2016 08 30:
      write(99,908)fname(1:8),trim(suffixused)
c - 2016 07 29:
c     write(99,908)fname(1:8),trim(suffixused),
c    *bw(ij),be(ij),bs(ij),bn(ij)
c - 2016 07 21:
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 - 2016 07 28 - Skip to here (skipping contours) for d3, 09 or "ete" plots
  807 continue

c - 2016 07 28 - Skip the coverage if this is a "d3" plot, 09 plot or "ete" plot *and*
c - we are at the overall map level of 1 ("entire")
      if(suffixused(ll-1:ll).eq.'d3' .and. ij.eq.1)goto 808
      if(suffixused(ll-1:ll).eq.'09' .and. ij.eq.1)goto 808
      if(fname(3:5).eq.'ete' .and. ij.eq.1)goto 808

c - GMT call to overlay the coverage
      write(99,909)trim(cvfname)

  909 format(
     *'psxy ',a,' -R -J -Sc0.02i -Gblack -O -K >> plot.ps')

  808 continue

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 - GMT call to convert the plot.ps file to plot.jpg.
      write(99,905)
  905 format(
     *'ps2raster plot.ps -Tj -P -A')

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 - Updated 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


  991 format(
     *'# -----------------------------------------------------',/,
     *'# color plots with coverage in ',a,1x,a,1x,a,1x,a,
     *   1x,' Grid:',a/,
     *'# -----------------------------------------------------',/,
     *'echo ...color plots with coverage in ',a,1x,a,1x,a,1x,a,
     *  1x,' Grid:',a)
c 991 format(
c    *'# -----------------------------------------------------',/,
c    *'# color plots with coverage in ',a,1x,a,1x,a,1x,a,/,
c    *'# -----------------------------------------------------',/,
c    *'echo ...color plots with coverage in ',a,1x,a,1x,a,1x,a)

      return

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine cpt	(	real*8	ave,
		real*8	std,
		real*8	csm,
		real*8	xlo,
		real*8	xhi,
		real*8	xin
	)

This subroutine generates the color pallette variables for a GMT color plot.

Parameters

[in]	ave	Average of the gridded data
[in]	std	Standard deviation of the gridded data
[in]	csm	Color Sigma Multiplier (how many sigmas on each side of the average do you want the colors to range?)
[out]	xlo	Low value
[out]	xhi	High value
[out]	xin	Interval

Changelog

2016 09 06:

Modified because the forcing of "scave" to be one non zero digit was throwing off the scalebar so far in Guam that the data in Guam wasn't even plotting. Change to make the interval still be one non zero digit, but then a simpler formula for the scaled average was put in.

2016 07 29:

Modified from original version to reflect "new math" invented this week that helps shrink the color bar and/or widen the color bar (see issues 14 and 15 in DRU-12, p. 48)

Definition at line 28 of file cpt.f.

Referenced by makeplotfiles02().

c - 2016 09 06
c   Modified because the forcing of "scave" to be one non zero
c   digit was throwing off the scalebar so far in Guam that
c   the data in Guam wasn't even plotting.  Change to make
c   the interval still be one non zero digit, but then a 
c   simpler formula for the scaled average was put in.

c - 2016 07 29
c   Modified from original version to reflect "new math" invented
c   this week that helps shrink the color bar and/or widen the
c   color bar (see issues 14 and 15 in DRU-12, p. 48)

      implicit none
      real*8 ave,std,csm,xlo,xhi,xin
      real*8 qv,qq,scave
      integer*4 iv,iq
      real*8 spread8,onzd2
c - This subroutine generates the color pallette variables
c - for a GMT color plot
c
c - Input:
c    ave = Average of the gridded data
c    std = Standard deviation of the gridded data
c    csm = Color Sigma Multiplier (how many sigmas on each
c          side of the average do you want the colors to range?)
c - Output:
c    xlo = Low value
c    xhi = High value
c    xin = Interval

      spread8 = 2 * csm * std / 8.d0
      xin = onzd2(spread8) 

c     iv = floor(dlog10(std))
c     qv = 10.d0**iv
c     qq = std / qv
c     iq = nint(qq)

c - 2016 09 06
c     write(6,*) ' cpt: ave (pre  onzd2) = ',ave
c     scave = onzd2(ave)
c     write(6,*) ' cpt: ave (post onzd2) = ',ave
      scave = xin * nint(ave / xin)

c     xin = iq * qv
c     scave = qv * nint(ave/qv)

c     xlo = scave - csm*xin
c     xhi = scave + csm*xin

      xlo = scave - 4*xin
      xhi = scave + 4*xin

      write(6,*) ' ---------------------'
      write(6,*) 'cpt: ave      = ',ave
      write(6,*) 'cpt: std      = ',std
      write(6,*) 'cpt: csm      = ',csm
      write(6,*) 'cpt: spread8  = ',spread8
      write(6,*) 'cpt: targ low = ',ave-csm*std
      write(6,*) 'cpt: targ high= ',ave+csm*std
      write(6,*) ' '
      write(6,*) 'cpt: xin      = ',xin
      write(6,*) 'cpt: scave    = ',scave
      write(6,*) 'cpt: xlo      = ',xlo  
      write(6,*) 'cpt: xhi      = ',xhi 

      return

Here is the caller graph for this function:

subroutine cpt2	(	real*8	med,
		real*8	csm,
		real*8	xlo,
		real*8	xhi,
		real*8	xin
	)

This subroutine generates the color pallette variables for a GMT color plot.

This particular routine is best for data that are all positive, but cluster near a small value while having a lot of outliers to the high-side. The color plot uses the MEDIAN (and a multiplier) to set the upper limit, while forcing the lower limit to be ZERO.

Parameters

[in]	med	Median of the gridded data
[in]	csm	Color Sigma Multiplier (The maximum color range will be based on csm*med. The minimum color range will be zero)
[out]	xlo	Low value
[out]	xhi	High value
[out]	xin	Interval

Definition at line 21 of file cpt2.f.

Referenced by makeplotfiles02(), and makeplotfiles03().

      implicit none
      real*8 med,csm,xlo,xhi,xin
      real*8 qv,qq,spread,std
      integer*4 iv,iq
c - This subroutine generates the color pallette variables
c - for a GMT color plot.  This particular routine is best
c - for data that are all positive, but cluster near a small
c - value while having a lot of outliers to the high-side.
c - The color plot uses the MEDIAN (and a multiplier) 
c - to set the upper limit, while forcing the lower limit
c - to be ZERO.
c
c - Input:
c    med = Median of the gridded data
c    csm = Color Sigma Multiplier (The maximum color range will
c          be based on csm*med.  The minimum color range will 
c          be zero)
c - Output:
c    xlo = Low value
c    xhi = High value
c    xin = Interval

      spread = csm*med
      std = spread/8.d0

      iv = floor(dlog10(std))
      qv = 10.d0**iv
      qq = std / qv
      iq = nint(qq)

      xin = iq * qv

      xlo = 0.d0
      xhi = 8 * xin
      return

Here is the caller graph for this function:

real function cubterp	(		x,
			f0,
			f1,
			f2,
			f3
	)

This function fits a cubic function through four points.

This function fits a cubic function through four equally spaced points along the x-axis at indices 0, 1, 2 and 3. The spacing along the x-axis is "dx"

Thus:

\begin{eqnarray*} f0 = f_0 &= y(x_0) \\ f1 = f_1 &= y(x_1) \\ f2 = f_2 &= y(x_2) \\ f3 = f_3 &= y(x_3) \end{eqnarray*}

Where:

\begin{eqnarray*} x_1 &= x_0 + dx \\ x_2 &= x_1 + dx \\ x_3 &= x_2 + dx \end{eqnarray*}

The input value is some value of "x" that falls between 0 and 3. The output value (cubterp) is the cubic function at x.

Parameters

[in]	x	Compute Interpolation at this positon, a value between 0 and 3 it is scaled relative to x_0 x_3 and dx. For example, a value of 1.5 is x_0 + 1.5*dx which falls between x1 and x2
[in]	f0	y value at x_0
[in]	f1	y value at x_1 = x_0 + dx
[in]	f2	y value at x_2 = x_0 + dx
[in]	f3	y value at x_3 = x_0 + dx`

This function uses Newton-Gregory forward polynomial

\begin{eqnarray*} \nabla f_0 &=& f_1 -f_0 \\ \nabla f_1 &=& f_2 -f_1 \\ \nabla^2 f_0 &=& \nabla f_1 - \nabla f_0 \\ cubterp(x, f_0, f_1, f_2, f_3) &=& f_0 + x \nabla f_0 + 0.5 x \left( x-1.0 \right) \nabla^2 f_0 \end{eqnarray*}

Definition at line 48 of file cubterp.f.

c - x = real*4
c - f0,f1,f2,f3 = real*4

c - This function fits a cubic function through
c - four points, *equally* spaced along the x-axis
c - at indices 0, 1, 2 and 3.  The spacing along the
c - x-axis is "dx"
c - Thus:
c - 
c -    f0 = y(x(0))
c -    f1 = y(x(1))
c -    f2 = y(x(2))
c -    f3 = y(x(3))
c -       Where
c -       x(1) = x(0) + dx
c -       x(2) = x(1) + dx
c -       x(3) = x(2) + dx

c - The input value is some value of "x" that falls
c - between 0 and 3.  The output value (cubterp) is
c - the cubic function at x.
c - 
c - This function uses Newton-Gregory forward polynomial

c     df0 =f1 -f0
c     df1 =f2 -f1
c     d2f0=df1-df0

c     qterp=f0 + x*df0 + 0.5*x*(x-1.0)*d2f0

      df0 = f1 - f0
      df1 = f2 - f1
      df2 = f3 - f2
     
      d2f0 = df1 - df0
      d2f1 = df2 - df1

      d3f0 = d2f1 - d2f0

      cubterp=f0 + x*df0 + 0.5*x*(x-1.0)*d2f0 
     *     +(1.0/6.0)*d3f0*x*(x-1.0)*(x-2.0)

      return

program decimate

(

)

Part of the NADCON5 NADCON5 Core Library , Extract a reduced (1 of n) dataset.

Decimate - Extract 1 of every n points

Program arguments

Arguments are newline terminated and read from standard input

When run from the command line, the program prints a prompt for each argument

They are enumerated here

Parameters

infile	Input File Name
outfile	Output File Name
ny	Latitude Decimation Ratio 1:ny
nx	Longitude Decimation Ratio 1:nx

Program Inputs:

• lin Input File A

Program Outputs:

• lout Output File to Write Decimated File

Definition at line 25 of file decimate.f.

References inouti(), and inoutr().

Here is the call graph for this function:

program gabs

(

)

Part of the NADCON5 NADCON5 Core Library , Convert values in a *.b grid to absolute value.

Belongs to the suite of ".b" file manipulators

This program will convert every value in a ".b" grid to its absolute value.

Program arguments

Arguments are newline terminated and read from standard input

When run from the command line, the program prints a prompt for each argument

They are enumerated here

Parameters

infile	Input File Name
outfile	Output File Name

Program Inputs:

• lin Input File (*.b grid)

Program Outputs:

• lout Output File (*.b grid)

Definition at line 26 of file gabs.f.

subroutine getgridbounds	(	character*10	region,
		real*8	xn,
		real*8	xs,
		real*8	xw,
		real*8	xe
	)

Subroutine to collect up the GRID boundaries for use in creating NADCON 5.

This CAN BE different than the MAP boundaries as such:

GRID boundaries will be just four values (n/s/w/e) for any region

MAP boundaries will allow multiple maps to be made and may or may not align with the GRID boundaries. Used to allow for more "close up" maps and such, without the need to screw up the MAP boundaries.

Parameters

[in]	region	Region Name
[out]	xn	north boundary for this region
[out]	xs	south boundary for this region
[out]	xw	west boundary for this region
[out]	xe	east boundary for this region

Subroutine Input Files:

• `Data/grid.parameters

Definition at line 26 of file getgridbounds.f.

Referenced by checkgrid(), makework(), mymedian5(), and myrms().

c
c - Subroutine to collect up the GRID boundaries
c - for use in creating NADCON 5
c
c - This CAN BE different than the MAP 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
c - Output:
c     xn,xs,xw,xe = N/S/W/E boundaries of the grid to be 
c                   created for this region

      character*10 region
      real*8    xn,xs,xw,xe
      character*80 card

      ifile = 3

      open(ifile,
     *file='Data/grid.parameters',
     *status='old',form='formatted')

c - Read 2 header lines
      read(ifile,'(a)')card
      read(ifile,'(a)')card

c - Now loop and find our region
    1 read(ifile,'(a)',end=2)card 
        if(trim(card(1:10)).eq.trim(region))then
          read(card(14:23),*)xn
          read(card(27:36),*)xs
          read(card(40:49),*)xw
          read(card(53:62),*)xe
          close(ifile)
          return
        endif
      goto 1

c - If we get here, the region sent to this subroutine wasn't found
c - in our file /Data/grid.parameters.  Crash and burn.
    2 write(6,100)trim(region)
      stop 10001 
  100 format(6x,'FATAL Subroutine getgridbounds: Unknown region: ',a)

Here is the caller graph for this function:

subroutine getmag	(		x,
			ix
	)

Subroutine to return the magnitude of a double precision value.

Parameters

[out]	x	result, magnitude of ix
[in]	ix	input douple precision

Definition at line 9 of file getmag.f.

c - Subroutine to return the magnitude of a double precision
c - value.
      implicit real*8 (a-h,o-z)
      y = dlog10(x)
      ix = floor(y) 
      write(6,*) ' getmag: x,y,ix = ',x,y,ix
      return

subroutine getmapbounds	(	character*1	mapflag,
			maxplots,
		character*10	region,
			nplots,
		character*10	olddtm,
		character*10	newdtm,
		real*8, dimension(maxplots)	bw,
		real*8, dimension(maxplots)	be,
		real*8, dimension(maxplots)	bs,
		real*8, dimension(maxplots)	bn,
		real*4, dimension(maxplots)	jm,
		real*4, dimension(maxplots)	b1,
		real*4, dimension(maxplots)	b2,
		character*10, dimension(maxplots)	fn,
		logical, dimension(maxplots)	lrv,
		real*8, dimension(maxplots)	rv0x,
		real*8, dimension(maxplots)	rv0y,
		real*8, dimension(maxplots)	rl0y
	)

Subroutine to collect up the MAP boundaries for use in creating NADCON 5.

This CAN BE different than the GRID boundaries as such:

GRID boundaries will be just four values (n/s/w/e) for any region

MAP boundaries will allow multiple maps to be made and may or may not align with the GRID boundaries. Used to allow for more "close up" maps and such, without the need to screw up the MAP boundaries.

Parameters

[in]	mapflag	Map Generation Flag
[in]	maxplots
[in]	region	region to get map bounds
[out]	nplots	number of plots generated
[in]	olddtm	source datum
[in]	newdtm	target datum
[out]	bw	western bound of plot (Array of length maxplots)
[out]	be	eastern bound of plot (Array of length maxplots)
[out]	bs	southern bound of plot (Array of length maxplots)
[out]	bn	northern bound of plot (Array of length maxplots)
[out]	jm	(Array of length maxplots)
[out]	b1	(Array of length maxplots)
[out]	b2	(Array of length maxplots)
[out]	fn	(Array of length maxplots)
[out]	lrv	(Array of length maxplots)
[out]	rv0x	(Array of length maxplots)
[out]	rv0y	(Array of length maxplots)
[out]	rl0y	(Array of length maxplots)

Version for NADCON 5 Built upon the original version used in GEOCON v2.0 Do not use with GEOCON v2.0

Broken down into sub-subroutines to make it easier to swap out when I make different choices.

Changelog

2016 08 29:

Taking in olddtm and newdtm now, and adding code to use that to filter out "Saint" regions in Alaska when plotting transformations not supported in those regions.

2016 08 26:

Used actual mercator projection math to compute the exact reference vector and label locations 1/2 inch and 3/4 inch respectively below the S/W corner of the plot.

2016 07 21:

Two new columns added to "map.parameters", which have the location of the reference vector. Return a logical "lrv" as true if there is an optional special location for the reverence vector. Return as false if not. If true, return lon/lat coords of ref vector origin in rv0x/rv0y. If false, return zeros in those fields.

Also, compute "jm" on the fly, ignoring what is in the table. All plots will now be forced PORTRAIT and forced no wider than 6" and no taller than 8", while maintaining proper X/Y ratios in a Mercator projection. That means, make the biggest plot possible, with the right ratio, that is neither wider than 6" nor taller than 8" and then, whatever the width of that largest plot is – return that width in the "jm" field.

Definition at line 73 of file getmapbounds.f.

Referenced by checkgrid(), makeplotfiles01(), makeplotfiles02(), and makeplotfiles03().

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='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

c - 2016 08 29:
c - Specific to the case where we are looking at region=alaska,
c - and NOT region=<some saint island>, then for the specific
c - case of transforming from NAD 27 forward, let's skip
c - the <saint islands>
        if(region.eq.'alaska' .and.
     *     iflag.eq.2 .and.
     *     iflag0.eq.2 .and.
     *     (card(16:25).eq.'stpaul' .or.
     *      card(16:25).eq.'stgeorge' .or.
     *      card(16:25).eq.'stlawrence' .or.
     *      card(16:25).eq.'stmatthew') .and.
     *      olddtm.eq.'nad27' )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('***************************',/,
     *  '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)

      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

Here is the caller graph for this function:

subroutine gridstats	(	character*200	fname,
			ave,
			std,
		real*8	med
	)

Subroutine to print grid statistics to stdout.

Parameters

[in]	fname	name of grid stat file
[out]	ave	average
[out]	std	standard deviatio
[out]	median

Definition at line 10 of file gridstats.f.

Referenced by makeplotfiles02(), and makeplotfiles03().

      implicit real*8(a-h,o-z)
      parameter(nmax = 1561*3541)
      real*8 med
      real*4 data(5000),s1,select2,arr(nmax)
      integer*4 nla,nlo,ikind
      character*200 fname
      open(11,file=fname,status='old',form='unformatted')

      ave = 0.d0
      rms = 0.d0

      read(11)glamn,glomn,dla,dlo,nla,nlo,ikind
      if(ikind.eq.0)stop 20304

      ikt = 0
      do 1 i=1,nla
        if(ikind.eq.1)read(11)(data(j),j=1,nlo)
        do 3 j=1,nlo
          ikt = ikt + 1
          arr(ikt) = data(j)
    3   continue
        do 2 j=1,nlo
          ave = ave + dble(data(j))
          rms = rms + dble(data(j))**2
    2   continue
    1 continue
      fact = dble(nla*nlo) / dble(nla*nlo - 1)

      ave = ave / dble(nla*nlo)
      rms = sqrt(rms / dble(nla*nlo))
      std = dsqrt(fact*(rms**2 - ave**2))

      npts = ikt
      nmed = 1+ ikt/2

      write(6,*) ' IN gridstats, npts = ',npts
      write(6,*) ' IN gridstats, nmed = ',nmed

c     nmed = ikt/3
      s1 = select2(nmed,npts,arr,nmax)

      write(6,*) ' IN gridstats, s1   = ',s1
      med = s1
      write(6,*) ' IN gridstats, med  = ',med

      write(6,100)trim(fname),ave,std,med
  100 format('gridstats for ',a,' = ',3f20.10)

      close(11)

      return

Here is the caller graph for this function:

program gscale

(

)

Part of the NADCON5 NADCON5 Core Library , Scales a grid by a factor.

Program arguments

Arguments are newline terminated and read from standard input

When run from the command line, the program prints a prompt for each argument

They are enumerated here

Parameters

infile	Input File Name
factor	Scaling Factor
outfile	Output File Name

Program Inputs:

• lin1 Input File

Program Outputs:

• lout Output File

Definition at line 23 of file gscale.f.

program gsqr

(

)

Part of the NADCON5 NADCON5 Core Library , Squares values in a *.b grid.

Belongs to the suite of ".b" file manipulators

This program will convert every value in a ".b" grid to its squared value.

Program arguments

Arguments are newline terminated and read from standard input

When run from the command line, the program prints a prompt for each argument

They are enumerated here

Parameters

infile	Input File Name
outfile	Output File Name

Program Inputs:

• lin Input File (*.b grid)

Program Outputs:

• lout Output File (*.b grid)

Definition at line 26 of file gsqr.f.

program gsqrt

(

)

Part of the NADCON5 NADCON5 Core Library , Square Root of values in a *.b grid.

Belongs to the suite of ".b" file manipulators

This program will convert every value in a ".b" grid to its square-root value.

Program arguments

Arguments are newline terminated and read from standard input

When run from the command line, the program prints a prompt for each argument

They are enumerated here

Parameters

infile	Input File Name
outfile	Output File Name

Program Inputs:

• lin Input File (*.b grid)

Program Outputs:

• lout Output File (*.b grid)

Definition at line 26 of file gsqrt.f.

subroutine indexxd	(	integer	n,
		integer	nd,
		real*8, dimension(nd)	arr,
		integer, dimension(nd)	indx
	)

Subroutine to perform ?? indexing on floating point data (double precision)

Parameters

[in]	n	number of iterations (rows?)
[in]	nd	array and index dimensions
[in]	arr	input data array
[out]	indx	index out

Changelog

1/8/2004:

Modified to allow indx and arr to be DIMENSIONED differently than the number of good values they contain

11/7/2003

Modified to REAL*8 by D. Smith,

Definition at line 20 of file indexxd.for.

Referenced by getmedian().

c - Modified to REAL*8 by D. Smith, 11/7/2003
c - Further modified 1/8/2004 to allow indx and arr to
c - be DIMENSIONED differently than the number of good
c - values they contain

      integer nd
      INTEGER n,indx(nd),m,nstack
      REAL*8 arr(nd)
      parameter(m=7,nstack=50)
      INTEGER i,indxt,ir,itemp,j,jstack,k,l,istack(nstack)
      REAL*8 a
      do 11 j=1,n
        indx(j)=j
11    continue
      jstack=0
      l=1
      ir=n
1     if(ir-l.lt.m)then
        do 13 j=l+1,ir
          indxt=indx(j)
          a=arr(indxt)
          do 12 i=j-1,1,-1
            if(arr(indx(i)).le.a)goto 2
            indx(i+1)=indx(i)
12        continue
          i=0
2         indx(i+1)=indxt
13      continue
        if(jstack.eq.0)return
        ir=istack(jstack)
        l=istack(jstack-1)
        jstack=jstack-2
      else
        k=(l+ir)/2
        itemp=indx(k)
        indx(k)=indx(l+1)
        indx(l+1)=itemp
        if(arr(indx(l+1)).gt.arr(indx(ir)))then
          itemp=indx(l+1)
          indx(l+1)=indx(ir)
          indx(ir)=itemp
        endif
        if(arr(indx(l)).gt.arr(indx(ir)))then
          itemp=indx(l)
          indx(l)=indx(ir)
          indx(ir)=itemp
        endif
        if(arr(indx(l+1)).gt.arr(indx(l)))then
          itemp=indx(l+1)
          indx(l+1)=indx(l)
          indx(l)=itemp
        endif
        i=l+1
        j=ir
        indxt=indx(l)
        a=arr(indxt)
3       continue
          i=i+1
        if(arr(indx(i)).lt.a)goto 3
4       continue
          j=j-1
        if(arr(indx(j)).gt.a)goto 4
        if(j.lt.i)goto 5
        itemp=indx(i)
        indx(i)=indx(j)
        indx(j)=itemp
        goto 3
5       indx(l)=indx(j)
        indx(j)=indxt
        jstack=jstack+2
        if(jstack.gt.nstack)pause 'NSTACK too small in indexx'
        if(ir-i+1.ge.j-l)then
          istack(jstack)=ir
          istack(jstack-1)=i
          ir=j-1
        else
          istack(jstack)=j-1
          istack(jstack-1)=l
          l=i
        endif
      endif
      goto 1

Here is the caller graph for this function:

subroutine indexxi	(	integer	n,
		integer	nd,
		integer*4, dimension(nd)	arr,
		integer, dimension(nd)	indx
	)

Subroutine to perform ?? indexing on integer data.

Parameters

[in]	n	number of iterations (rows?)
[in]	nd	array and index dimensions
[in]	arr	input data array
[out]	indx	index out

Changelog

2/5/2013:

Modified by D. Smith. Arr has been changed to integer*4. And like other versions of "indexx" which I've modified, I allow indx and arr to be DIMENSIONED differently than the number of good values they contain

Parameters

nd	Modified by D. Smith, 2/5/2013 . Arr has been changed to integer*4. And like other versions of "indexx" which I've modified, I allow indx and arr to be DIMENSIONED differently than the number of good values they contain

Definition at line 20 of file indexxi.for.

Referenced by mymedian5(), and myrms().

c> Modified by D. Smith, 2/5/2013 .  Arr has been
c> changed to integer*4.  And like other versions
c> of "indexx" which I've modified, 
c> I allow indx and arr to
c> be DIMENSIONED differently than the number of good
c> values they contain

      integer nd
      INTEGER n,indx(nd),m,nstack
cDAS  REAL*4 arr(nd)
      integer*4 arr(nd)
      parameter(m=7,nstack=50)
      INTEGER i,indxt,ir,itemp,j,jstack,k,l,istack(nstack)
cDAS  REAL*4 a     
      integer*4 a 

      do 11 j=1,n
        indx(j)=j
11    continue
      jstack=0
      l=1
      ir=n
1     if(ir-l.lt.m)then
        do 13 j=l+1,ir
          indxt=indx(j)
          a=arr(indxt)
          do 12 i=j-1,1,-1
            if(arr(indx(i)).le.a)goto 2
            indx(i+1)=indx(i)
12        continue
          i=0
2         indx(i+1)=indxt
13      continue
        if(jstack.eq.0)return
        ir=istack(jstack)
        l=istack(jstack-1)
        jstack=jstack-2
      else
        k=(l+ir)/2
        itemp=indx(k)
        indx(k)=indx(l+1)
        indx(l+1)=itemp
        if(arr(indx(l+1)).gt.arr(indx(ir)))then
          itemp=indx(l+1)
          indx(l+1)=indx(ir)
          indx(ir)=itemp
        endif
        if(arr(indx(l)).gt.arr(indx(ir)))then
          itemp=indx(l)
          indx(l)=indx(ir)
          indx(ir)=itemp
        endif
        if(arr(indx(l+1)).gt.arr(indx(l)))then
          itemp=indx(l+1)
          indx(l+1)=indx(l)
          indx(l)=itemp
        endif
        i=l+1
        j=ir
        indxt=indx(l)
        a=arr(indxt)
3       continue
          i=i+1
        if(arr(indx(i)).lt.a)goto 3
4       continue
          j=j-1
        if(arr(indx(j)).gt.a)goto 4
        if(j.lt.i)goto 5
        itemp=indx(i)
        indx(i)=indx(j)
        indx(j)=itemp
        goto 3
5       indx(l)=indx(j)
        indx(j)=indxt
        jstack=jstack+2
        if(jstack.gt.nstack)pause 'NSTACK too small in indexx'
        if(ir-i+1.ge.j-l)then
          istack(jstack)=ir
          istack(jstack-1)=i
          ir=j-1
        else
          istack(jstack)=j-1
          istack(jstack-1)=l
          l=i
        endif
      endif
      goto 1

Here is the caller graph for this function:

integer*2 function iselect2	(	integer	k,
		integer	n,
		integer*2, dimension(nmax)	arr,
		integer	nmax
	)

Function to select an element of a partially filled, but packed multi dimensional array, integer*2

Finds the "kth" element of an array, "arr", which is dimensioned to be "nmax" values long, but which only has data in the first "n" cells.

Changelog

1/14/2016:

Like "select2" but modified by D. Smith to allow an "nmax" array given, but which only has values in elements 1-n, and to have "arr" be Integer*2

Definition at line 18 of file iselect2.for.

c - Like "select2" but modified by D. Smith on 1/14/2016
c - to allow an "nmax" array given, but which only
c - has values in elements 1-n, and to have "arr"
c - be Integer*2

      INTEGER k,n,nmax
c     REAL*8 select,arr(nmax)
      integer*2 iselect2,arr(nmax)
      INTEGER i,ir,j,l,mid
c     REAL*8 a,temp
      integer*2 a,temp
      l=1
      ir=n
1     if(ir-l.le.1)then
        if(ir-l.eq.1)then
          if(arr(ir).lt.arr(l))then
            temp=arr(l)
            arr(l)=arr(ir)
            arr(ir)=temp
          endif
        endif
        iselect2=arr(k)
        return
      else
        mid=(l+ir)/2
        temp=arr(mid)
        arr(mid)=arr(l+1)
        arr(l+1)=temp
        if(arr(l+1).gt.arr(ir))then
          temp=arr(l+1)
          arr(l+1)=arr(ir)
          arr(ir)=temp
        endif
        if(arr(l).gt.arr(ir))then
          temp=arr(l)
          arr(l)=arr(ir)
          arr(ir)=temp
        endif
        if(arr(l+1).gt.arr(l))then
          temp=arr(l+1)
          arr(l+1)=arr(l)
          arr(l)=temp
        endif
        i=l+1
        j=ir
        a=arr(l)
3       continue
          i=i+1
        if(arr(i).lt.a)goto 3
4       continue
          j=j-1
        if(arr(j).gt.a)goto 4
        if(j.lt.i)goto 5
        temp=arr(i)
        arr(i)=arr(j)
        arr(j)=temp
        goto 3
5       arr(l)=arr(j)
        arr(j)=a
        if(j.ge.k)ir=j-1
        if(j.le.k)l=i
      endif
      goto 1

real*4 function onzd

(

real*4

x

)

Function to round a digit to one significant figure (one non zero digit), single precision.

Function "onzd" stands for "One Non Zero Digit"

It takes a Real*4 number as input, and rounds that number to the closest number containing only 1 non-zero digit. The list of such numbers is infifinite, but contain these, in order:

0.7 , 0.8 , 0.9 , 1 , 2 , 3 , 4 , 5 , 6 , 7 , 8 , 9, 10 , 20 , 30 , etc etc

Parameters

[in]

x

input value

Returns

real*4 rounded value of x to one non zero digit

Examples of input/output are:

   0.000019      =>    0.000020
   0.007432      =>    0.007000
   1.7           =>    2.000000
   9.143         =>    9.000000
  17.4           =>   20.000000
 947.3           =>  900.000000
 987.432         => 1000.000000
1014.8           => 1000.000000
1502.7           => 2000.000000

Definition at line 29 of file onzd.f.

      implicit none
      real*4 onzd,x,y
      real*8 q,qten
      integer*4 imag,iq,isign
      
c - Function "onzd" stands for "One Non Zero Digit"  

c - It takes a Real*4 number as input, and rounds that 
c - number to the closest number containing only 1 non-zero digit.
c - The list of such numbers is infifinite, but contain
c - these, in order:
c -    0.7 , 0.8 , 0.9 , 1 , 2 , 3 , 4 , 5 , 6 , 7 , 8 , 9, 10 , 20 , 30 , etc etc
c - 
c - Examples of input/output are:
c -       0.000019      =>    0.000020
c -       0.007432      =>    0.007000
c -       1.7           =>    2.000000
c -       9.143         =>    9.000000
c -      17.4           =>   20.000000
c -     947.3           =>  900.000000
c -     987.432         => 1000.000000
c -    1014.8           => 1000.000000
c -    1502.7           => 2000.000000

         isign = +1
         y = x
         if(x.lt.0)then
           isign = -1
           y = -x
         endif

c - 1) Determine magnitude of x, in terms of integer exponent of ten
        imag = floor(log10(y))

c - 1a) Get the multiplier
        qten = (10.d0**imag)

c - 2) Get into a range between 0.0 and 10.0
        q = dble(y) / qten

c - 3) Round to the closest integer (0 through 10)
        iq = nint(q)

c - 4) Scale back to the original size 
        onzd = isign * dble(iq) * qten

c       write(6,100)x,imag,qten,q,iq,onzd
c 100   format(f30.15,1x,i10,1x,f30.15,1x,f30.15,1x,i10,1x,f30.15)
    
        return

real*8 function onzd2

(

real*8

x

)

Function to round a digit to one significant figure (one non zero digit), double precision.

Function "onzd" stands for "One Non Zero Digit"

Version 2 operates just like version 1 (onzd()), only the input and output will be real*8 values, not real*4.

It takes a Real*8 number as input, and rounds that number to the closest number containing only 1 non-zero digit. The list of such numbers is infifinite, but contain these, in order:

0.7 , 0.8 , 0.9 , 1 , 2 , 3 , 4 , 5 , 6 , 7 , 8 , 9, 10 , 20 , 30 , etc etc

Parameters

[in]

x

input value

Returns

real*8 rounded value of x to one non zero digit

Examples of input/output are:

   0.000019      =>    0.000020
   0.007432      =>    0.007000
   1.7           =>    2.000000
   9.143         =>    9.000000
  17.4           =>   20.000000
 947.3           =>  900.000000
 987.432         => 1000.000000
1014.8           => 1000.000000
1502.7           => 2000.000000

Definition at line 32 of file onzd2.f.

Referenced by checkgrid(), makeplotfiles01(), makeplotfiles02(), makeplotfiles03(), and myrms().

      implicit none
      real*8 onzd2,x,y
      real*8 q,qten
      integer*4 imag,iq,isign
      
c - Function "onzd2" stands for "One Non Zero Digit, version 2"  

c - Version 2 operates just like version 1 (onzd.f), only the input
c - and output will be real*8 values, not real*4.

c - It takes a Real*8 number as input, and rounds that 
c - number to the closest number containing only 1 non-zero digit.
c - The list of such numbers is infifinite, but contain
c - these, in order:
c -    0.7 , 0.8 , 0.9 , 1 , 2 , 3 , 4 , 5 , 6 , 7 , 8 , 9, 10 , 20 , 30 , etc etc
c - 
c - Examples of input/output are:
c -       0.000019      =>    0.000020
c -       0.007432      =>    0.007000
c -       1.7           =>    2.000000
c -       9.143         =>    9.000000
c -      17.4           =>   20.000000
c -     947.3           =>  900.000000
c -     987.432         => 1000.000000
c -    1014.8           => 1000.000000
c -    1502.7           => 2000.000000

        y = x
        isign = +1
        if(x.lt.0)then
          isign = -1
          y = -x
        endif

c - 1) Determine magnitude of x, in terms of integer exponent of ten
        imag = floor(dlog10(y))

c - 1a) Get the multiplier
        qten = (10.d0**imag)

c - 2) Get into a range between 0.0 and 10.0
        q = dble(y) / qten

c - 3) Round to the closest integer (0 through 10)
        iq = nint(q)

c - 4) Scale back to the original size 
        onzd2 = isign * dble(iq) * qten

c       write(6,100)x,imag,qten,q,iq,onzd2
c 100   format(f30.15,1x,i10,1x,f30.15,1x,f30.15,1x,i10,1x,f30.15)
    
        return

Here is the caller graph for this function:

subroutine plotcoast	(	character*10	region,
			ifnum
	)

Subroutine to write GMT-based commands to create a shoreline Write GMT-based commands to create a shoreline based on region.

Use GMT-default coastline for:

• conus
• alaska

Use Dru's custome coastline for:

• hawaii
• prvi
• as
• guamcnmi
• stlawrence
• stmatthew
• stgeorge
• stpaul

Parameters

[in]	region	The Region to create coastline for
[in]	ifnum	the file descriptor of the output file to write GMT commands to

Changelog

2016 01 07:

Forced the Alaska region to plot the islands of St. George, St. Matthew and St. Paul (St. Lawrence is already plotted), as well as 35 missing Aleutian Islands

2015 09 23:

Added four new regions:

• St. Lawrence Island, Alaska
• St. Matthew Island, Alaska
• St. George Island, Alaska
• St. Paul Island, Alaska

Definition at line 40 of file plotcoast.f.

Referenced by bwplotcv(), bwplotvc(), coplot(), and coplotwcv().

c - 2016 01 07
c     Forced the Alaska region to plot the
c     islands of St. George, St. Matthew and St. Paul
c     (St. Lawrence is already plotted), as well
c     as 35 missing Aleutian Islands

c - 2015 09 23:  Added four new regions:
c      St. Lawrence Island, Alaska
c      St. Matthew Island, Alaska
c      St. George Island, Alaska
c      St. Paul Island, Alaska

      character*10 region

c - Write GMT-based commands to create a shoreline
c - based on region.
c - Use GMT-default coastline for:
c      conus
c      alaska
c - Use Dru's custome coastline for:
c      hawaii
c      prvi
c      as
c      guamcnmi
c      stlawrence
c      stmatthew 
c      stgeorge  
c      stpaul    


      if    (trim(region).eq.'conus')then
        write(ifnum,100) 
      elseif(trim(region).eq.'alaska')then
        write(ifnum,109) 
      elseif(trim(region).eq.'hawaii')then
        write(ifnum,101) 
      elseif(trim(region).eq.'prvi')then
        write(ifnum,102) 
      elseif(trim(region).eq.'guamcnmi')then
        write(ifnum,103) 
      elseif(trim(region).eq.'as')then
        write(ifnum,104) 
      elseif(trim(region).eq.'stlawrence')then
        write(ifnum,105) 
      elseif(trim(region).eq.'stmatthew')then
        write(ifnum,106) 
      elseif(trim(region).eq.'stgeorge')then
        write(ifnum,107) 
      elseif(trim(region).eq.'stpaul')then
        write(ifnum,108) 
     
      else
        write(6,200)region
        stop 10001
      endif
      return

  200 format('FATAL in plotcoast.f: Unknown region: ',a10)

c - CONUS 
  100 format(
     *'pscoast -Df -R -JM -W0.25p -N1 -N2 -A1200 -O >>  plot.ps')

c - HAWAII
  101 format(
     *'psxy Boundaries/HI_Hawaii.gmt -R -JM -B -O -K',
     *' >> plot.ps',/,
     *'psxy Boundaries/HI_Kahoolawe.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/HI_Kauai.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/HI_Lanai.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/HI_Maui.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/HI_Molokai.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/HI_Nihau.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/HI_Oahu.gmt -R -JM -B -O ',
     *' >> plot.ps')

c - PRVI
  102 format(
     *'psxy Boundaries/VQ_StCroix.gmt -R -JM -B -O -K',
     *' >> plot.ps',/,
     *'psxy Boundaries/VQ_StJohn.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/VQ_StThomas.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/PR_Culebra.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/PR_Desecheo.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/PR_Mona.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/PR_Viequez.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/PR_PuertoRico.gmt -R -JM -B -O ',
     *' >> plot.ps')
 
c - GUAM/CNMI
  103 format(
     *'psxy Boundaries/GQ_Guam.gmt -R -JM -B -O -K',
     *' >> plot.ps',/,
     *'psxy Boundaries/CQ_Rota.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/CQ_Saipan.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/CQ_Tinian.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/CQ_Anatahan.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/CQ_Sarigan.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/CQ_Guguan.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/CQ_Alamagan.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/CQ_Pagan.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/CQ_Agrihan.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/CQ_Asuncion.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/CQ_MaugWest.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/CQ_MaugEast.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/CQ_MaugNorth.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/CQ_Pajaros.gmt -R -JM -B -O',
     *' >> plot.ps')

c - AS  
  104 format(
     *'psxy Boundaries/AS_Aunuu.gmt -R -JM -B -O -K',
     *' >> plot.ps',/,
     *'psxy Boundaries/AS_Nuutele.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/AS_Ofu.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/AS_Olosega.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/AS_Rose.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/AS_Swains.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/AS_Tau.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/AS_Tutila.gmt -R -JM -B -O ',
     *' >> plot.ps')

c - St. Lawrence Island, Alaska
  105 format(
     *'psxy Boundaries/AK_StLawrence.gmt -R -JM -B -O -K',
     *' >> plot.ps',/,
     *'psxy Boundaries/AK_NorthPunuk.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/AK_MiddlePunuk.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/AK_SouthPunuk.gmt -R -JM -B -O ',
     *' >> plot.ps')

c - St. Matthew Island, Alaska
  106 format(
     *'psxy Boundaries/AK_StMatthew.gmt -R -JM -B -O -K',
     *' >> plot.ps',/,
     *'psxy Boundaries/AK_Hall.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/AK_Pinnacle.gmt -R -JM -B -O ',
     *' >> plot.ps')

c - St. George Island, Alaska
  107 format(
     *'psxy Boundaries/AK_StGeorge.gmt -R -JM -B -O',
     *' >> plot.ps')

c - St. Paul Island, Alaska
  108 format(
     *'psxy Boundaries/AK_StPaul.gmt -R -JM -B -O -K',
     *' >> plot.ps',/,
     *'psxy Boundaries/AK_Otter.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/AK_Walrus.gmt -R -JM -B -O ',
     *' >> plot.ps')

c - ALASKA
  109 format(
     *'pscoast -Df -R -JM -W0.25p -N1 -N2 -A1200 -O -K ',
     *' >>  plot.ps',/,
     *'psxy Boundaries/AK_NorthPunuk.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/AK_MiddlePunuk.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/AK_SouthPunuk.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/AK_StMatthew.gmt -R -JM -B -O -K',
     *' >> plot.ps',/,
     *'psxy Boundaries/AK_Hall.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/AK_Pinnacle.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/AK_StGeorge.gmt -R -JM -B -O -K',
     *' >> plot.ps',/,
     *'psxy Boundaries/AK_StPaul.gmt -R -JM -B -O -K',
     *' >> plot.ps',/,
     *'psxy Boundaries/AK_Otter.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/AK_Walrus.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/AK_Adak.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/AK_Agattu.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/AK_Amchitka.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/AK_Amlia.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/AK_Amukta.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/AK_Andreanof.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/AK_Atka.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/AK_Attu.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/AK_Buldir.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/AK_Carlisle.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/AK_Chagulak.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/AK_Chugul.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/AK_Fenimore.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/AK_FourMountains.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/AK_GreatSitkin.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/AK_Herbert.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/AK_Igitkin.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/AK_Ikiginak.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/AK_Kagalaska.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/AK_Kagami.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/AK_Kanaga.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/AK_Kasatochi.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/AK_Kiska.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/AK_Koniuji.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/AK_Oglodak.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/AK_Seguam.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/AK_Semichi.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/AK_Semisopochnoi.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/AK_Shemya.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/AK_Tagalak.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/AK_Tanaga.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/AK_Ulak.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/AK_Uliaga.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/AK_Umak.gmt -R -JM -B -O -K ',
     *' >> plot.ps',/,
     *'psxy Boundaries/AK_Yunaska.gmt -R -JM -B -O ',
     *' >> plot.ps')

Here is the caller graph for this function:

real function qterp	(		x,
			f0,
			f1,
			f2
	)

This function fits a quadratic function through 3 points.

This function fits a parabola (quadratic) function through three equally spaced points along the x-axis at indices 0, 1, and 2. The spacing along the x-axis is "dx"

Thus:

\begin{eqnarray*} f0 = f_0 &= y(x_0) \\ f1 = f_1 &= y(x_1) \\ f2 = f_2 &= y(x_2) \end{eqnarray*}

Where:

\begin{eqnarray*} x_1 &= x_0 + dx \\ x_2 &= x_1 + dx \\ x_3 &= x_2 + dx \end{eqnarray*}

The input value is some value of "x" that falls between 0 and 2. The output value (qterp) is the quadratic function at x.

Parameters

[in]	x	Compute Interpolation at this positon, a value between 0 and 3 it is scaled relative to x_0 x_2 and dx. For example, the value of 1.5 is x_0 + 1.5*dx which falls between x1 and x2
[in]	f0	y value at x_0
[in]	f1	y value at x_1 = x_0 + dx
[in]	f2	y value at x_2 = x_0 + dx

Returns

real quadratically interpolated value of f(x*) where x* = x_0 + x*dxx

This function uses Newton-Gregory forward polynomial

\begin{eqnarray*} \nabla f_0 &=& f_1 -f_0 \\ \nabla f_1 &=& f_2 -f_1 \\ \nabla^2 f_0 &=& \nabla f_1 - \nabla f_0 \\ qterp(x, f_0, f_1, f_2) &=& f_0 + x \nabla f_0 + 0.5 x \left( x-1.0 \right) \nabla^2 f_0 \end{eqnarray*}

Definition at line 47 of file qterp.f.

c - x = real*4
c - f0,f1,f2 = real*4

c - This function fits a parabola (quadratic) through
c - three points, *equally* spaced along the x-axis
c - at indices 0, 1 and 2.  The spacing along the
c - x-axis is "dx"
c - Thus:
c - 
c -    f0 = y(x(0))
c -    f1 = y(x(1))
c -    f2 = y(x(2))
c -       Where
c -       x(1) = x(0) + dx
c -       x(2) = x(1) + dx

c - The input value is some value of "x" that falls
c - between 0 and 2.  The output value (qterp2) is
c - the parabolic function at x.
c - 
c - This function uses Newton-Gregory forward polynomial

      df0 =f1 -f0
      df1 =f2 -f1
      d2f0=df1-df0

      qterp=f0 + x*df0 + 0.5*x*(x-1.0)*d2f0

      return

program regrd2

(

)

Part of the NADCON5 NADCON5 Core Library , regrid data.

Regrid gridded data using biquadratic interpolation

Program arguments

Arguments are newline terminated and read from standard input

When run from the command line, the program prints a prompt for each argument

They are enumerated here

Parameters

infile	Input Master Grid File Name
outfile	Output Regrid File Name
nrow	Number of rows in new Grid (Lat)
ncol	Number of cols in new Grid (Lon)

Program Inputs:

• lin Input File

Program Outputs:

• lout Output File

Definition at line 25 of file regrd2.f.

References bquad1(), bquad2(), r1(), and r2().

Here is the call graph for this function:

real*8 function select2	(	integer	k,
		integer	n,
		real*8, dimension(nmax)	arr,
		integer	nmax
	)

Function to select an element of a partially filled, but packed multi dimensional array, double precision.

Finds the "kth" element of an array, "arr", which is dimensioned to be "nmax" values long, but which only has data in the first "n" cells.

Changelog

7/17/2008:

Like "select2" but modified by D. Smith to allow an "nmax" array given, but which only has values in elements 1-n, and to have "arr" be Integer*2

Definition at line 18 of file select2_dbl.for.

c - Like "select" but modified by D. Smith on 7/17/2008
c - to allow an "nmax" array given, but which only
c - has values in elements 1-n

      INTEGER k,n,nmax
      REAL*8 select2,arr(nmax)
      INTEGER i,ir,j,l,mid
      REAL*8 a,temp
      l=1
      ir=n
1     if(ir-l.le.1)then
        if(ir-l.eq.1)then
          if(arr(ir).lt.arr(l))then
            temp=arr(l)
            arr(l)=arr(ir)
            arr(ir)=temp
          endif
        endif
        select2=arr(k)
        return
      else
        mid=(l+ir)/2
        temp=arr(mid)
        arr(mid)=arr(l+1)
        arr(l+1)=temp
        if(arr(l+1).gt.arr(ir))then
          temp=arr(l+1)
          arr(l+1)=arr(ir)
          arr(ir)=temp
        endif
        if(arr(l).gt.arr(ir))then
          temp=arr(l)
          arr(l)=arr(ir)
          arr(ir)=temp
        endif
        if(arr(l+1).gt.arr(l))then
          temp=arr(l+1)
          arr(l+1)=arr(l)
          arr(l)=temp
        endif
        i=l+1
        j=ir
        a=arr(l)
3       continue
          i=i+1
        if(arr(i).lt.a)goto 3
4       continue
          j=j-1
        if(arr(j).gt.a)goto 4
        if(j.lt.i)goto 5
        temp=arr(i)
        arr(i)=arr(j)
        arr(j)=temp
        goto 3
5       arr(l)=arr(j)
        arr(j)=a
        if(j.ge.k)ir=j-1
        if(j.le.k)l=i
      endif
      goto 1

real*4 function select2	(	integer	k,
		integer	n,
		real*4, dimension(nmax)	arr,
		integer	nmax
	)

Function to select an element of a partially filled, but packed multi dimensional array, single precision.

Finds the "kth" element of an array, "arr", which is dimensioned to be "nmax" values long, but which only has data in the first "n" cells.

Changelog

7/17/2008:

Like "select2" but modified by D. Smith to allow an "nmax" array given, but which only has values in elements 1-n, and to have "arr" be Integer*2

Definition at line 18 of file select2_mod.for.

c - Like "select" but modified by D. Smith on 7/17/2008
c - to allow an "nmax" array given, but which only
c - has values in elements 1-n.  Maintain the Real*4 nature of arr.

c - Finds the "kth" element of an array, "arr", which
c - is dimensioned to be "nmax" values long, but which
c - only has data in the first "n" cells.

      INTEGER k,n,nmax
      REAL*4 select2,arr(nmax)
      INTEGER i,ir,j,l,mid
      REAL*4 a,temp
      l=1
      ir=n
1     if(ir-l.le.1)then
        if(ir-l.eq.1)then
          if(arr(ir).lt.arr(l))then
            temp=arr(l)
            arr(l)=arr(ir)
            arr(ir)=temp
          endif
        endif
        select2=arr(k)
        return
      else
        mid=(l+ir)/2
        temp=arr(mid)
        arr(mid)=arr(l+1)
        arr(l+1)=temp
        if(arr(l+1).gt.arr(ir))then
          temp=arr(l+1)
          arr(l+1)=arr(ir)
          arr(ir)=temp
        endif
        if(arr(l).gt.arr(ir))then
          temp=arr(l)
          arr(l)=arr(ir)
          arr(ir)=temp
        endif
        if(arr(l+1).gt.arr(l))then
          temp=arr(l+1)
          arr(l+1)=arr(l)
          arr(l)=temp
        endif
        i=l+1
        j=ir
        a=arr(l)
3       continue
          i=i+1
        if(arr(i).lt.a)goto 3
4       continue
          j=j-1
        if(arr(j).gt.a)goto 4
        if(j.lt.i)goto 5
        temp=arr(i)
        arr(i)=arr(j)
        arr(j)=temp
        goto 3
5       arr(l)=arr(j)
        arr(j)=a
        if(j.ge.k)ir=j-1
        if(j.le.k)l=i
      endif
      goto 1

program subtrc

(

)

Part of the NADCON5 NADCON5 Core Library , Subtract one grid from another.

Program arguments

Arguments are newline terminated and read from standard input

When run from the command line, the program prints a prompt for each argument

They are enumerated here

Parameters

infileA	First Input File Name
infileB	Second Input File Name
outfile	Output File Name of A*B

Program Inputs:

• lin1 Input File A
• lin2 Input File B
• lout Output File to Write A-B

Definition at line 21 of file subtrc.f.

subroutine vecstats	(	character*200	fname,
		integer*4	n
	)

Subroutine to tell us how many thinned vectors were used to make a grid.

Parameters

[in]	fname	vector filename to read
[out]	n	number of thinned vectors

Definition at line 8 of file vecstats.f.

Referenced by makeplotfiles02().

c - Subroutine to tell us how many thinned vectors were
c - used to make a grid
      character*200 fname
      integer*4 n
      character*80 card
      open(90,file=fname,status='old',form='formatted')
      n = 0
    1 read(90,'(a)',end=2)card
        n = n + 1
      goto 1
    2 return

Here is the caller graph for this function:

program xyz2b

(

)

Part of the NADCON5 NADCON5 Core Library , Converts GMT *.grd to a *.b NADCON style grid file.

Turn gmt/netcdf grd dump into my grid file (real number version) assumes grd dump is longitude/latitude/real (binary s.p.)

Program arguments

Arguments are newline terminated and read from standard input

When run from the command line, the program prints a prompt for each argument

They are enumerated here

Parameters

infile	Input File Name
outfile	Output File Name

Program Inputs:

• lin Input File (*.grd)

Program Outputs:

• lout Output File (*.b)

Definition at line 25 of file xyz2b.f.

References put1(), put2(), w1(), and w2().

Here is the call graph for this function: