c  utility.for
c
c             2-d grid utility program
c
c  format:
c  header record contains:
c    56 charater title, ( 14 * 4 character words)
c     8 character creation program or operation identifier
c    number of columns,  the number of words in a row record minus one
c    number of rows,     the number of records after the header.
c    number of z values associated with each grid point ( always one ).
c    xo,  coordinate (in real data units ) along profile of the rows.
c    dx   spacing    (d.u.) along a row.
c    yo   coordinate (d.u.) for the transverse position of a row.
c    dy   spacing    (d.u.) of the rows.

c   followed by row records consisting of a y coordinate (unused if
c   dely is nonzero) and grid value array ncolumns in length.
c
c     Mike Webring, USGS, branch of geophysics
c     Converted to PC with minor changes by Lin Cordell Aug 90
c     The plug option doesn't work and has been disabled, use MEGAPLUG
c     Converted to ask and iosys by Jeff Phillips Aug 96
c
	common /work/ z(2000), zs(2000), spare(26000)
	common /ogrid/ jf, iproj, cmerid, baslat
        dimension cmbl(2),ibie(2),idbe(3), colz(2),xycr(4)
        character*56 id, id2, if, jf, kf, blank
        character ians*1, ifunc*20, p*8, prjnam*80, prompt*80
        data dv/1.e38/,
     1    itty,idv,jdv,kdv,ldv/5,12,13,14,15/
	dval  = dv
	big   = 1.e29
	blank = ' '
c
        call askin
        call pfinit('utility')
c
        ifunc='?'
1       call askc('function',ifunc,ierr)
        if ( ierr .eq. -2 ) stop
        if(ifunc.eq.'help')ifunc='zilch'
c  convert ifunc to lower case
	call cvc( ifunc, 2 )
	ians = ifunc(1:1)
        ifunc='s'
        if ( ians .eq. 's' ) stop
        if(ians.eq.'p')then
          print*,'The plug option is disabled; use PLUGGRID'
          go to 1
        endif
20      if=' '
        if(ians.eq.'x' .or. ians.eq.'m' .or. ians.eq.'e' .or.
     1  ians.eq.'h' .or. ians.eq.'c' .or. ians.eq.'l' .or.
     1  ians.eq.'r') go to 30
        print  25
25	format(' available operations',/,
     1 ' r reorder elements',/,
     1 ' x extract subset',/,
     1 ' l look at grid values',/,
     1 ' h change header info',/,
     1 ' e edit grid points',/,
     1 ' m merge subset- xo,yo in number of cols,rows offset',/,
     1 ' c binary xyz coordinate output',/,
     1 ' s stop')
	go to 1
c
30      call askc('input grid',if,ierr)
        if(ierr.eq.-2) go to 1
        call gopen(idv,if,'old','read',ierr)
        if(ierr.ne.0) then
        print *, ' could not open grid, please reenter filename'
	  go to 30
	endif
        call gheader('r',idv,id,nc,nr,xo,dx,yo,dy,ierr)
c  get the old creation program name
        call gpgm(p)
50    print 55, id, p, nc, nr, xo, yo, dx, dy
55      format(' id=',a56, 'pgm=',a8/, ' size= ',2i5,
     1 /,' xo,yo =', 1p2g14.7,'  dx&dy=',2g14.5)
c  get the projection name
        call mpsnam('r',1,prjnam,nchp,ierr)
        if(ierr.eq.0) then
          print 56, prjnam(1:nchp)
56        format(' projection = ',a)
          if(prjnam.ne.'not specified') then
c  get the projection parameters
            call mpsref('r',1,cmbl,ncmbl,ierr)
            if(ierr.eq.0) then
              print *, 'central meridian & base latitude =', cmbl
            endif
          endif
        endif
c
	if ( ians .ne. 'l' ) then
c  open output grid
          jf=' '
