c   USGS Computer Program "MAGMRG"; version 1.000
c   Technical Contact:  Ron Sweeney
c
      dimension id1(14),id2(14),pgm1(2),pgm2(2),data1(2001),data(2000)
      dimension dvals(2001),idvals(2001),dum(2001),iscsp4(2,4,2000)
      dimension rows(2000,8),spl1(14,100),spl2(14,100),a(4,10)
      dimension iscspl(2,2000),iscsp2(8,2000)
      character*20 infl1,infl2,itmpf,iofil,insw1,insw2,in1,in2
      character*20 iout1,iout2,iosw1
      character*10 idisp,idel
      character*4 ie,iw
      character*1 dir,idiag,idedge
      logical normal
      equivalence (data1(2),data(1)),(dvals(1),idvals(1))
      equivalence (iscsp2(1,1),iscsp4(1,1,1),rows(1,1))
      equivalence (spl1(1,1),iscspl(1,1)),(spl2(1,1),iscspl(1,701))
      data dvals(1)/0.0/,(dvals(i),i=2,2001)/2000*.1701412e39/
c
c   Write the ID-Stamp to default output device.
c
      write(6,15)
   15 format(' USGS Computer Program "MAGMRG"; version 1.000',/,
     1       ' Technical Contact:  Ron Sweeney')
c
      j4=4
      odval=1.0e30
      data1(1)=0.0
      ie='e(n)'
      iw='w(s)'
      maxd=8000
      insw1='swtch1.tmp'
      insw2='swtch2.tmp'
      iosw1='swtch3.tmp'
      idisp='keep'
      idel='delete'
c     itmpf='scr.tmp'
      normal=.true.
      inscl=0
      insc=0
      write(6,7)
    7 format(2x,'Enter direction for splining, ''n-s'' or ''e-w'':  ',$)
      read(5,2) dir
      write(6,1)
    1 format(2x,'Enter input filename # 1 (west-south):  ',$)
      read(5,2) in1
    2 format(a)
      write(6,3)
    3 format(2x,'Enter input filename # 2 (east-north):  ',$)
      read(5,2) in2
      write(6,4)
    4 format(2x,'Enter output filename:  ',$)
      read(5,2) iofil
c
      itmpf=iofil
c
      open(unit=12,file=itmpf,status='unknown',form='unformatted')
