      complex function zquint (z,m,n,x,y,zdx,zdy)
cc      subroutine zquintc (z,m,n,x,y,zquint,zdx,zdy)
      complex*8  z,zquint_f,zdx,zdy
      integer    m,n
      real       x,y
      dimension  z(m,n)
      save

      i17=17
cc      call zquintc_f (z,i17,m,n,x,y,zquint_f,zdx,zdy)
cc      zquint=zquint_f
      zquint=zquint_f (z,i17,m,n,x,y,zdx,zdy)
      return
      end

      complex function zquint_f (z,flz,m,n,x,y,zdx,zdy)
cc      subroutine zquintc_f (z,flz,m,n,x,y,zquint_f,zdx,zdy)
C Interpolation * 2-D Quadratic polyn., full array.
c
c And another interpolation routine.
c This time the complex array z(m,n) has data everywhere.
c Returned value is interpolated quadratically at 1 < x < m, 1 < y < n.
c Additionally returned: The gradient at (x,y), (zdx,zdy)
c
c UNIX-version GEO mode:
c zdx and zdy are in physical units; length unit=rad. 
c If this is not desired, call Quint_Opt ('+G-P')
c
c UNIX-version P/S calls:
c Reading a global tide model should set common area /crsar/.
c /crsar/ parameters can alternately be set using 
c     Call Set_Got_Ar (...)
c c.f. oteu54.f If /crsar/ thus describes the geography of Z(M,N), 
c     Call Quint_Opt ('+G') 
c signifies GEO mode. Zquint will interpret X,Y as longitude,latitude. 
c Circular connection is determined by the GLOBAL parameter in common /crsar/.
c
c N.B. in GEO mode Z(.,1) is the northernmost row of data. 
c
c N.B. there is no automatic /crsar/ validisation; default values fit the
c Schwiderski grid. 
c
c In non-GEO mode, 
c     Call Quint_Opt ('+C') 
c will declare Z(M,N) as circularly connected in the first array subscript. 
c Default option is '-C-G'.
c
c To process Catesian arrays with north row first 
c     Call Quint_Opt ('-G+R')
c
c To process spherical grids with the south row first:
c    Call Quint_Opt ('+G+R') 
c
c To enter Trace mode,
c     Call Quint_Opt ('+T') 
c
c More in the header of Subroutine Quint_Qpt below
c
      complex*8 z(m,n),zr(6),zp(6),zv,zdx,zdy
     &,         zt,ztx,zty,zvv,zvx,zvy
      logical   qi,qgeo,qphys,qcc,qg,qt,global,qrv,qdp
     &,         qflg,qall,qignore_flags,qsterp,qgric
      integer   flz(m,n)
      character type*1
      dimension ik(6),jk(6), a(6,6), iia(6),jja(6)
      common /cquint/ re,qgeo,qphys,qcc,qi,qt,qrv,qignore_flags
     &,               qsterp,qgric,type
      save   /cquint/, qg, a, qi

      data ik/0,1,2,0,0,1/, jk/0,0,0,1,2,1/
     &,     a/6*1.0
     &,       0.,1.,2.,0.,0.,1.
     &,       0.,1.,4.,0.,0.,1.
     &,       0.,0.,0.,1.,2.,1.
     &,       0.,0.,0.,1.,4.,1.
     &,       0.,0.,0.,0.,0.,1./
      data qg/.true./