39        if(ians.ne.'c') then
            if(ians.eq.'x') then
              jf=if(1:index(if,'.'))//'ext'
            else if(ians.eq.'e') then
              jf=if(1:index(if,'.'))//'edt'
            else if(ians.eq.'r') then
              jf=if(1:index(if,'.'))//'rot'
            else if(ians.eq.'m') then
              jf=if(1:index(if,'.'))//'mrg'
            else if(ians.eq.'h') then
              jf=if(1:index(if,'.'))//'hdr'
            endif
            call askc('output grid',jf,ierr)
            call gopen(jdv,jf,'new','write',ierr)
          else
            jf=if(1:index(if,'.'))//'xyz'
            call askc('output xyz file',jf,ierr)
            call xopen(jdv,jf,'new','write',itype,ierr)
          endif
          if(ierr.eq.-2) then
            close(idv)
            go to 30
          endif
          if(ierr.ne.0) then
            print*,'Error - try again'
            go to 39
          endif
        endif
	length = 26
	if ( nc + 1 .gt. length ) length = nc + 1
c
c       if(ians.eq.'p') go to 100
	if(ians.eq.'m') go to 200
	if(ians.eq.'e') go to 300
	if(ians.eq.'h') go to 400
	if(ians.eq.'c') go to 500
	if(ians.eq.'l') go to 600
	if(ians.eq.'x') go to 700
	if(ians.eq.'r') go to 800
	go to 1
c
c  plug holes using minimum curvature
c100     open( unit=jdv, file=jf, form='unformatted', status='new',
c     1       recl=length )
c        call utp( dval, itty, idv, jdv )
c        go to 999
c
c  merge function, secondary must be a subset of the primary grid
200     call askc('secondary filename',kf,ierr)
        if(ierr.eq.-2) go to 999
        call gopen(kdv,kf,'old','read',ierr)
        if(ierr.ne.0) go to 200
        call gheader('r',kdv,id2,nc2,nr2,xo2,dx2,yo2,dy2,ierr)
	if(dx2.ne.dx) stop ' dx1 not= dx2'
	if(dy2.ne.dy) stop ' dy1 not= dy2'
c
        if(ierr.ne.0) go to 999
	gx=(xo2-xo)
	gy=(yo2-yo)
	ix=1+int(gx/dx)
	iy=1+int(gy/dy)
	ix2=ix-1+nc2
	iy2=iy-1+nr2
	if(ix.le.0 .or. iy.le.0 .or. ix2.gt.nc .or. iy2.gt.nr) then
       print *,' secondary grid is not a subset of the primary grid'
       print *,' starting x,y index, ending x,y index'
       print *,ix,iy,ix2,iy2
	 go to 1
	endif
	gx2=amod(gx,dx)
	gy2=amod(gy,dy)
	if(gx2.gt.0.1 .and. gx2.lt.0.9) stop ' x direction not registered'
	if(gy2.gt.0.1 .and. gy2.lt.0.9) stop ' y direction not registered'
	ir=0
        call gheader('w',jdv,id,nc,nr,xo,dx,yo,dy,ierr)
	if(iy.eq.1) go to 211
	do 210 i=1,iy-1
        call grow('r',idv,i,z,nc,ierr)
	ir=ir+1
210     call grow('w',jdv,i,z,nc,ierr)
211	do 220 i=1,nr2
        call grow('r',idv,i,z,nc,ierr)
	ir=ir+1
        call grow('r',kdv,i,zs,nc2,ierr)
	do 230 j=1,nc2
	k=ix-1+j
230	z(k)=zs(j)
220     call grow('w',jdv,i,z,nc,ierr)
	if(ir.eq.nr) go to 299
	do 240 i=ir+1,nr
        call grow('r',idv,i,z,nc,ierr)
	ir=ir+1
240     call grow('w',jdv,i,z,nc,ierr)
299     call gclose(kdv,'keep')
	go to 999
c
c  edit function
300	ir=0
        call gheader('w',jdv,id,nc,nr,xo,dx,yo,dy,ierr)
        if(ierr.ne.0) then
          print*,'Error reading header'
          go to 999
        endif
        ibie(1)=1
        ibie(2)=min(10,nc)
        idbe(1)=0
        idbe(2)=0
        idbe(3)=0
        print*,'edit function'
        print*,' set first parameter =  0 or "mu" if done with question'
        print*,'                     = -1 to write row'
        jr=1