c
c * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
c     write(6,6)
c   6 format(2x,'enter # of data pts to eliminate from w(s) grid,',/,
c    12x,'gap for splining (normally = 0 2),',/,2x,'''x'' if ',
c    2'secondary boundary > 100 pts or jagged,',/,2x,'and ''d'' if ',
c    3'pts to be eliminated from rt data edge.',/,2x,'(2i2,2a1):')
c     read(5,9) ingapl,ingap,idiag,idedge
c   9 format(2i2,2a1)
c     if(idiag.eq.'x'.or.idiag.eq.'X') normal=.false.
c
c * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
      write(6,6)
    6 format(2x,'Enter # of data pts to eliminate from w(s) grid,',/,
     12x,'& gap for splining (e.g. 0 2):  ',$)
      read(5,*) ingapl,ingap
      idiag='x'
      idedge='d'
      normal=.false.
      if(normal) go to 14
c
c * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
c     write(6,12)
c  12 format(2x,'enter secondary splining parms (normally = 0 2) (2i2):'
c    1)
c     read(5,13) inscl,insc
c  13 format(2i2)
c
c * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
      write(6,12)
   12 format(2x,'Enter secondary splining parms (e.g. 0 2):  ',$)
      read(5,*) inscl,insc
   14 iout1=iofil
      iout2=itmpf
      infl1=in1
      infl2=in2
      if(dir.eq.'e'.or.dir.eq.'w'.or.dir.eq.'E'.or.dir.eq.'W') go to 5
      call switch(infl1,insw1,maxd,rows(1,1),rows(1,5),10,11,idisp,*960)
      infl1=insw1
      call switch(infl2,insw2,maxd,rows(1,1),rows(1,5),10,11,idisp,*970)
      infl2=insw2
      iout1=iosw1
      iout2=iofil
      idisp=idel
    5 open(unit=10,file=infl1,status='old',form='unformatted')
      read(10,end=900) id1,pgm1,nc1,nr1,nz1,xo1,dx1,yo1,dy1
      open(unit=11,file=infl2,status='old',form='unformatted')
      read(11,end=902) id2,pgm2,nc2,nr2,nz2,xo2,dx2,yo2,dy2
c
c   initialize program parameters.
c
      do 11 j=1,2000
      do 11 i=1,8
   11 iscsp2(i,j)=-1
      isc=0
      irspl=0
      nspl=10
      iro=0
      ir1=0
      ir2=0
      tol=0.00005
c     rndp=0.00005
c     rndm=-0.00005
c
c   isav1, isav2 are the min and max lats (n-s) or longs (e-w)
c    for secondary splining.
c
      isav1=0
      isav2=0
      ic1a=0
      ic1b=0
      ic2a=0
      ic2b=0
      xe1=xo1+(nc1-1)*dx1
      xe2=xo2+(nc2-1)*dx2
      yn1=yo1+(nr1-1)*dy1
      yn2=yo2+(nr2-1)*dy2
c     rnd=rndp
c   ********************************************************************
cc     if(xo1.lt.xo2) go to 60
c      if(xo1-xo2.lt.tol) go to 60
c
c   here xo2 .le. xo1.  initialize unit, cols, and rows for leftmost
c    grid.
c
c     iul=11
c     icl=nc2
c     irl=nr2
c     iur=10
c     icr=nc1
c     irr=nr1
c     xout=xo2
c     yol=yo2
c     yor=yo1
c     ynl=yn2
c     ynr=yn1
c     dyl=dy2
c     dyr=dy1
cc    if(xe2.gt.xo1) rnd=rndm
cc    idfx=(xo1-xe2+rnd)/dx2
c     dfx=(xo1-xe2)/dx2
c     hlf=sign(0.5,dfx)
c     idfx=dfx+hlf
cc    if(xe2.gt.xe1) go to 55
c     if(xe1-xe2.lt.tol) go to 55
c     ncout=nc1+nc2+idfx-1
c     go to 70
c  55 ncout=nc2
c     go to 70
c
c   here xo1 .lt. xo2.  initialize unit, cols, and rows for rightmost
c    grid.
c
c  60 iul=10
      iul=10
      icl=nc1
      irl=nr1
      iur=11
      icr=nc2
      irr=nr2
      xout=xo1
      yol=yo1
      yor=yo2
      ynl=yn1
      ynr=yn2
      dyl=dy1
      dyr=dy2
      iovx=0
      iovx1=1
c     if(xo1.le.xo2) go to 61
      if(xo1-xo2.le.tol) go to 61
      xout=xo2
c     iovx=(xo1-xo2+rnd)/dx1
      ovx=(xo1-xo2)/dx1
      hlf=sign(0.5,ovx)
      iovx=ovx+hlf
      iovx1=iovx+1
      icl=icl+iovx
c  61 if(xe1.gt.xo2) rnd=rndm
c     idfx=(xo2-xe1+rnd)/dx1
   61 dfx=(xo2-xe1)/dx1
      hlf=sign(0.5,dfx)
      idfx=dfx+hlf
c     if(xe1.gt.xe2) go to 65
      if(xe2-xe1.lt.tol) go to 65
      ncout=icl+nc2+idfx-1
      go to 70
   65 ncout=icl
c
c   here the size of output array calculated in x-direction.  now
c    calculate in y-direction.
c
c  70 rnd=rndp
c     if(yo1.lt.yo2) go to 80
   70 if(yo1-yo2.lt.tol) go to 80
c
c   here yo2 .le. yo1.
c
      yout=yo2
c     if(yn2.gt.yo1) rnd=rndm
c     idfy=(yo1-yn2+rnd)/dy2
      dfy=(yo1-yn2)/dy2
      hlf=sign(0.5,dfy)
      idfy=dfy+hlf
c     if(yn2.gt.yn1) go to 75
      if(yn1-yn2.lt.tol) go to 75
      nrout=nr1+nr2+idfy-1
      go to 90
   75 nrout=nr2
      go to 90
c
c   here yo1 .lt. yo2.
c
   80 yout=yo1
c     if(yn1.gt.yo2) rnd=rndm
c     idfy=(yo2-yn1+rnd)/dy1
      dfy=(yo2-yn1)/dy1
      hlf=sign(0.5,dfy)
      idfy=dfy+hlf
c     if(yn1.gt.yn2) go to 85
      if(yn2-yn1.lt.tol) go to 85
      nrout=nr1+nr2+idfy-1
      go to 90
   85 nrout=nr1
c
c   now output the header record on unit 12 for output data set.
c
   90 iposl=icl+idfx
      iposr=icr+idfx
      iclp1=icl+1
      iposl1=iposl-1
      icout=ncout-iposl1
      icrl1=icr+iposl1
      icrp1=icr+1
      write(12) id1,pgm1,ncout,nrout,nz1,xout,dx1,yout,dy1
      ys2=yor-dyr
      ys1=yol-dyl
c     if(yol.le.yor) go to 101
      if(yol-yor.le.tol) go to 101
      ipsn=iposl1
c
c   here the rightmost grid is lowest.
c
  108 ys2=ys2+dyr
c     if(yol.le.ys2) go to 102
      if(yol-ys2.le.tol) go to 102
      call rdbin(iur,dum,1,data,icr,*800)
      ir2=ir2+1
      if(ncout.eq.icrl1) go to 128
      do 125 i=icrp1,icout
  125 data(i)=dvals(2)
  128 call iorow(12,dvals,iposl,data,icout)
      iro=iro+1
      go to 108
c
c   here grid 2 (rt) lower boundary .ge. (higher) than grid 1 (left)
c    lower boundary.
c
  101 ys2=ys2+dyr
      ipsn=iclp1
  102 ys1=ys1+dyl
c
c   ys1 & ys2 now represent the next rows to be read.
c
c     if(ys1.ge.ys2) go to 106
      if(ys2-ys1.le.tol) go to 106
      call rdbin(iul,dum,1,data(iovx1),nc1,*802)
      if(iovx.eq.0) go to 105
c
c   pad left side of grid 1 with dvals, if necessary.
c
      do 104 i=1,iovx
  104 data(i)=dvals(2)
  105 ir1=ir1+1
      if(iposr.gt.0) then
      call iorow(12,data1,iclp1,dvals(2),iposr)
      else
      call iorow1(12,data1,iclp1)
      end if
      iro=iro+1
      go to 102
c
c   here we have 2 rows at same lat, with or without overlap.
c
  106 if(normal) go to 107
      isav1=iro
      isc=isc+1
      iscspl(1,isc)=ipsn
      iscspl(2,isc)=ipsn
  107 call rdbin(iul,dum,1,data(iovx1),nc1,*804)
      if(iovx.eq.0) go to 121
c
c   pad left side of grid 1 with dvals, if necessary.
c
      do 120 i=1,iovx
  120 data(i)=dvals(2)
  121 ir1=ir1+1
c
c   search for first good pt at rightmost boundary of data set.
c
      ist=icl
      if(idedge.ne.'d'.and.idedge.ne.'D') ist=ist-ingapl
      istr=ist+1
      do 112 i1r=1,ist
      i1=istr-i1r
      if(data(i1).lt.odval) go to 115
  112 continue
c
c   here all dvals found for leftmost data set.  output rightmost row
c    and continue.
c
  126 if(normal) go to 119
      isc=isc+1
      iscspl(1,isc)=iposl1
      iscspl(2,isc)=iposl1
  119 call rdbin(iur,dum,1,data,icr,*808)
      ir2=ir2+1
      l1=iposl
      l2=icr
      if(ncout.eq.icrl1) go to 142
      do 122 i=icrp1,icout
  122 data(i)=dvals(2)
      l2=icout
      go to 142
c
c   here calc posn to start inputting good data values from rightmost
c    grid row.
c
  115 if((idedge.ne.'d'.and.idedge.ne.'D').or.ingapl.eq.0) go to 127
      ist=i1-ingapl
      if(ist.le.0) go to 126
      istr=ist+1
      do 129 i1r=1,ist
      i1=istr-i1r
      if(data(i1).lt.odval) go to 127
  129 continue
      go to 126
  127 igap=ingap
      irst=i1-iposl+igap+1
      ilst=i1-1
      ilpos=i1
      if(irst.ge.1) go to 114
c
c   here the posn of the first pt in rightmost grid to leave gap of 2
c    intervals is < 1.  must set posn = 1 and adjust gap.
c
      igap=igap-irst+1
      irst=1
c
c   now read in data from rightmost grid, starting with pt irst.  this
c    point belongs in posn ilpos+igap.
c
  114 ipos=ilpos+igap
      idim=icr-irst+1
      irspl=iro+2
      if(idim) 109,109,111
c
c   here row is complete (row from rt grid is contained in row from
c    left grid + the splining gap) - no splining necessary.  let npts =
c    igap + idim - 1.  if npts > 0, then read these pts from rightmost
c    end of grid 2 row into data array.  if npts .le. 0, read grid 2
c    row into dummy array.
c
  109 npts=igap+idim-1
      ipos=ilpos+1
      if(npts.le.0) go to 110
      irst=irst-igap+1
      if(irst.ge.1) go to 113
c
c   here irst < 1.  adjust ipos, npts, and irst accordingly, and insert
c    dvals in data array between left and right grid rows.
c
      iend=ipos-irst
      npts=npts+irst-1
      irst=1
      do 118 idv=ipos,iend
  118 data(idv)=dvals(2)
      ipos=iend+1
  113 call rdbin(iur,dum,irst,data(ipos),npts,*808)
      ir2=ir2+1
      if(normal) go to 147
      isc=isc+1
      ilpo1=ilpos+1
      iscspl(1,isc)=ilpo1
      iscspl(2,isc)=ilpo1
      go to 141
c
c   here read grid 2 row into dummy array.
c
  110 call rdbin(iur,dum,1,dum(2),icr,*808)
      ir2=ir2+1
      if(normal) go to 147
      isc=isc+1
      iscspl(1,isc)=ipos
      iscspl(2,isc)=ipos
      go to 141
  111 call rdbin(iur,dum,irst,data(ipos),idim,*808)
      ir2=ir2+1
c
c   insert dvals in gap between left and right rows.
c
      ilpo1=ilpos+1
      l=ipos-1
      do 103 i=ilpo1,l
  103 data(i)=dvals(2)
      do 116 i2=ipos,ncout
      if(data(i2).lt.odval) go to 117
  116 continue
c
c   here data all dvals beyond irst.  no splining can be done.  output
c    data for leftmost grid and continue.
c
      go to 144
c
c   set irpos = first non-dval to right of gap.
c
  117 irpos=i2
      igap=igap+i2-ipos
      if(igap.gt.nspl) go to 144
      irst=i2+1
      ipdv=0
      call splln(data,ncout,igap,ilpos,irpos,ilst,irst,a,ipdv)
c
c   now output the splined data row.
c
      ii1=ilpos+1
      ii2=irpos-1
      if(.not.normal) go to 137
      ispn=1
c   ********************************************************************
      if(ic1a.ne.0) go to 137
      isav1=iro
c     if(yol.gt.yor) go to 135
      if(yor-yol.lt.tol) go to 135
      ic1a=ii1
      ic1b=ii2
      if(ic1b.lt.icl) ic1b=icl
      go to 137
  135 ic1a=iposl
      ic1b=ii2
c   ********************************************************************
  137 if(ipdv.ne.0) go to 143
      is=0
      do 140 i=ii1,ii2
      is=is+1
      val=a(4,5)
      do 138 jr=1,3
      j=j4-jr
  138 val=val*is+a(j,5)
      data(i)=val
  140 continue
  143 if(normal) go to 141
      isc=isc+1
      iscspl(1,isc)=ii1
      iscspl(2,isc)=ii2
      go to 141
c
c   here no splining to be done.  gap too large.
c
  144 ilpo1=ilpos+1
      if(normal) go to 141
      isc=isc+1
      iscspl(1,isc)=ilpo1
      iscspl(2,isc)=ilpo1
      go to 141
  147 ispn=0
      if(ic1a.ne.0) go to 141
      isav1=iro
c     if(yol.gt.yor) go to 148
      if(yor-yol.lt.tol) go to 148
      ic1a=ilpos+1
      ic1b=icl
      if(ic1a.gt.ic1b) isav1=0
      go to 141
  148 if(ilpos.gt.iposl) go to 150
      ic1a=ilpos
      go to 151
  150 ic1a=iposl
  151 ic1b=ilpos
  141 l1=1
      l2=ncout
  142 call iorow(12,dvals,l1,data,l2)
      iro=iro+1
      if(ir1.eq.irl) go to 175
      if(ir2.eq.irr) go to 180
      go to 107
c
c   here pad data rows with data from rightmost grid.
c
  175 isav2=irspl
      ipsn=iposl1
  176 if(ir2.ge.irr) go to 178
      call rdbin(iur,dum,1,data,icr,*800)
      ir2=ir2+1
      if(ncout.eq.icrl1) go to 166
      do 163 i=icrp1,icout
  163 data(i)=dvals(2)
  166 call iorow(12,dvals,iposl,data,icout)
      iro=iro+1
      go to 176
c
c   here pad data rows with data from leftmost grid.
c
  180 isav2=irspl
      ipsn=iclp1
  179 call rdbin(iul,dum,1,data(iovx1),nc1,*802)
      if(iovx.eq.0) go to 187
c
c   pad left side of grid 1 with dvals, if necessary.
c
      do 186 i=1,iovx
  186 data(i)=dvals(2)
  187 ir1=ir1+1
      if(iposr.gt.0) then
      call iorow(12,data1,iclp1,dvals(2),iposr)
      else
      call iorow1(12,data1,iclp1)
      end if
      iro=iro+1
      if(ir1.lt.irl) go to 179
c
c   here check for jagged secondary splining boundary.
c
  178 if(normal) go to 181
c   ********************************************************************
c
c   here set last value in iscspl array .  then calc new array (iscsp2)
c    to determine n-s splining.
c
      isc=isc+1
      iscspl(1,isc)=ipsn
      iscspl(2,isc)=ipsn
      isplin=0
      lsc=1
      lst1=iscspl(1,1)
      lst2=iscspl(2,1)
      do 250 ipos=2,isc
      isc1=iscspl(1,ipos)
      isc2=iscspl(2,ipos)
c
c   check for -1 pad, set if gap > 10.
c
c     if(lst1.eq.-1.and.lst2.eq.-1) go to 249
c     if(isc1.eq.-1.and.isc2.eq.-1) go to 249
      if(isc1-lst1) 201,225,210
c
c   check for 1st row; if so, spline down.  all other splines are up.
c
  201 if(ipos.gt.2) go to 202
      iscsp2(1,1)=isc1
      iscsp2(2,1)=lst1-1
      if(isav1.le.0) go to 225
      go to 227
c
c   here calc 1st & 2nd values for spline on last row.
c
  202 iscsp2(5,lsc)=isc1
      if(isc2.lt.lst1) go to 204
      iscsp2(6,lsc)=lst1-1
      go to 224
  204 iscsp2(6,lsc)=isc2
      go to 224
c
c   check for 1st row; if so, spline down.  all other splines are up.
c
  210 if(ipos.gt.2) go to 212
      iscsp2(1,1)=lst2+1
      iscsp2(2,1)=isc2
      if(isav1.le.0) go to 249
      go to 248
c
c   here calc 1st & 2nd values for spline on current row.
c
  212 iscsp2(1,ipos)=lst1
      iscsp2(2,ipos)=isc1-1
  224 if(ipos.eq.isc.and.isav2.gt.nrout) go to 225
  227 isplin=1
c
c   here set secondary spline values based on end spline pts from last
c    2 rows.
c
  225 if(isc2-lst2) 226,249,240
c
c   note:  don't have to check for ipos=2 (1st row to spline.)
c          already handled above.
c
  226 iscsp2(3,ipos)=isc2+1
      iscsp2(4,ipos)=lst2
      go to 248
c
c   calc 1st & 2nd values for last row.
c
  240 if(ipos.gt.2) go to 246
      iscsp2(7,1)=lst2+1
      iscsp2(8,1)=isc2
      if(isav1.le.0) go to 249
      go to 247
  246 if(lst2.lt.isc1) go to 241
      iscsp2(7,lsc)=lst2+1
      go to 244
  241 iscsp2(7,lsc)=isc1
  244 iscsp2(8,lsc)=isc2
  248 if(ipos.eq.isc.and.isav2.gt.nrout) go to 249
  247 isplin=1
  249 lsc=ipos
      lst1=isc1
      lst2=isc2
  250 continue
c     write(6,252) ((iscspl(i,j),i=1,2),j=1,isc)
c 252 format(2x,'****iscspl:',(/,2x,10i6))
c     write(6,253) ((iscsp2(i,j),i=1,8),j=1,isc)
c 253 format(2x,'****iscsp2:',(/,2x,16i4))
      if(isplin.eq.0) go to 810
      write(6,251)
  251 format(2x,'n-s splining necessary.  switching rows to columns ',
     1'to accommodate.')
      close(10,status=idisp)
      close(11,status=idisp)
        close(12)
c
c   save iscsp2 array on scratch disk until switch requiring rows buffer
c    has completed.
c
      open(unit=12,file=insw1,status='unknown',form='unformatted')
      write(12) iscsp2
      close(12)
c
c   note:  here splined data set is file scr.tmp; switch to swtch3.tmp.
c
      call switch(itmpf,iosw1,maxd,rows(1,1),rows(1,5),10,11,idel,*990)
c
c   now restore iscsp2 array.
c
      open(unit=12,file=insw1,status='old',form='unformatted')
      read(12) iscsp2
      close(12,status='delete')
c
c   here secondary splining is necessary for e-w or n-s data set.  saved
c    row coordinates are now translated column coordinates.
c
      open(unit=10,file=iosw1,status='old',form='unformatted')
      open(unit=11,file=iout2,status='unknown',form='unformatted')
      read(10) id1,pgm1,ncout,nrout,nz1,xout,dxout,yout,dyout
      write(11) id1,pgm1,ncout,nrout,nz1,xout,dxout,yout,dyout
      do 280 jcol=1,nrout
      call rdbin(10,dum,1,data,ncout,*995)
c
c   now scan entire iscsp2 array for each row.  -1 indicates no splining
c    necessary.  if 2 consecutive columns require splining, set gap = 3;
c    otherwise, perform splines one at a time.
c
      ione=0
      do 278 jrow=1,isc+1
      jro=isav1+jrow-1
      if(jrow.gt.isc) go to 265
      do 277 i4=1,4
      if(iscsp4(1,i4,jrow).eq.-1.or.iscsp4(2,i4,jrow).eq.-1) go to 277
      if(iscsp4(1,i4,jrow).gt.jcol.or.jcol.gt.iscsp4(2,i4,jrow))
     1go to 277
c
c   here jro'th column requires splining.  check for two consecutive
c    values.
c
      if(jro.ge.ncout) go to 278
      if(ione.eq.1) go to 268
      ione=1
      go to 278
  277 continue
c
c   here have 1 value, jro-1, that requires splining.
c
  265 if(ione.eq.0) go to 278
      isrch1=jro-2
      isrch2=jro
      go to 269
c
c   here have two values, jro-1 & jro, to spline across.
c
  268 isrch1=jro-2
      isrch2=jro+1
  269 jro1=jro-1
      if(jro1.ne.isav1) go to 281
c
c   here we are splining along min row.  allow change to secondary
c    splining parms.
c
      isrch2=isrch2+inscl
      isrch1=min0(isrch2-insc,isrch1)
      go to 282
  281 if(jro.ne.isav2.and.jro1.ne.isav2) go to 282
c     write(6,283) jro,isc,isav2,isav1,isrch1,isrch2,jrow
c 283 format(2x,'******',7i4)
c
c   here we are splining along max row.  allow change to secondary
c    splining parms.
c
      isrch1=isrch1-inscl
      isrch2=max0(isrch1+insc,isrch2)
  282 ione=0
      if(isrch1.le.0.or.isrch2.gt.ncout) go to 278
      isrchr=isrch1+1
      do 270 i1r=1,isrch1
      i1=isrchr-i1r
      if(data(i1).lt.odval) go to 272
  270 continue
c
c   no spline to be performed.
c
      go to 278
  272 ilpos=i1
      ilst=i1-1
      do 273 i2=isrch2,ncout
      if(data(i2).lt.odval) go to 274
  273 continue
c
c   no spline to be performed.
c
      go to 278
  274 irpos=i2
      irst=i2+1
      igap=irpos-ilpos
      if(igap.gt.nspl) go to 278
      ipdv=0
      call splln(data,ncout,igap,ilpos,irpos,ilst,irst,a,ipdv)
c
c   now insert splined values.
c
      ii1=ilpos+1
      ii2=irpos-1
      if(ipdv.ne.0) go to 278
      is=0
      do 276 i=ii1,ii2
      is=is+1
      val=a(4,5)
      do 275 jr=1,3
      j=j4-jr
  275 val=val*is+a(j,5)
      data(i)=val
  276 continue
  278 continue
      call iorow(11,dvals,1,data,ncout)
  280 continue
c
c   here, if e-w, switch back from scr.tmp file to iofil.
c
      go to 955
c   ********************************************************************
c
c   here all data has been output on unit 12.  now do n-s splines.
c    note:  for rows isav1 & isav2, the splines must be done from
c    columns ic1a-ic1b and ic2a-ic2b, respectively.
c
  181 if(isav2.le.1.or.isav2.ge.iro) go to 172
      if(ispn) 152,152,158
c 152 if(ynl.lt.ynr) go to 155
  152 if(ynl-ynr.lt.tol) go to 155
      ic2a=ilpos+1
      ic2b=icl
      if(ic2a.gt.ic2b) isav2=0
      go to 177
  155 if(ilpos.gt.iposl) go to 156
      ic2a=ilpos
      go to 157
  156 ic2a=iposl
  157 ic2b=ilpos
      go to 177
c 158 if(ynl.lt.ynr) go to 162
  158 if(ynl-ynr.lt.tol) go to 162
      ic2a=ii1
      ic2b=ii2
      if(ic2b.lt.icl) ic2b=icl
      go to 177
  162 ic2a=iposl
      ic2b=ii2
      go to 177
  172 if(isav1.gt.1.and.isav1.lt.iro) go to 177
c
c   rename file and print message.
c
      write(6,174)
  174 format(2x,'no n-s splining necessary.  good end of job.')
      go to 810
  177 close(10,status=idisp)
      close(11,status=idisp)
      close(12)
      open(unit=10,file=itmpf,status='old',form='unformatted')
      open(unit=11,file=iout1,status='unknown',form='unformatted')
      read(10) id1,pgm1,ncout,nrout,nz1,xout,dx1,yout,dy1
      write(11) id1,pgm1,ncout,nrout,nz1,xout,dx1,yout,dy1
      i1=isav1-7
      i2=isav2-7
      ie1=isav1+7
      ie2=isav2+7
      ir=0
      ir1=0
      ir2=0
c
c   irout is needed to determine if iorow called.  irow = # of values in
c    rows array.
c
      irout=0
      irow=0
  182 call rdbin(10,dum,1,data,ncout,*850)
      ir=ir+1
      if(isav1.le.1) go to 190
      if(ir.lt.i1) go to 199
      if(ir.gt.ie1) go to 189
      if(ir.eq.isav1) go to 184
c
c   ir1 = posn for current row stored for s spline.
c
      ir1=ir1+1
      ipos=0
      do 183 i=ic1a,ic1b
      ipos=ipos+1
  183 spl1(ir1,ipos)=data(i)
      if(ir.gt.isav1) go to 184
      call iorow(11,dvals,1,data,ncout)
      irout=irout+1
      go to 189
c
c   save row in rows array.
c
  184 irow=irow+1
      do 185 i=1,ncout
  185 rows(i,irow)=data(i)
      irout=irout+1
      if(ir.ne.ie1) go to 189
c
c   here perform s spline and output all rows .lt. isav2.  then transfer
c    those values .ge. isav2 to beginning of rows array.
c
      call rowspl(ic1a,ic1b,isav1,nrout,ir1,spl1,rows)
c
c   now output all rows .lt. isav2.
c
      imov=0
      ipos=isav1
      do 171 i=1,irow
      if(ipos.ge.isav2) go to 167
      call iorow(11,dvals,1,rows(1,i),ncout)
      go to 171
  167 imov=imov+1
      do 168 j=1,ncout
  168 rows(j,imov)=rows(j,i)
  171 ipos=ipos+1
      ir1=0
      irow=imov
c
c   now do n spline.
c
  189 if(isav2.ge.nrout) go to 198
  190 if(ir.lt.i2) go to 198
      if(ir.gt.ie2) go to 199
      if(ir.eq.isav2) go to 196
c
c   save points in spl2 array.
c
      ir2=ir2+1
      ipos=0
      do 191 i=ic2a,ic2b
      ipos=ipos+1
  191 spl2(ir2,ipos)=data(i)
c     write(6,1005) ir2,irow,ir,irout,ic2a,ic2b,ic1a,ic1b
c1005 format(2x,'****',8i4)
      if(ir.lt.ie2) go to 194
c
c   here save row in rows array.
c
      irow=irow+1
      do 195 i=1,ncout
  195 rows(i,irow)=data(i)
      irout=irout+1
c
c   here perform n spline and output all rows.
c
      call rowspl(ic2a,ic2b,isav2,nrout,ir2,spl2,rows)
c
c   now output all rows.
c
      do 193 i=1,irow
      call iorow(11,dvals,1,rows(1,i),ncout)
  193 continue
      ir2=0
      irow=0
      go to 182
  194 if(ir.lt.isav2) go to 198
c
c   save row in rows array.
c
  196 if(irout.eq.ir) go to 182
      irow=irow+1
      do 197 i=1,ncout
  197 rows(i,irow)=data(i)
      irout=irout+1
      go to 182
c
c   if iorow already called for current row, or current row already
c    stored in rows array, do not call iorow again.
c
  198 if(irout.eq.ir) go to 182
  199 call iorow(11,dvals,1,data,ncout)
      irout=irout+1
      go to 182
c
c   now handle end of file at 182.
c
  850 if(irow-1) 950,870,880
c
c   here output the one saved row, as is, from the rows array.
c
  870 call iorow(11,dvals,1,rows(1,1),ncout)
      go to 950
  880 if(ir1.eq.0) go to 890
c
c   here rows are stored in rows array and points are ready for
c    splining.  spline and output and set rows array for n spline.
c    (ir should be > isav1).
c
c   now do n spline; output all data.
c
      call rowspl(ic1a,ic1b,isav1,nrout,ir1,spl1,rows)
      imov=0
      ipos=isav1
      do 884 i=1,irow
      if(ipos.ge.isav2) go to 882
      call iorow(11,dvals,1,rows(1,i),ncout)
      go to 884
  882 imov=imov+1
      do 883 j=1,ncout
  883 rows(j,imov)=rows(j,i)
  884 ipos=ipos+1
      if(imov.eq.0) go to 950
      irow=imov
  890 if(ir2.eq.0) go to 950
      call rowspl(ic2a,ic2b,isav2,nrout,ir2,spl2,rows)
      do 892 i=1,irow
      call iorow(11,dvals,1,rows(1,i),ncout)
  892 continue
      go to 950
  900 write(6,901) infl1
  901 format(2x,'end of file reading header on file:',a20)
      go to 999
  902 write(6,901) infl2
      go to 999
  800 write(6,801) iur,ie,ir2
  801 format(2x,'****unexpected end of file on unit ',i2,' for ',
     1a4,' grid.',/,6x,i4,' data records read.',/,6x,'execution ',
     2'terminating....')
  807 if(dir.eq.'e'.or.dir.eq.'w'.or.dir.eq.'E'.or.dir.eq.'W') go to 803
      close(10,status='delete')
      close(11,status='delete')
      close(12)
c
      iout1=itmpf 
      call switch(itmpf,iofil,maxd,rows(1,1),rows(1,5),10,11,idisp,*980)
      go to 999
  803 close(12)
c     call rename(itmpfr,iofilr)
      go to 999
  802 write(6,801) iul,iw,ir1
      go to 807
  804 write(6,805) iul,iw
  805 format(2x,'end of file on unit ',i2,' for ',a4,' grid.',/,6x,
     1'continue reading other grid to complete output grid.')
  806 call rdbin(iur,dum,1,data,icr,*810)
      ir2=ir2+1
      call iorow(12,dvals,iposl,data,icr)
      iro=iro+1
      if(ir2.lt.irr) go to 806
      go to 810
c
c   output row from left grid, then read and output data from the left
c    grid till the end of data set.
c
  808 write(6,805) iur,ie
  809 if(iposr.gt.0) then
      call iorow(12,data1,iclp1,dvals(2),iposr)
      else
      call iorow1(12,data1,iclp1)
      end if
      iro=iro+1
      if(ir1.eq.irl) go to 810
      call rdbin(iul,dum,1,data(iovx1),nc1,*810)
      ir1=ir1+1
      if(iovx.eq.0) go to 809
c
c   pad left side of grid 1 with dvals, if necessary.
c
      do 811 i=1,iovx
  811 data(i)=dvals(2)
      go to 809
  810 close(12)
c     call rename(itmpfr,iout1r)
      if(dir.eq.'e'.or.dir.eq.'w'.or.dir.eq.'E'.or.dir.eq.'W') go to 952
      close(10,status='delete')
      close(11,status='delete') 
      call switch(itmpf,iofil,maxd,rows(1,1),rows(1,5),10,11,idisp,*980)
      go to 949
  950 close(10,status='delete')
      close(11)
c
      if(dir.eq.'e'.or.dir.eq.'w'.or.dir.eq.'E'.or.dir.eq.'W') go to 952
  948 call switch(iout1,iofil,maxd,rows(1,1),rows(1,5),10,11,idisp,*980)
  949 write(6,951) id1,pgm1,nr1,nc1,yo1,dy1,xo1,dx1,id2,pgm2,nr2,nc2,
     1yo2,dy2,xo2,dx2,id1,pgm1,nrout,ncout,yout,dy1,xout,dx1
      go to 999
  952 write(6,951) id1,pgm1,nc1,nr1,xo1,dx1,yo1,dy1,id2,pgm2,nc2,nr2,
     1xo2,dx2,yo2,dy2,id1,pgm1,ncout,nrout,xout,dx1,yout,dy1
  951 format(2x,'end of job.',/,2x,'input grids =:',/,2(/,6x,14a4,/,
     16x,2a4,/,6x,2i5,4f14.2),/,2x,'output grid =:',//,6x,14a4,/,6x,
     22a4,/,6x,2i5,4f14.2)
      go to 999
  955 close(10,status='delete')
      close(11)
      if(dir.ne.'e'.and.dir.ne.'w'.and.dir.ne.'E'.and.dir.ne.'W')
     1go to 957
      call switch(iout2,iofil,maxd,rows(1,1),rows(1,5),10,11,idel,*990)
      write(6,951) id1,pgm1,nc1,nr1,xo1,dx1,yo1,dy1,id2,pgm2,nc2,nr2,
     1xo2,dx2,yo2,dy2,id1,pgm1,nrout,ncout,yout,dyout,xout,dxout
      go to 999
  957 write(6,951) id1,pgm1,nr1,nc1,yo1,dy1,xo1,dx1,id2,pgm2,nr2,nc2,
     1yo2,dy2,xo2,dx2,id1,pgm1,ncout,nrout,xout,dxout,yout,dyout
      go to 999
  960 write(6,961) infl1
  961 format(2x,'****unexpected end of file converting rows to columns',
     1/,6x,'for file ',a20,/,6x,'  execution terminating....')
      close(11,status='delete')
      go to 999
  970 write(6,961) infl2
      close(11,status='delete')
      go to 999
  980 write(6,961) iofil
      close(10,status='delete')
      close(11)
      go to 999
  990 write(6,961) itmpf
      close(11,status='delete')
      go to 999
  995 write(6,996) iosw1
  996 format(2x,'****unexpected end of file on file ',a20,/,6x,'while ',
     1'performing secondary splining.',/,6x,'execution terminating....')
  999 stop
      end

      subroutine iorow(io,dat1,n1,dat2[huge],n2)
      dimension dat1(n1),dat2(n2)
      write(io) dat1,dat2
      return
      end
      subroutine iorow1(io,dat1,n1)
      dimension dat1(n1)
      write(io) dat1
      return
      end
      subroutine rdbin(iu,dum,i1,data[huge],i2,*)
      dimension dum(i1),data(i2)
      read(iu,end=100) dum,data
      return
  100 return 1
      end
      subroutine switch(infil,iofil,maxd,dati[huge],
     1 dato[huge],in,iout,idisp,*)
      dimension dati(1),dato(1),id(14),pgm(2)
      character*20 infil,iofil
      character*10 idisp
        kout=22
      open(unit=in,file=infil,status='old',form='unformatted')
        if (infil .eq. iofil) then
      open(unit=kout,status='scratch',form='unformatted')
        else
        kout=iout
      open (kout,file=iofil,status='unknown',form='unformatted')
        endif
      read(in,end=100) id,pgm,nc,nr,nz,xo,dx,yo,dy
      write(kout) id,pgm,nr,nc,nz,yo,dy,xo,dx
      ncol=maxd/nr
      inc=min0(ncol,nc)
      nrow=maxd/nc
      inr=min0(nrow,nr)
      call rdrow(in,dati,nc,inr,kout,dato,nr,inc,*100)
      close(in,status=idisp)
        if (infil .ne. iofil) goto 86
        open(iout,file=iofil,status='unknown',form='unformatted')
c       now copy grid from scratch unit to iofil
        rewind (kout)
        read(kout)id,pgm,nc,nr,nz,xo,dz,yo,dy
        write(iout)id,pgm,nc,nr,nz,xo,dx,yo,dy
        do 85 k=1,nr
        call rdbin(kout,dati,1,dato,nc,*100)
        call iorow(iout,dati,1,dato,nc)
85      continue
        close(kout)
86      close(iout)
      return
  100 return 1
      end
      subroutine splln(data,n1,igap,i1pos,i2pos,i1st,i2st,a,ipdv)
c
c   this subroutine calculates spline coefficients (a) for gap between
c    s-n or w-e data sets.  igap = gap in grid units between non-dval
c    points of left and right of gap.
c
      dimension data(n1),s(10),sx(10),a(4,10)
      data odval/1.0e38/
      ispos=7
      istrt=i1st
      i1=i1pos
   18 ispos=ispos-1
      s(ispos)=igap
      sx(ispos-1)=data(i1)
      if(ispos.eq.2) go to 27
      if(istrt.ge.1) go to 25
      go to 27
   25 do 26 i1=istrt,1,-1
      if(data(i1).lt.odval) go to 30
   26 continue
c
c   here we have .le. 5 pts on left of boundary for splining.
c
   27 nptlf=7-ispos
      i1pt=ispos-1
      s(i1pt)=0.0
      go to 31
   30 igap=istrt-i1+1
      istrt=i1-1
      go to 18
c
c   here there are 5 points to spline on left of boundary.  now get 5 on
c    right of boundary.
c
   31 ispos=7
      sx(6)=data(i2pos)
      if(i2st.gt.n1) go to 34
      istrt=i2st
   29 do 32 i2=istrt,n1
      if(data(i2).lt.odval) go to 33
   32 continue
      go to 34
c
c   at 33, we have good spline value on right of gap.
c
   33 igap=i2-istrt+1
      s(ispos)=igap
      sx(ispos)=data(i2)
      ispos=ispos+1
      istrt=i2+1
      if(ispos.le.10.and.i2.lt.n1) go to 29
c
c   here there are .le. 5 points on right side of boundary for splining.
c
   34 nptrt=ispos-6
      npts=nptlf+nptrt
      if(npts.le.2) go to 51
c
c   now spline the 10 values.
c
      call splnon(sx(i1pt),s(i1pt),a(1,i1pt),npts)
c     iend=i1pt+npts-1
c     write(6,50) i1pt,npts,(sx(i),s(i),(a(j,i),j=1,4),i=i1pt,iend)
c  50 format(2x,'spline coeffs:',/,2x,2i10,/,(2x,6f12.2))
      return
   51 ipdv=1
      return
      end
       subroutine splnon(g,s,a,im)
c      this program finds the spline coefficients for irregularly-
c      spaced data
       dimension g(50),s(50),xc(50),xd(50),yx(50),a(4,50)
       ixx=im-1
       st=s(2)*s(2)
       xc(1)=0.5
       xd(1)=1.5*(g(2)-g(1))/s(2)
       do 1 il=2,ixx
       dnm=2.0*(s(il)+s(il+1))-s(il+1)*xc(il-1)
       xcnt=1.0/dnm
       xc(il)=s(il)*xcnt
       st2=s(il)*s(il)
       st=s(il+1)*s(il+1)
       axx=3.0/(s(il)*s(il+1))
    1  xd(il)=xcnt*(axx*(st*(g(il)-g(il-1))+st2*(g(il+1)-g(il)))-
     1 s(il+1)*xd(il-1))
       xs=(3.0*(g(im)-g(ixx))/s(im))-xd(ixx)
       xt=2.0-xc(ixx)
       yx(im)=xs/xt
       k=im
       do 2 j=2,im
       k=k-1
    2  yx(k)=xd(k)-xc(k)*yx(k+1)
       do 3 j=1,ixx
       a(1,j)=g(j)
       a(2,j)=yx(j)
       st=g(j+1)-g(j)-yx(j)*s(j+1)
       st2=(yx(j+1)-yx(j))/(2.0*s(j+1))
       axx=s(j+1)*s(j+1)
       a(4,j)=2.0*(st2-st/axx)/s(j+1)
    3  a(3,j)=(st/axx)-a(4,j)*s(j+1)
       return
       end
      subroutine rowspl(ist,iend,ispl,nr,ir,spl,rows[huge])
c   this subroutine performs spline interpolation between columns ist
c    and iend to get interpolated values for row ispl.  spl array
c    contains ir rows needed for the interpolation.  rows array contains
c    irow rows .ge. ispl row for later outputting.  nr = no. of rows
c    in overall grid.
c
      dimension spl(14,100),a(4,10),rows(1000,8)
      data odval/1.0e38/
      irpos=ist
      il=min0(7,ispl-1)
      iu=min0(7,nr-ispl)
      ifin=iend-ist+1
      do 188 i=1,ifin
c
c   now set up spline for i'th column.  note that at least one point
c    must exist on each side of row ispl.
c
      do 186 j=il,1,-1
      if(spl(j,i).lt.odval) go to 177
  186 continue
c
c   here only dvals found to south of boundary.  pad with dval.
c
      go to 173
  177 i1pos=j
      i1st=j-1
      ibeg=il+1
      do 187 j=ibeg,ir
      if(spl(j,i).lt.odval) go to 178
  187 continue
c
c   here all dvals found to north of boundary.
c
      go to 173
  178 i2pos=j
      i2st=j+1
      igap=i2pos-i1pos+1
      ipdv=0
      call splln(spl(1,i),ir,igap,i1pos,i2pos,i1st,i2st,a,ipdv)
      if(ipdv.ne.0) go to 173
c
c   now add splined values to row ispl.
c
      is=il-i1pos+1
      val=a(4,5)
      do 175 j=3,1,-1
  175 val=val*is+a(j,5)
      rows(irpos,1)=val
  173 irpos=irpos+1
c
c   note that irpos always points to the next posn in the current row
c    that requires a splined value.
c
  188 continue
      return
      end
      subroutine rdrow(in,data[huge],nc,inr,iout,outd[huge],nr,inc,*)
      dimension data(nc,inr),outd(nr,inc),id(14),pgm(2),dum(7)
        dimension zero(1)
      data zero/0.0/
c     write(6,4) inr,inc
c   4 format(2x,'inr = ',i4,', inc = ',i4)
      itotc=inc
      iec=inc
      numc=0
    1 numr=0
      idc=numc*inc
      ier=inr
      itotr=inr
    2 do 5 ir=1,ier
      call rdbin(in,dum,1,data(1,ir),nc,*50)
    5 continue
c
c   now output data.
c
      idr=numr*inr
      do 10 i=1,iec
      iposc=idc+i
      do 10 j=1,ier
      iposr=idr+j
   10 outd(iposr,i)=data(iposc,j)
c     write(6,11) itotr
c  11 format(2x,'itotr = ',i4)
      if(itotr.ge.nr) go to 20
      ier=nr-itotr
      ier=min0(ier,inr)
      itotr=itotr+ier
      numr=numr+1
      go to 2
c
c   output iec columns of data.
c
   20 do 25 i=1,iec
      call iorow(iout,zero,1,outd(1,i),nr)
   25 continue
c     write(6,26) itotc
c  26 format(2x,'itotc = ',i4)
      if(itotc.ge.nc) go to 30
      iec=nc-itotc
      iec=min0(iec,inc)
      itotc=itotc+iec
      numc=numc+1
      rewind in
      read(in,end=30) id,pgm,dum
      go to 1
   30 return
   50 return 1
      end