cccc      qall=.true.
cccc      goto 1
cccc      
cccc      entry zquint_f (z,flz,m,n,x,y,zdx,zdy)
      qall=.false.
      if (flz(1,1).eq.17) qall=.true.

 1    if (qi) then
         if (qg) then
            call gaussj (a,6,6,dummy,0,1)
         endif
         resx=1.
         resy=1.
         if (qgeo) then
            call enq_got_ar (xw,yn,xe,ys,resx,resy,global)
            call enq_got_rnd (xrnd,yrnd)
            qcc=global
            if (qt) then
               print *,'<zquint>t> after enq_got_ar:'
               print *,'<zquint>t> ',xw,yn,xe,ys,resx,resy,global
               print *,'<zquint>t> ',xrnd,yrnd
            endif
         elseif (qsterp) then
            call enqsca (scale,qdp)
            if (type.eq.'Z') then
               hx=float(m+1)/2.
               hy=float(n+1)/2.
            else
               hx=float(m)/2.
               hy=float(n)/2.
            endif
         else
            cylat=1.
         endif
         rad=atan(1.)/45.
         qi=.false.
         qg=.false.
         if (qt) print *,'<zquint>>> init done'
      endif
      if (qt) print *,a
      if (qgeo) then
         yj=(yn-y)*resy-yrnd+1.
         xi=(x-xw)*resx-xrnd
         cylat=cos(y*rad)
      elseif (qsterp.and..not.qgric) then
         call sterp (x,y,xkm,ykm)
         xi=xkm/scale+hx-1.
         yj=ykm/scale+hy-1.
      else
         yj=y-1.
         xi=x-1.
      endif
      jl=max0(1,min0(n-2,ifix(yj)))
      if (qcc) then
         xi=amod(xi+float(m),float(m))
         il=ifix(xi)
      else
         if (xi.gt.m.or.yj.lt.-.5.or.yj.gt.n+.5) then
            if (qt) print *,'<zquint>>> Outside: xi,yj=',xi,yj
            zquint_f=(0.,0.)
            zdx=(0.0,0.0)
            zdy=(0.0,0.0)
            return
         endif
         il=max0(1,min0(m-2,ifix(xi)))
      endif
      dy=yj-jl
      dx=xi-il
      m1=m-1
      jjl=jl
      if (qrv) jjl=n+1-jl
      if (qt) print *,'<zquint>>> loc at '
     &,         x,y,il,jl,jjl,xi,yj,dx,dy,z(il,jjl)

      nav=6
      qflg=qall.or.qignore_flags
      do 101 k=1,6
      if (qcc) then
         ii=mod(il+m1+ik(k),m)+1
      else
         ii=il+ik(k)
      endif
      jj=jl+jk(k)
      iia(k)=ii
      jja(k)=jj
      if (qrv) jj=n+1-jj
      if (.not.qall) then
         qflg=flz(ii,jj).ne.0
      endif
      if (qflg) then
         zr(k)=z(ii,jj)
      else
         zr(k)=0.0d0
         nav=nav-1
      endif
 101  continue
      
      if (nav.le.0) then
c this only happens if there is a flag file involved
         zquint_f=(0.0,0.0)
         zdx=(0.0,0.0)
         zdy=(0.0,0.0)
         return
      endif

      if (nav.ne.6) then
c this only happens if there is a flag file involved
         zv=(0.0,0.0)
         vv=0.0
         do k=1,6
            vzr=abs(zr(k))
            if (vv.le.1.e-10.and.vzr.gt.1.e-10) then
               vv=vzr
               zdx=zr(k)
            endif
            zv=zv+zr(k)
         enddo
         zquint_f=zv/nav
         zdy=(0.0,0.0)
         return
      endif
         
      if (qt) then
         print *,'<zquint>>> Grid values: ' 
         print '(2i5,1p,2e16.8)', (iia(k),jja(k),zr(k),k=1,6)
      endif
      do 1022 i=1,6
      zp(i)=(0.0,0.0)
      do 102 j=1,6
 102     zp(i)=zp(i)+a(i,j)*zr(j)
 1022 continue
      
      zv=zp(1)
      ztx=zp(3)*dx*resx
      ztx=ztx+zp(2)
      ztx=ztx*dx*resx
      zty=zp(5)*dy*resy
      zty=zty+zp(4)
      zty=zty*dy*resy
      zvv=zp(6)*(dx*dy)*resx*resy
      zv=zv+ztx
      zv=zv+zty
      zv=zv+zvv
cc      if (qall) then
cc         zquint=zv
cc      else
         zquint_f=zv
cc      endif
      ztx=zp(3)*dx*resx
      zty=zp(5)*dy*resy
      zvx=zp(6)*dx*resx
      zvy=zp(6)*dy*resy
      zdx=zp(2)*resx
      zdx=zdx+2.*ztx
      zdx=zdx+zvy
      zdy=zp(4)*resy
      zdy=zdy+2.*zty
      zdy=zdy+zvx
      if (qt) print '(a,2f9.4,2x,1p,3(2e12.4,2x))'
     &,  ' <zquint>>> RAW:   ',resx,resy,zv,zdx,zdy
c      if (qt) print *,'<zquint>p> ',zp 
      if (qphys) then
         if (qsterp) then
            zdx=-zdx*re/scale
            zdy=-zdy*re/scale
            if (qt) print '(a,1p,3(2e12.4,2x))'
     &,             ' <zquint>>> PHY/S: ',zv,zdx,zdy 
         else
cerr            zdx=zdx*resx/cylat/rad/re
cerr            zdy=-zdy*resy/rad/re
            zdx=zdx/cylat/rad/re
            zdy=-zdy/rad/re
            if (qt) print '(a,0p,f10.4,a,14x,1p,2(2e12.4,2x))'
     &,             ' <zquint>>> PHY/G @',cylat,': ',zdx,zdy 
         endif
      endif   

      return
      end

      subroutine quint_opt (opt)