301     call aski4('row to edit',jr,ierr)
        if(ierr.eq.-2) jr=nr+1
        if(jr.le.0) jr=nr+1
303     ir=ir+1
        call grow('r',idv,ir,z,nc,ierr)
        if(ierr.ne.0) stop 'error reading row'
	if(jr.le.ir .and. jr.le.nr) go to 309
        call grow('w',jdv,ir,z,nc,ierr)
        if(ierr.ne.0) stop 'error writing row'
	if(ir.ge.nr) go to 999
	go to 303
309     write(prompt,340) ir
        iedit=0
340     format('preview row',i5,'?')
        iyn=0
        call aski4l(prompt,iyn,ierr)
        if(iyn.eq.0) go to 310
         ibie(1)=1
         ibie(2)=10
345      call aski4a('col: from to',ibie,2,ierr)
         if(ierr.eq.-2) then
           if(iedit.eq.0) go to 310
           jans=1
           call aski4l('write row?',jans,ierr)
           if(ierr.eq.-2) go to 399
           if(jans.eq.0) go to 310
           go to 399
         endif
        if(ibie(1).lt.0) go to 399
        if(ibie(1).eq.0) go to 310
      print  332,(z(i),i=ibie(1),ibie(2))
        inc=ibie(2)-ibie(1)
        ibie(1)=ibie(2)+1
        ibie(2)=ibie(1)+inc
        go to 345
310     idbe(1)=0
        idbe(2)=0
        idbe(3)=0
        iedit=1
        call aski4a('dval: from to number-of-rows',idbe,3,ierr)
        if(ierr.eq.-2) go to 320
        if(idbe(1).lt.0) go to 399
        if(idbe(1).eq.0) go to 320
        if(idbe(3).eq.0) idbe(3)=1
        idbe1sav=idbe(1)
313     do 312 i=idbe(1),idbe(2)
312     z(i)=dval
         if(idbe(3).lt.1) then
          print*,'editing row',jr
          go to 310
        endif
        call grow('w',jdv,ir,z,nc,ierr)
        if(ierr.ne.0) stop 'error writing row'
        call grow('r',idv,ir,z,nc,ierr)
        if(ierr.ne.0) stop 'error reading row'
	ir=ir+1
        jr=jr+1
        idbe(3)=idbe(3)-1
        go to 313
320     colz(1)=0.
        colz(2)=0.
        call askf4a('column_number new_value',colz,2,ierr)
        if(ierr.eq.-2) go to 330
        ndx=int(colz(1))
        if(ndx.eq.0) go to 330
        if(ndx.lt.0) go to 399
        z(ndx)=colz(2)
	go to 320
330     iyn=0
        call aski4l('print row?',iyn,ierr)
        ibie(1)=1
        ibie(2)=10
        if(iyn.eq.1) go to 345
        jans=1
        call aski4l('write row?',jans,ierr)
        if(ierr.eq.-2) go to 399
        if(jans.eq.0) go to 310
332     format(1x,1p5e12.4)
399        call grow('w',jdv,jr,z,nc,ierr)
        if(ierr.ne.0) stop 'error writing row'
        jr=jr+1
	go to 301
c
c  change header info
400	ios = 0
401     call askalt
        call askc('enter new title',id,ierr)
        call askalt
        if(ierr.eq.-2) go to 999
402     call askc('enter program id',p,ierr)
        if(ierr.eq.-2) go to 401
        call pfinit(p)
403     call askf4('enter new x origin',xo,ierr)
        if(ierr.eq.-2) go to 402
404    call askf4('enter new y origin',yo,ierr)
       if(ierr.eq.-2) go to 403
405     call askf4('enter new x increment',dx,ierr)
406     call askf4('enter new y increment',dy,ierr)
        call mphnam()
        call mpxnam(2,prjnam,iproj,ierr)
407     call aski4('enter new projection code',iproj,ierr)
        if(ierr.eq.-2) go to 406
        call mpxnam(1,prjnam,iproj,ierr)
        call mpsnam('w',1,prjnam,nchp,ierr)

