c       TMERGER
c	Merges two overlapping data sets gridded using the same origin
c	and grid parameters.  The data sets are weighted and summed in
c	the area of overlap.
c
c       written by Jeff Phillips and Tien Grauch
c
c      basic code from program JMRG by Phillips.
c      modified for oddf by Tien Grauch March 1994
c      Grauch added code to construct the large-size matching grids before 
c       merging and plug only in the overlap region Oct 1994
c      modified to restrict overlap areas by Tien Grauch Nov 1994
c
c      parameter nwork is the max rows X cols allowed for overlap region
c
        parameter ( nwork=50000, maxcol=1000 )
        dimension iwork(nwork), work(nwork)
c        equivalence(iwork,work)
        dimension d1(maxcol),d2(maxcol)
      common /specs/xo,yo,nc,nr,dx,dy,xo2,yo2,nc2,nr2
      common /rcspecs/ics1,ics2,ice1,ice2,jrs1,jrs2,jre1,jre2,
     1 icso,iceo,jrso,jreo
	character*50 ifil1,ifil2,ofile
cu      character state*6, ans*1, id*80,ido*80
      character state*6, ans*1, id*56,ido*56
      real*8 dsum
        data dval/1.0e+38/
        data dvtest/1.0e+36/
c state is status of temp files upon closing.  change to 'keep'
c  for debugging purposes
c     state='keep'
      state='delete'
      call pfinit( 'tmerger')
      print 800