c
c Options: The system wakes up with -G -S -P -C -R -T, ER=6.371e6 [m]
c
c +C      - Assume array is circular connected.
c
c +I      - Ignore flags, eventhough zquint_f is called
c
c -G      - Assume plane grid, Cartesian relation between array row and y
c           coordinate. That means the first array row is south of the rest.
c +G      - Assume circular connection according to "/crsar/ ...,global".
c           Assume northernmost row has number 1, and
c           use area bounds according to "/crsar/..." 
c           Imply +P (qphys=.true.), i.e. gradients are per radian.
c
c           To set common /crsar/ use  Set_Got_Ar (...)  or read a tide
c           model with a global reading routine. 
c
c +S:Z or +S:M  Stereographic grid. Specify grid kind Z or M
c +SGZ or +SGM  Stereographic grid, calling with coordinates on the grid
c               (ploating point, 0<=X<=M, 0<=Y<=N)
c               You can code "=" or "g" for "G".
c
c +R      - Reverse order north/south. Obs! Question: What's normal ?
c           With -G, +R means that the first row of the array is the north
c           of the others. With +G, +R means that the first row of the  
c           array is south of the others. 
c
c -P      - Gradients are per unit step in the Cartesian grid. 
c +P      - Physical units of the gradients ON THE EARTH SURFACE. 
c +Punit  - The following units can be given:
c            +P/Mm    - gradient per 1000 km
c            +P/km    - gradient per km
c            +P/m     - gradient per m
c            +P/deg   - gradient per degree
c
c ER=val, - Numeric value of the planet radius in m.
c
c +T      - Switch on trace option = extensive print.
c
      character opt*(*), typ*1
      logical quch_in, q_index

      common /cquint/ re,qgeo,qphys,qcc,qi,qt,qrv,qignore_flags
     &,               qsterp,qgric,typ
      save   /cquint/
      logical qgeo, qphys, qcc, qi, qt, qrv, qsterp
     &,       qignore_flags, qgric

      data    qgeo/.false./, qphys/.false./, qcc/.false./, 
     &        qi/.true./, qt/.false./, qrv/.false./, qsterp/.false./
     &,       qignore_flags/.false./, qgric/.false./

      data    re/1.0/, er/6.371e6/, typ/' '/

      entry zquint_opt (opt)
      entry rquint_opt (opt)

      if (quch_in('-S',opt,i)) then
         qsterp=.false.
         qgric=.false.
      endif
      if (quch_in('+S',opt,i)) then
         qsterp=.true.
         qgric=index('Gg=',opt(i+2:i+2)).gt.0
         typ=opt(i+3:i+3)
      endif
      if (quch_in('+I',opt,i)) qignore_flags=.true.
      if (quch_in('-I',opt,i)) qignore_flags=.false.
      if (quch_in('+C',opt,i)) qcc=.true.
      if (quch_in('-C',opt,i)) qcc=.false.
      if (quch_in('+G',opt,i)) qgeo=.true.
      if (quch_in('-G',opt,i)) qgeo=.false.
      if (quch_in('+G',opt,i)) qphys=.true.
      if (quch_in('-G',opt,i)) qphys=.false.
      if (quch_in('+T',opt,i)) qt=.true.
      if (quch_in('-T',opt,i)) qt=.false.

      if (quch_in('+P',opt,i)) qphys=.true.
      if (q_index(opt,'ER=',i)) read (opt(i+3:),*) er
      if (q_index(opt,'+P1000km',i).or.
     &    q_index(opt,'+PMm',i).or.
     &    q_index(opt,'+P/Mm',i).or.
     &    q_index(opt,'+P/1000km',i)) then
         re=er*1e-6
      elseif (q_index(opt,'+Pkm',i).or.
     &        q_index(opt,'+P/km',i)) then
         re=er*1e-3
      elseif (q_index(opt,'+Pm',i).or.
     &        q_index(opt,'+P/m',i)) then
         re=er
      elseif (q_index(opt,'+Pdeg',i).or.
     &        q_index(opt,'+P/deg',i)) then
         re=45./atan(1.)
      else
         re=1.
      endif

      if (quch_in('-P',opt,i)) qphys=.false.
      if (quch_in('-R',opt,i)) qrv=.false.
      if (quch_in('+R',opt,i)) qrv=.true.
      qi=.true.
      print '(a,5l2,1p,2e12.4)'
     &,' <QuintO>>> CirC, RevO, Geo, Ster, Phys, R_unit, ER ='
     &,      qcc,qrv,qgeo,qsterp,qphys,re,er
      if (qsterp) print '(a,l2,2a)'
     &,' <QuintO>>> STEREOGRAPHIC mode: coordinates on grid =',qgric 
     &,', TYPE=',typ
      if (qt) print '(a)',' <QuintO>>> TRACE MODE'
      return
      end

C:Provides  subroutine quint_opt F=zquint.f
C:Provides  complex function zquint F=zquint.f
C:Provides  entry zquint_f F=zquint.f
C:Provides  entry zquint_opt F=zquint.f
C:Provides  entry rquint_opt F=zquint.f