408     call askf4('enter new central meridian',cmbl(1),ierr)
        if(ierr.eq.-2) go to 407
409     call askf4('enter new base latitude',cmbl(2),ierr)
        if(ierr.eq.-2) go to 408
        call mpsref('w',1,cmbl,ncmbl,ierr)
        if(ierr.ne.0) go to 999
        call gheader('w',jdv,id,nc,nr,xo,dx,yo,dy,ierr)
	do j = 1, nr
          call grow('r',idv,i,z,nc,ierr)
          call grow('w',jdv,i,z,nc,ierr)
	enddo
	go to 999
c
c  coordinate output
500	nxb=1
	nxe=nc
	nyb=1
	nye=nr
c        call xopen(jdv,jf,'new','write',itype,ierr)
c        if(ierr.ne.0) go to 999
        iyn=1
410     call aski4l( 'output whole grid ?', iyn, ierr )
        if(ierr.eq.-2) go to 999
        if(iyn.eq.0) then
          xycr(1)=xo
          xycr(2)=yo
          xycr(3)=nc
          xycr(4)=nr
          call askf4a('xo,yo (data units), nc,nr',xycr,4,ierr)
          if(ierr.eq.-2) go to 410
          xo2=xycr(1)
          yo2=xycr(2)
          nc2=xycr(3)
          nr2=xycr(4)
	  xo2=xo2 + dx*0.499
	  nxb=int((xo2-xo)/dx)+1
	  nxe=nxb+nc2-1
	  if(nxb.lt.1)  nxb=1
	  if(nxe.gt.nc) nxe=nc
	  yo2=yo2 + dy*0.499
	  nyb=int((yo2-yo)/dy)+1
	  nye=nyb+nr2-1
	  if(nyb.lt.1) nyb=1
	  if(nye.gt.nr) nye=nr
	endif
	if((nye-nyb).le.0 .or. (nxe-nxb).le.0) then
        print *,' specified nc or nr <1'
	  go to 500
	endif
c
	xo3=xo + dx*(nxb-1)
	yo3=yo + dy*(nyb-1)
	if(nyb.gt.1) then
	  do j=1,nyb-1
	  read(idv)
	  enddo
	endif
c
	y=yo3
	do 550 j=nyb,nye
        call grow('r',idv,j,z,nc,ierr)
        if(ierr.ne.0) then
        print  *,' eof while reading'
	  go to 999
	endif
	x=xo3
	do 540 i=nxb,nxe
	if(z(i).gt.1.e30) go to 540
        write(jdv) x,y,z(i)
540	x=x+dx
550	y=y+dy
	go to 999
c
c  look mode
600	ir=0
        jr=0
610     jr=jr+1
        call aski4('row:',jr,ierr)
        if(ierr.eq.-2) go to 699
	if(jr.lt.ir) go to 699
        if(jr.eq.0) go to 610
	if(jr.eq.ir) go to 605
	if(jr.eq.ir+1) go to 603
	do 602 i=ir+1,jr-1
	read(idv,end=699)
602	ir=ir+1
603     call grow('r',idv,jr,z,nc,ierr)
        if(ierr.ne.0) go to 699
	ir=ir+1
         ibie(1)=1
         ibie(2)=10
605      call aski4a('col: from to',ibie,2,ierr)
         if(ierr.eq.-2) go to 610
        if(ibie(1).lt.0) go to 610
        if(ibie(1).eq.0) go to 699
        print  332,(z(i),i=ibie(1),ibie(2))
        inc=ibie(2)-ibie(1)
        ibie(1)=ibie(2)+1
        ibie(2)=ibie(1)+inc
	go to 605
699	close(idv)
	go to 1
c
700	maxcol = 4000
	call utx( dval, maxcol, z, itty, idv, jdv )
	go to 999
c
800	call utr( itty, idv, jdv, kdv, ldv )
c
999	close(idv)
	close(jdv)
	go to 1
	end
c***********************************************************************
	subroutine utx( dval, maxcol, z, itty, idv, jdv )