800   format(/' Program gdmerger, version 1.0'//
     1' TMERGER merges two overlapping data sets by weighting'/
     2' then summing in the area of overlap.'/
     3' To accommodate the temporary files generated',
     4' during this program,'/' you will need 3 times the size of',
     5' the output grid + 2 times the'/' size of the overlap between',
     6' grids.'//' Technical contacts - Jeff Phillips and Tien Grauch'/)
1      write(6,100)
100	format(' first grid: '$)
	read(5,105) ifil1
105	format(a)
        call gopen(10, ifil1, 'old', 'read', ierr)
        if(ierr.ne.0) go to 1
cu      call gh1f4('r', 10, id, nc,nr,xo,dx,yo,dy,ierr)
        call gheader('r',10,id,nc,nr,xo,dx,yo,dy,ierr)
      print 812, id,xo,yo,nc,nr,dx,dy
cu812   format(/'title: ',a80/5x,'xo=',g10.4,15x,'yo=',g10.4,/
812   format(/' title: ',a56/5x,'xo=',g10.4,15x,'yo=',g10.4,/
     1 5x,'nc=',i4,21x,'nr=',i4/5x,'dx=',g10.4,15x,'dy=',g10.4/)
2      write(6,110)
110	format(' second grid: '$)
	read(5,105) ifil2
        call gopen(11, ifil2, 'old', 'read', ierr)
        if(ierr.ne.0) go to 2
cu      call gh1f4('r', 11, id, nc2,nr2,xo2,dx2,yo2,dy2,ierr)
        call gheader('r',11,id,nc2,nr2,xo2,dx2,yo2,dy2,ierr)
      print 812, id,xo2,yo2,nc2,nr2,dx2,dy2
c
4       print 801
801   format('Merge options:'/5x,'1 - treat both grids with equal',
     1' priority'/5x'2 - restrict overlap region, priority to grid 1',
     2/5x,'3 - restrict overlap region, priority to grid 2'/
     3/10x,'enter option: ',$)
      read*,iopt
      if(iopt.lt.1.or.iopt.gt.3) go to 4
      nband=0
      if(iopt.gt.1) then
        print 802
802     format(' Restrict overlap region to how many grid pts? ',$)
        read*,nband
        if(nband.le.0) nband=1
      endif
c
	write(6,130)
130	format(' output grid: '$)
	read(5,105) ofile
      print*,'Enter output title'
      read 106, ido
cu106   format(a80)
106   format(a)
c
c check grid nodes  -- this section inserted by Tien from gdbiharm2.f
      if(dx.ne.dx2.or.dy.ne.dy2) stop 'grid intervals unequal'
      print*
cu      type*,'xo,xo2,dx=',xo,xo2,dx
cu      type*,'yo,yo2,dy=',yo,yo2,dy
      amodx=abs(amode(xo-xo2,dx))
      amody=abs(amode(yo-yo2,dy))
      print*,'xo,xo2,dx,amodx=',xo,xo2,dx,amodx
      print*,'yo,yo2,dy,amody=',yo,yo2,dy,amody
cu      if(amod((xo-xo2),dx)).gt.0.00001) then
      if(amodx.gt.0.0001) then
        print 814
814   format('grid nodes do not match (roundoff?).'/
     1' Continue anyway? ',$)
        read 105, ans
        if(ans.ne.'y'.and.ans.ne.'Y') stop
      else
cu         if(amod((yo-yo2),dy).gt.0.00001) then
         if (amody.gt.0.0001) then
          print 814
          read 105, ans
          if(ans.ne.'y'.and.ans.ne.'Y') stop
         endif
      endif
c
c find specs of output and specs of overlap work area 
c
      call find_specs(xst,yst,nclap,nrlap,xout,yout,ncout,nrout,
     1 iEW,iNS)
      if(ncout.gt.maxcol) then
         print 815, maxcol
815      format('too many output cols, max is ',i5)
         stop
      endif
      if(nclap*nrlap.gt.nwork) then
        print 816, nwork
816     format(' merge-boundary work area too big, max nc*nr =',i8)
c        stop
      endif
      if(nband.ge.nclap. or.nband.ge.nrlap) then
        print*,'restriction on overlap region too big'
        stop
      endif
c
c  output grid1 & grid2 data into big grids , units 14,15, resp.
c  non-defined grid nodes are set to dvals
c
cu      call gopen(14, 'jmrg1.tmp', 'new', ' ', ierr)
cu      call gopen(15, 'jmrg2.tmp', 'new', ' ', ierr)
cu      call gh1f4('w', 14, id,ncout,nrout,xout,dx,yout,dy,ierr)
cu      call gh1f4('w', 15, id,ncout,nrout,xout,dx,yout,dy,ierr)
      call gopen(14, ' ', 'scratch', 'readwrite', ierr)
      call gopen(15, ' ', 'scratch', 'readwrite', ierr)
      call gheader('w',14,id,ncout,nrout,xout,dx,yout,dy,ierr)
      call gheader('w',15,id,ncout,nrout,xout,dx,yout,dy,ierr)
c construct large grid1
      jr1=0
      do 220 j=1,nrout

      do 211 k=1,maxcol
      d2(k)=dval
211   continue

      jrow=j
      if(jrow.ge.jrs1.and.jrow.le.jre1) then
         jr1=jr1+1
cu         call growf4('r', 10, jr1, d1, nc, ierr)
         call grow('r', 10, jr1, d1, nc, ierr)
         ic1=0
         do 215 i=1,ncout 
         icol=i
         if(icol.ge.ics1.and.icol.le.ice1) then
            ic1=ic1+1
            d2(icol)=d1(ic1)
         endif
215      continue
      endif

cu      call growf4('w', 14, jrow, d2, ncout, ierr)
      call grow('w', 14, jrow, d2, ncout, ierr)

220   continue
c
c construct large grid2
      jr2=0
      do 320 j=1,nrout

      do 311 k=1,maxcol
      d1(k)=dval
311   continue

      jrow=j
      if(jrow.ge.jrs2.and.jrow.le.jre2) then
         jr2=jr2+1
cu         call growf4('r', 11, jr2, d2, nc2, ierr)
         call grow('r', 11, jr2, d2, nc2, ierr)
         ic2=0
         do 315 i=1,ncout 
         icol=i
         if(icol.ge.ics2.and.icol.le.ice2) then
            ic2=ic2+1
            d1(icol)=d2(ic2)
         endif
315      continue
      endif

cu      call growf4('w', 15, jrow, d1, ncout, ierr)
      call grow('w', 15, jrow, d1, ncout, ierr)

320   continue
      call gclose(10, 'keep')
      call gclose(11, 'keep')
cu      call gclose(14, 'keep')
cu      call gclose(15, 'keep')
cu      call gopen(14, 'jmrg1.tmp', 'old', ' ', ierr)
cu      call gopen(15, 'jmrg2.tmp', 'old', ' ', ierr)
cu      call gh1f4('r', 14, id,ncout,nrout,xout,dx,yout,dy,ierr)
cu      call gh1f4('r', 15, id,ncout,nrout,xout,dx,yout,dy,ierr)
      call gheader('r',14,id,ncout,nrout,xout,dx,yout,dy,ierr)
      call gheader('r',15,id,ncout,nrout,xout,dx,yout,dy,ierr)
c following added for pc version
      if(jrso.eq.1) go to 330
      do 325 j=1,jrso-1
      call grow('r', 14, jrow, d1, ncout, ierr)
      call grow('r', 15, jrow, d2, ncout, ierr)
  325 continue
  330 continue

c  construct weighting grid for overlap area
c
cu      call gopen(12, 'jmrg.tmp', 'new', ' ', ierr)
      call gopen(12, ' ', 'scratch', 'readwrite', ierr)
cu      call gh1f4('w', 12, id,nclap,nrlap,xst,dx,yst,dy,ierr)
      call gheader('w',12,id,nclap,nrlap,xst,dx,yst,dy,ierr)
        print *,'...assigning weights of 0, 1, and unknown...'
        print*

c
c      weight of 1.0  - data available in grid1 only
c                0.0  - data available in grid2 only
c                dval - data available in both grids
c                0.5  - data available in neither grid     
c
c
c ngrid1 and ngrid2 are temporary flags to tell if grid1, grid2
c  have valid data in them or not.  ialldvals is a flag that
c  checks to see if there is any valid data
c    mncol,mnrow and mxcol,mxrow are values that will be used
c    to restrict the search when readjusting weights (if applicable)
c    they describe approx. where the area of dvals (overlap) are.
c
      mncol=999999
      mnrow=999999
      mxcol=-999999
      mxrow=-999999
      ialldvals=1
      jlap=0
      ndvals=0
      dsum=0.d0
      do 400 j=jrso, jreo
      jrow=j
cu      call growf4('r', 14, jrow, d1, ncout, ierr)
cu      call growf4('r', 15, jrow, d2, ncout, ierr)
      call grow('r', 14, jrow, d1, ncout, ierr)
      call grow('r', 15, jrow, d2, ncout, ierr)
      jlap=jlap+1

      ilap=0
      nrowpt= (jlap-1) * nclap

      do 350 i=icso, iceo
      ilap=ilap+1
      npt= nrowpt + ilap
      ngrid1=0 
        if(d1(i).lt.dvtest) ngrid1=1 
c
      ngrid2=0
        if(d2(i).lt.dvtest) ngrid2=1
c
c note:
c    ngrid1   ngrid2   wt 
c       0       0      .5
c       1       0       1
c       1       1      dval
c       0       1       0  
c
c assemble the work array with all the weights 
      if(ngrid1.eq.0) then
        if(ngrid2.eq.0) then
           work(npt)=0.5e0
        else
           work(npt)=0.0e0
           ialldvals=0
        endif
      else
        if(ngrid2.eq.0) then
           work(npt)=1.0e0
           ialldvals=0
        else
           work(npt)=dval
           ialldvals=0
           mncol=min(mncol,ilap)
           mxcol=max(mxcol,ilap)
           mnrow=min(mnrow,jlap)
           mxrow=max(mxrow,jlap)
           ndvals=ndvals+1
           dsum=dsum+d1(i)-d2(i)
        endif
      endif 
c
c
350   continue
cu      call growf4('w', 12, jlap, work(nrowpt+1), nclap, ierr)
      call grow('w', 12, jlap, work(nrowpt+1), nclap, ierr)
400   continue
c
      if(ialldvals.ne.0) stop 'no data in overlap area'
      diff=dsum/dfloat(ndvals)
      print*,'Average difference (grid1 - grid2) in overlap area = ',
     1diff
      print*,'Do you want to add this or some other non-zero value to gr
     1id1?'
      read 105, ans
      if(ans.eq.'y'.or.ans.eq.'Y') then
        print*,'Enter the value [',diff,']:'
        read(*,'(a)') ifil1
        leng=len_trim(ifil1)
        if(leng.ne.0) then
          lenm=leng-index(ifil1(1:leng),'.')
          if(lenm.eq.leng) lenm=0
          write(ifil2,'(a,i3,a,i3,a)') '(f',leng,'.',lenm,')'
          read(ifil1(1:leng),ifil2) diff
        endif
        print*,'Adding ',diff,' to grid1'
      else
        diff=0.d0
      endif
      mncol=mncol-1
      mxcol=mxcol+1
      mnrow=mnrow-1
      mxrow=mxrow+1
      if(mncol.lt.1.or.mncol.gt.nclap) mncol=1
      if(mnrow.lt.1.or.mnrow.gt.nrlap) mnrow=1
      if(mxcol.gt.nclap.or.mxcol.lt.1) mxcol=nclap
      if(mxrow.gt.nrlap.or.mxrow.lt.1) mxrow=nrlap
c
c if needed, readjust weights to reflect overlap restriction
c
      if(iopt.eq.1) go to 440
      print*
      print*,'readjusting weights to restrict overlap region...'
      print*
      if(iopt.eq.2) then
        ghipr=1.0e0
        glopr=0.0e0
      else
        ghipr=0.0e0
        glopr=1.0e0
      endif

c            
c find a value indicating data only in the low priority grid 
c  and replace all dvals within an nband radius with a flag 
c  value of 2.0
c
      ireadjust=0
      do 420 jlap=mnrow,mxrow
      nrowpt=(jlap-1)*nclap

      do 415 ilap=mncol,mxcol
      npt= nrowpt + ilap

      if(work(npt).eq.glopr) then
             do k=1, nband
             if(ilap+k.le.mxcol) then 
                nptk=npt + k
                if(work(nptk).eq.dval) work(nptk)=2.0e0
             endif
             if(ilap-k.ge.mncol) then 
                nptk=npt - k
                if(work(nptk).eq.dval) work(nptk)=2.0e0
             endif
             if(jlap+k.le.mxrow) then 
                nptk= npt + k*nclap
                if(work(nptk).eq.dval) work(nptk)=2.0e0
             endif
             if(jlap-k.ge.mnrow) then 
                nptk= npt - k*nclap
                if(work(nptk).eq.dval) work(nptk)=2.0e0
             endif
             end do
      endif
415   continue 
420   continue
c
c  replace all values of 2 with dvals and all dvals with the value
c   of the high priority grid and write out new weights
c
      do 430 jlap=1,nrlap
      nrowpt= (jlap-1)*nclap 
      do 425 ilap=1,nclap
      npt= nrowpt + ilap
      if(work(npt).eq.dval) then
         work(npt)=ghipr
         ireadjust=1
         go to 425
      endif
      if(work(npt).eq.2.0) then
         work(npt)=dval
         ireadjust=1
      endif
425   continue
cu      call growf4('w', 12, jlap, work(nrowpt+1), nclap, ierr)
      call grow('w', 12, jlap, work(nrowpt+1), nclap, ierr)
430   continue
      if(ireadjust.eq.0) print*,' no readjustment made'
c
c reopen unit 12 in order to read
c new unit 13 will be jmrg.wgt = plugged weighting grid, ready to use in
c       overlap area.
c
cu440   call gclose(12, 'keep')
cu      call gopen(12, 'jmrg.tmp', 'old', ' ', ierr)
cu      call gh1f4('r', 12, id, nclap, nrlap, xst, dx, yst, dy,ierr)
440   call gheader('r',12,id,nclap,nrlap,xst,dx,yst,dy,ierr)
cu      call gopen(13, 'jmrg.wgt', 'new', 'w', ierr)
      call gopen(13,'jmrg.wgt','new','readwrite',ierr)
cu      call gh1f4('w', 13, id,nclap,nrlap,xst,dx,yst,dy,ierr)
      call gheader('w',13,id,nclap,nrlap,xst,dx,yst,dy,ierr)
c
        nwrk=nclap*nrlap
        print*,ier,nwrk,nclap,nrlap
          call mcplug(20,12,13,ier,work,iwork,nwrk,nclap,nrlap)
c
	if(ier.eq.1) go to 91
c
      call gclose(12,state)
cu      call gclose(13,'keep')
cu      call gopen(13, 'jmrg.wgt', 'old', ' ', ierr)
cu      call gh1f4('r', 13, id,nclap,nrlap,xst,dx,yst,dy,ierr)
      call gheader('r',13,id,nclap,nrlap,xst,dx,yst,dy,ierr)
c unit 18 is final output grid
c note units 14 & 15 are still open
cu      call gopen(18, ofile, 'new', ' ', ierr)
      call gopen(18, ofile, 'new', 'write', ierr)
cu      call gh1f4('w', 18,ido,ncout,nrout,xout,dx,yout,dy,ierr)
      call gheader('w',18,ido,ncout,nrout,xout,dx,yout,dy,ierr)
      print 812, ido,xout,yout,ncout,nrout,dx,dy
c
c start making the final output grid.
c reuse the array work for output
c
c following added for pc version
c      rewind(14)
c      read(14) id
c      rewind(15)
c      read(15) id
      call gheader('r',14,id,ncout,nrout,xout,dx,yout,dy,ierr)
      call gheader('r',15,id,ncout,nrout,xout,dx,yout,dy,ierr)
      if(jrso.eq.1) go to 455
      do 450 j= 1,jrso-1
      jrow=j
cu      call growf4('r', 14, jrow, d1, ncout, ierr)
cu      call growf4('r', 15, jrow, d2, ncout, ierr)
      call grow('r', 14, jrow, d1, ncout, ierr)
      call grow('r', 15, jrow, d2, ncout, ierr)

      if(iNS.eq.1.or.iNS.eq.2) then
         do i= 1, ncout
         if(d1(i).lt.dvtest) d1(i)=d1(i)+diff
         work(i) = d1(i)
         enddo
      else
         do i=1, ncout
         work(i) = d2(i)
         enddo
      endif

cu      call growf4('w', 18, jrow, work, ncout, ierr)
      call grow('w', 18, jrow, work, ncout, ierr)
450   continue
      
c apply weighting in the overlap region

455   jlap=0
      do 600 j=jrso,jreo
      jrow=j
      jlap=jlap+1
cu      call growf4('r', 13, jlap, work(icso), nclap, ierr)
cu      call growf4('r', 14, jrow, d1, ncout, ierr)
cu      call growf4('r', 15, jrow, d2, ncout, ierr)
      call grow('r', 13, jlap, work(icso), nclap, ierr)
      call grow('r', 14, jrow, d1, ncout, ierr)
      call grow('r', 15, jrow, d2, ncout, ierr)

      if(icso.eq.1) go to 480
        if(iEW.eq.1.or.iEW.eq.2) then
           do i=1,icso-1
           if(d1(i).lt.dvtest) d1(i)=d1(i)+diff
           work(i)= d1(i)
           end do
        else
           do i=1,icso-1
           work(i)= d2(i)
           end do
        endif
c
c overlap region
c
480     ilap=0
        do 490 i=icso, iceo
        if(d1(i).lt.dvtest) then
         if(d2(i).lt.dvtest) then
            wt0= work(i)
            wt1= 1.e0 - wt0
            work(i) = (d1(i)+diff)*wt0 + d2(i)*wt1
         else
            work(i)=d1(i)+diff
         endif
        else
         if(d2(i).lt.dvtest) then
            work(i) = d2(i)
         else
            work(i)=dval
         endif
        endif
490     continue
c
c rest of columns if applicable
c
        if(iceo.eq.ncout) go to 560
        if(iEW.eq.0.or.iEW.eq.2) then
           do i=iceo+1, ncout
           if(d1(i).lt.dvtest) d1(i)=d1(i)+diff
           work(i)= d1(i)
           end do
        else
           do i=iceo+1, ncout
           work(i)= d2(i)
           end do
        endif

cu560   call growf4('w', 18, jrow, work, ncout, ierr)
560   call grow('w', 18, jrow, work, ncout, ierr)

600   continue
c
c rest of rows if not in overlap region
c
      if(jreo.eq.nrout) go to 660
      do 650 j= jreo+1, nrout
      jrow=j
cu      call growf4('r', 14, jrow, d1, ncout, ierr)
cu      call growf4('r', 15, jrow, d2, ncout, ierr)
      call grow('r', 14, jrow, d1, ncout, ierr)
      call grow('r', 15, jrow, d2, ncout, ierr)

      if(iNS.eq.0.or.iNS.eq.2) then
         do i= 1, ncout 
         if(d1(i).lt.dvtest) d1(i)=d1(i)+diff
         work(i) = d1(i)
         end do
      else
         do i= 1, ncout 
         work(i) = d2(i)
         end do
      endif
cu      call growf4('w', 18, jrow, work, ncout, ierr)
      call grow('w', 18, jrow, work, ncout, ierr)
650   continue

660   call gclose(14, 'delete')
      call gclose(15, 'delete')
      call gclose(13, 'keep')
      call gclose(18, 'keep')
	stop
91	print *,' error in mcplug '
c        stop
	end
c***********************************************************************
cc
cu      subroutine mcplug(nim,idv,jdv,ier,zg,iqd,nwrk,nc,nr)
      subroutine mcplug(nim,idv,jdv,ier,zg[huge],iqd[huge],nwrk,nc,nr)
c      dimension n(4)
      dimension zg(nwrk),iqd(nwrk)
c      common /plug1/ zg(100000)
c      common /plug2/ iqd(100000)
c  assuming dimension of wz is maxcol Tien
      dimension wz(1000)
cv      data dv/'37777677777'o/,p2/'m-c ','plug'/
      data dv/1.0e+38/
c      data dvtest/1.0e+36/
c      nwrk=100000
      ier=0
c      read(idv)id,p,nc,nr,nz,xo,dx,yo,dy
      nn=nc*nr
      if(nn.gt.nwrk) go to 99
      do 1 i=1,nn
1        iqd(i)=-1
      ndx=1
      do 2 i=1,nr
      jrow=i
cu      call growf4('r', idv, jrow, zg(ndx), nc, ie)
      call grow('r', idv, jrow, zg(ndx), nc, ie)
c         call rowio(nc,zg(ndx),-1,idv,jdv,ie)
         if(ie.eq.1) go to 99
2        ndx=ndx+nc
      do 3 i=1,nn
         if(zg(i).ge.1.e29) iqd(i)=0
3        continue
c
c  gridr and curvmn are from the minc program
c
20    call gridr(nc,nr,zg)
c      go to 5
      na=0
      do 21 i=1,nn-nc
      if(iqd(i).eq.0) then
        if(zg(i).gt.0.49.and.zg(i).lt.0.51) then
c          iqd(i)=-1
          na=na+1
        else if (zg(i).ne.zg(i+1)) then
          temp=0.5*(zg(i)+zg(i+1))
          if(temp.gt.0.49.and.temp.lt.0.51) then
            zg(i)=temp
c            iqd(i)=-1
            na=na+1
          else
            zg(i)=dv
          endif
        else if (zg(i).ne.zg(i+nc)) then
          temp=0.5*(zg(i)+zg(i+nc))
          if(temp.gt.0.49.and.temp.lt.0.51) then
            zg(i)=temp
c            iqd(i)=-1
            na=na+1
          else
            zg(i)=dv
          endif
        else
          zg(i)=dv
        endif
      endif
21    continue
      print *,na,' unknown weights replaced by 0.5'
c      go to 5
      call gridr(nc,nr,zg)
      na=0
      do 22 i=1,nn-nc
      if(iqd(i).eq.0) then
        if((zg(i).gt.0.24.and.zg(i).lt.0.26).or.(zg(i).gt.0.74.and.
     1  zg(i).lt.0.76)) then
c          iqd(i)=-1
          na=na+1
        else if (zg(i).ne.zg(i+1)) then
          temp=0.5*(zg(i)+zg(i+1))
          if((temp.gt.0.24.and.temp.lt.0.26).or.(temp.gt.0.74.and.
     1    temp.lt.0.76)) then
            zg(i)=temp
c            iqd(i)=-1
            na=na+1
          else
            zg(i)=dv
          endif
        else if (zg(i).ne.zg(i+nc)) then
          temp=0.5*(zg(i)+zg(i+nc))
          if((temp.gt.0.24.and.temp.lt.0.26).or.(temp.gt.0.74.and.
     1    temp.lt.0.76)) then
            zg(i)=temp
c            iqd(i)=-1
            na=na+1
          else
            zg(i)=dv
          endif
        else
          zg(i)=dv
        endif
      endif
22    continue
      print *,na,' unknown weights replaced by 0.25 or 0.75'
      call gridrold(nc,nr,zg)
      if(nim.eq.0) go to 5
23      eps=0.
      call curvmn(zg,iqd,wz,nc,nr,eps,nim,a,b,ni)
c5     write(jdv)id,p2,nc,nr,nz,xo,dx,yo,dy
5      ndx=1
      do 4 i=1,nr
        jrow=i
cu         call growf4('w', jdv, jrow, zg(ndx), nc, ie)
         call grow('w', jdv, jrow, zg(ndx), nc, ie)
c         call rowio(nc,zg(ndx),0,idv,jdv,ie)
4        ndx=ndx+nc
      return
99    ier=1
      return
      end
c****************************************************************************
c function amode works like amod should
      function amode(x,y)
      if(y.eq.0.0) then
        amode=0.0
      else
        amode=x-anint(x/y)*y
      endif
      return
      end