c  grid utility extract function
c  no restriction on output grid size or location (subject to maxcol)
	common /ogrid/ ofile, iproj, cmerid, baslat
        dimension      z(maxcol)
        dimension ixiy(2), xoyo(2), ncnr(2), ixye(2), xeye(2)
        character      ofile*56, id*56, prompt*80
	irow=0
	jrow=0
        call gheader('r',idv,id,nc,nr,xo,dx,yo,dy,ierr)
	go to 5
1     print *,' conversion error...'
5      call askalt
400    call askc('output title',id,ierr)
       call askalt
c
      print *,' >> new origin'
15    print *,' enter either 0 for integer starting location'
c      print *,'           or 1 for lower-left coordinate'
      prompt='           or 1 for lower-left coordinate'
      iswt=0
401   call aski4(prompt,iswt,ierr)
      if(ierr.eq.-2) go to 400
      if(iswt.lt.0 .or. iswt.gt.1) go to 15
      if(iswt.eq.0) then
        print *,' enter starting column and row number'
        print *,' 1,1 is the the start of the parent grid'
c        print *,' negative numbers ok'
        prompt=' negative numbers ok'
        ixiy(1)=1
        ixiy(2)=1
        call aski4a(prompt,ixiy,2,ierr)
        if(ierr.eq.-2) go to 15
      else
c        print *,' enter lower-left x,y coordinate in data units'
        prompt=' enter lower-left x,y coordinate in data units'
        xoyo(1)=xo
        xoyo(2)=yo
        call askf4a(prompt,xoyo,2,ierr)
        if(ierr.eq.-2) go to 15
c        print*,xoyo
        gx=(xoyo(1)-xo)/dx
        gy=(xoyo(2)-yo)/dy
        print*,gx,gy
        ixiy(1)=1+int(gx+sign(0.5,gx))
        ixiy(2)=1+int(gy+sign(0.5,gy))
c        print*,ixiy
      endif
c      xo2=xo+dx*float(ixiy(1)-1)
c      yo2=yo+dy*float(ixiy(2)-1)
      xo2=(ixiy(1)-1)
      yo2=(ixiy(2)-1)
c      print*,xo2,yo2
      xo2=xo2*dx
      yo2=yo2*dy
c      print*,xo2,yo2
      xo2=xo+xo2
      yo2=yo+yo2
c      print*,xo2,yo2,xo,yo,dx,dy
c
      print *,' >> output size'
20    print *,' enter either 0 for ending column and row number'
      print *,'           or 1 for upper-right coordinate'
c      print *,'           or 2 for specified ncol, nrow'
      prompt='           or 2 for specified ncol, nrow'
      iswt=2
      call aski4(prompt,iswt,ierr)
      if(ierr.eq.-2) go to 15
      if(iswt.lt.0 .or. iswt.gt.2) go to 20
      if(iswt.eq.2) then
        ncnr(1)=nc
        ncnr(2)=nr
c        print *,' enter new dimensions (integer ncol,nrow)'
        prompt=' enter new dimensions (integer ncol,nrow)'
        call aski4a(prompt,ncnr,2,ierr)
        if(ierr.eq.-2) go to 20
      else if(iswt.eq.1) then
c        print *,' enter upper-right x,y coordinate in data units'
        prompt=' enter upper-right x,y coordinate in data units'
        xeye(1)=xo+float(nc-1)*dx
        xeye(2)=yo+float(nr-1)*dy
        call askf4a(prompt,xeye,2,ierr)
        if(ierr.eq.-2) go to 20
        ncnr(1)=1+int(0.5+(xeye(1)-xo2)/dx)
        ncnr(2)=1+int(0.5+(xeye(2)-yo2)/dy)
      else
c        print *,' enter ending column and row number in the parent grid'
        prompt=' enter ending column and row number in the parent grid'
        ixye(1)=nc
        ixye(2)=nr
        call aski4a(prompt,ixye,2,ierr)
        if(ierr.eq.-2) go to 20
        xe2=xo+dx*float(ixye(1)-1)
        ye2=yo+dy*float(ixye(2)-1)
        ncnr(1)=1+int(0.5+(xe2-xo2)/dx)
        ncnr(2)=1+int(0.5+(ye2-yo2)/dy)
      endif
      if(ncnr(1).le.0 .or. ncnr(2).le.0) go to 20
c
      print *,' output grid ncol,nrow will be ', ncnr
      length = 26
      if ( ncnr(1) + 1 .gt. length ) length = ncnr(1) + 1
c  ix1,ix2 are starting addresses for the input and output rows
      if(ixiy(1).gt.0) then
        ix1=1
        ix2=ixiy(1)
      else
        ix1=2-ixiy(1)
        ix2=1
      endif
      ic =ix1+nc-1
      ic2=ix2+ncnr(1)-1
      if(ic2.gt.ic) ic=ic2
      if(ic.gt.maxcol) then
        print *,maxcol,' column temp storage exceeded'
        return
      else
        do 50 i=1,ic
50      z(i)=dval
      endif
      call gheader('w',jdv,id,ncnr(1),ncnr(2),xo2,dx,yo2,dy,ierr)
c
      if ( ixiy(2) .eq. 1 ) go to 90
c  starting condition
      if(ixiy(2).le.0) then
        do 70 i=ixiy(2),0
        call grow('w',jdv,i,z(ix2),ncnr(1),ierr)
        jrow=jrow+1
        if(jrow.ge.ncnr(2)) go to 200
70      continue
      else
        do 80 i=1,ixiy(2)-1
        read(idv,end=120)
        irow=irow+1
        if(irow.ge.nr) go to 120
80      continue
      endif
c
c  transfer grid1 to grid2
90    input=1
100   if(jrow.ge.ncnr(2)) go to 200
      go to (110,150) input
110   call grow('r',idv,jrow,z(ix1),nc,ierr)
      irow=irow+1
120   if(irow.gt.nr .or. ierr.ne.0) then
        input=2
        do 130 i=1,ic
130     z(i)=dval
      endif
150   call grow('w',jdv,irow,z(ix2),ncnr(1),ierr)
      jrow=jrow+1
      if(input.eq.2) go to 100
      y1=yo+ dy*float(irow-1)
      y2=yo2+dy*float(jrow-1)
      ychk=2.0*abs(y2-y1)
      if(ychk.gt.dy .and. irow.lt.nr)
     1 print *,' utx: y error irow,jrow,yi,yj =',irow,jrow,y1,y2
      go to 100
200   return
      end
c***********************************************************************
	subroutine utr( itty, isw, jsw, ksw, lsw )
c  grid reordering function
	common /work/  buf(30000)
	common /ogrid/ ofile, iproj, cmerid, baslat
c        character ofile*56, p*8, p2*8, id*56, id2*56, blank*56
        character ofile*56, p2*8, id*56, id2*56, blank*56
	data nbuf/ 30000 /
	blank = ' '
c
        call gheader('r',isw,id,nc,nr,xo,dx,yo,dy,ierr)
	x2  = float( nc - 1 ) * dx + xo
	y2  = float( nr - 1 ) * dy + yo
	dx2 = -dx
	dy2 = -dy
c
      id2=id
      call askalt
      call askc('output title',id2,ierr)
      call askalt
      print  15
15	format(' enter 1 for transpose',/,
     1      '       2 for column reverse',/,
     1      '       3 for row    reverse',/,
     1      '       4 for 90  deg counter-clockwise rotation',/,
     1      '       5 for 90  degree      clockwise rotation ',/,
     1      '       6 for 180 degree rotation')
1	iop = 0
	read( itty, *, err=2, end=2 ) iop
2	if ( iop .lt. 1 .or. iop .gt. 6 ) then
        print *, iop, ' not recognized'
	  go to 1
	endif
c
	length = 26
	if ( iop.eq.1 .or. iop.eq.4 .or. iop.eq.5 ) then
c  flip col and rows, record length = nr + 1
	  if ( nr + 1 .gt. length ) length = nr + 1
	  else
	  if ( nc + 1 .gt. length ) length = nc + 1
	endif
        call gopen(jsw,ofile,'new','write',ierr)
        if(ierr.ne.0) stop 'error opening file'
c	
	if ( iop .eq. 1 ) then
c  transpose
	p2 = 'transpos'
        call pfinit('transpos')
        call gheader('w',jsw,id2,nr,nc,yo,dy,xo,dx,ierr)
	call xpose( nc, nr, buf, nbuf, isw, jsw, ksw )
	go to 999
	endif
c
	if ( iop .eq. 2 ) then
c  reverse columns
	p2 = 'rev_col'
        call gheader('w',jsw,id2,nc,nr,x2,dx2,yo,dy,ierr)
	call revcol( nc, nr, buf, isw, jsw )
	go to 999
	endif
c
c  reverse rows
	if ( iop .eq. 3 ) then
        call pfinit('rev_row')
        call gheader('w',jsw,id2,nc,nr,xo,dx,y2,dy2,ierr)
	call revrow( nc, nr, buf, isw, jsw, ksw )
	go to 999
	endif
c
c  rotate 90 degree counter-clockwise
	if ( iop .eq. 4 ) then
      open( unit=lsw, status='scratch',
     1     form='unformatted' )
        call gheader('w',lsw,id2,nr,nc,yo,dy,xo,dx,ierr)
	call xpose( nc, nr, buf, nbuf, isw, lsw, ksw )
        call pfinit('-90clock')
        call gheader('w',jsw,id2,nr,nc,y2,dy2,xo,dx,ierr)

	call revcol( nr, nc, buf, lsw, jsw )
	close( lsw )
	go to 999
	endif
c
c  rotate 90 degree clockwise
	if ( iop .eq. 5 ) then
      open( unit=lsw, status='scratch',
     1     form='unformatted' )
        call pfinit('rev-col')
        call gheader('w',lsw,id2,nc,nr,x2,dx2,yo,dy,ierr)
	call revcol( nc, nr, buf, isw, lsw)
        call pfinit('+90clock')
        call gheader('w',jsw,id2,nr,nc,yo,dy,x2,dx2,ierr)
	call xpose( nc, nr, buf, nbuf, lsw, jsw, ksw )
	close( lsw )
	go to 999
	endif
c
c  rotate 180 degree
	if ( iop .eq. 6 ) then
      open( unit=lsw, status='scratch',
     1     form='unformatted' )
        call pfinit('rev-col')
        call gheader('w',lsw,id2,nc,nr,x2,dx2,yo,dy,ierr)
	call revcol( nc, nr, buf, isw, lsw)
        call pfinit('180clock')
        call gheader('w',jsw,id2,nc,nr,x2,dx2,y2,dy2,ierr)
	call revrow( nc, nr, buf, lsw, jsw, ksw )
	close( lsw )
	endif
c
999	return
	end
c***********************************************************************
	subroutine xpose( nc, nr, buf, nbuf, isw, jsw, ksw )
c  transpose using a random access file
	dimension buf(nbuf)
	if ( nc .gt. nbuf / 20 ) stop ' 20 * ncol > nbuf'
	if ( nr .gt. nbuf / 20 ) stop ' 20 * nrow > nbuf'
      open( ksw, status='scratch',
     1     form='unformatted', access='direct', recl=1600 )
	m   = 1 + ( nc - 1 ) / 20
	n   = 1 + ( nr - 1 ) / 20
	n20 = 20 * n
	m20 = 20 * m
	call rio( buf, m, n, 1, 1, m20, -1, isw, ksw )
	ic = 1
c
	do 100 ibuf = 1, m
	loc = ibuf
	call rio( buf, n, 0, loc, m, n20, 1, isw, ksw )
c
	ndex = 1
	do irow = 1, 20
c          call rowio (nr, buf(ndex), 0, jsw, jsw, iend )
          call grow('w',jsw,irow,buf(ndex),nr,ierr)
	  if ( ic .eq. nc ) go to 999
	  ic   = ic   + 1
	  ndex = ndex + n20
	enddo
100	continue
c
999	close(ksw)
	return
	end
c***********************************************************************
	subroutine revcol( nc, nr, z, isw, jsw )
	dimension ihdr(16), z(nc)
	rewind(isw)
	read(isw) ihdr, nc1, nr1
	if ( nc1 .ne. nc  .or.  nr1 .ne. nr ) then
        print *, ' passed dimensions (ncol, nrow):', nc, nr
        print *, ' header dimensions             :', nc1, nr1
	  stop
	endif
	do j = 1, nr
	  read( isw, end=10, iostat=ios ) xx, z
10	  if ( ios .ne. 0 ) then
          print *, ' premature EOF at row', j
	    return
	  endif
	  i2 = nc
	  do i = 1, nc / 2
	    tmp   = z(i2)
	    z(i2) = z(i)
	    z(i)  = tmp
	    i2    = i2 - 1
	  enddo
	  write( jsw ) xx, z
	enddo
	return
	end
c***********************************************************************
	subroutine revrow( nc, nr, z, isw, jsw, ksw )
	dimension z(nc)
      open( ksw, status='scratch',
     1     form='unformatted', access='direct', recl=4*(nc+1) )
	rewind(isw)
	read(isw)
	do j = 1, nr
	  read(isw) xx, z
	  write(ksw,rec=j) xx, z
	enddo
	j2 = nr
	do j = 1, nr
	  read(ksw,rec=j2) xx, z
	  write(jsw) xx, z
	  j2 = j2 - 1
	enddo
	close(ksw)
	return
	end
c***********************************************************************
	subroutine rio(z,nblk,nbuf,loc2,idloc,lenz,iop,isw,ksw)
c iop=-1 read grid, pad with dval's, and output with implied iop=1
c    =0 output 20x20 tranposed blocks
c    =1 input  20x20 block
c buffer z is dimensioned lenz x 20
c random access is nblk x nbuf, where 'loc' increments by 'idl'
c        dimension hdr(16),w(400),z(1)
        dimension w(400),z(1)
        character id*56
	logical more
        data dval/0.1e+39/
	loc=loc2
	assign 999 to iret
c
	if(iop)1,2,3
1	assign 99 to iret
	more=.true.
	itot=0
        rewind(isw)
        call gheader('r',isw,id,nc,nr,xo,dx,yo,dy,ierr)
	do 100 ibuf=1,nbuf
	do i=1,lenz*20
	  z(i)=dval
	enddo
	ireadz=1
	do irow=1,20
	  if(more) then
             call grow('r',isw,irow,z(ireadz),nc,ierr)
	    ireadz=ireadz+lenz
            if(ierr.eq.0) itot=itot+1
            if(ierr.ne.0 .or. itot.ge.nr) more=.false.
          endif
	enddo
	  go to 2
99	  continue
100	continue
      if(itot.ne.nr) print *,' found only ',itot,' of ',nr,' rows'
	return
c
c  extract 20x20 block, tranpose, and output
2	do iblk=1,nblk
	iws=1
	izs=(iblk-1)*20+1
	do j=1,20
	  iw=iws
	  iz=izs
	  do i=1,20
	    w(iw)=z(iz)
	    iw=iw+20
	    iz=iz+1
	  enddo
	  iws=iws+1
	  izs=izs+lenz
	enddo
      write(ksw,rec=loc) w
	loc=loc+idloc
	enddo
	go to iret,(99,999)
c
c  insert 20x20 block into buffer
3	do iblk=1,nblk
      read(ksw,rec=loc) w
	iw=1
	izs=(iblk-1)*20+1
	do j=1,20
	  iz=izs
	  do i=1,20
	    z(iz)=w(iw)
	    iw=iw+1
	    iz=iz+1
	  enddo
	izs=izs+lenz
	enddo
	loc=loc+idloc
	enddo
c
999	return
	end
c***********************************************************************
	subroutine cvc(a,i2)
c  i2=1 for lower case to upper
	character a*(*)
	n=len(a)
	if(i2.lt.1 .or. i2.gt.2) stop ' cvc: invalid convert parameter'
	do 40 i=1,n
	j=ichar(a(i:i))
	go to (10,20) i2
10	if(j.gt.122 .or. j.lt.97) go to 40
	j=j-32
	go to 30
20	if(j.gt.90 .or. j.lt.65) go to 40
	j=j+32
30	a(i:i)=char(j)
40	continue
	return
	end

