cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
cc                                                                    cc
cc                  mnspgr : spgrid test program                      cc
cc                                                                    cc
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c  ..local scalars..
      real del,ermax,erf,exr0,exr1,fp,pi,sum,r0,r1,one,ai,s
      integer i,ier,is,j,k,kwrk,lwrk,m,mu,mv,nc,nuest,nu,nvest,nv
c  ..local arrays..
      integer ider(4),iopt(3),iwrk(70),iw(25)
      real u(11),v(14),r(154),c(300),tu(25),tv(25),f(154),wk(100),
     * exact(154),err(14),sp(14),wrk(1500)
c  ..function references..
      real abs,atan2,tesspg
c  ..
c  set constants
      one = 1
      pi = atan2(0.,-one)
      del = pi*0.05
c we set up the number of u (latitude)-values of the grid.
      mu = 11
c we set up the u-coordinates of the grid.
      read(5,905)(iw(i),i=1,mu)
      do 10 i=1,mu
         ai = iw(i)
         u(i) = ai*del
  10  continue
c we set up the number of v (longitude)-values of the grid
      mv = 14
c we set up the v-coordinates of the grid.
      read(5,905)(iw(i),i=1,mv)
      do 20 i=1,mv
         ai = iw(i)
         v(i) = ai*del
  20  continue
c we fetch the data values at the grid points.
      m = mu*mv
      read(5,900) (r(i),i=1,m)
c we print the data values at the grid points. we also compute and print
c the exact value of the test function underlying the data.
      write(6,910)
      write(6,915) (j,j=1,mv,2)
      write(6,920)
      exr0 = tesspg(0.,0.)
      exr1 = tesspg(pi,0.)
      ermax = 0.
      sum = 0.
      l = 0
      do 40 i=1,mu
        l = (i-1)*mv+1
        k = 1
        do 30 j=1,7
          exact(l) = tesspg(u(i),v(k))
          erf = abs(exact(l)-r(l))
          sum = sum+erf
          if(erf.gt.ermax) ermax = erf
          sp(j) = exact(l)
          err(j) = r(l)
          l = l+2
          k = k+2
  30    continue
        write(6,925) i,(err(j),j=1,7)
        write(6,930) (sp(j),j=1,7)
  40  continue
      write(6,915) (j,j=2,mv,2)
      write(6,920)
      do 60 i=1,mu
        l = (i-1)*mv+2
        k = 2
        do 50 j=1,7
          exact(l) = tesspg(u(i),v(k))
          erf = abs(exact(l)-r(l))
          sum = sum+erf
          if(erf.gt.ermax) ermax = erf
          sp(j) = exact(l)
          err(j) = r(l)
          l = l+2
          k = k+2
  50    continue
        write(6,925) i,(err(j),j=1,7)
        write(6,930) (sp(j),j=1,7)
  60  continue
      ai = m
      sum = sum/ai
      write(6,935) sum,ermax
      write(6,940) exr0,exr1
c  we set up the dimension information
      nuest = 19
      nvest = 21
      kwrk = 70
      lwrk = 1500
c main loop for the different spline approximations
      do 300 is=1,6
        go to (110,120,130,140,150,160),is
c  we start computing a set of spline approximations with
c  only c0-continuity at the poles,
 110    iopt(2) = 0
        ider(2) = 0
        iopt(3) = 0
        ider(4) = 0
c  with no data values at the poles.
        ider(1) = -1
        ider(3) = -1
c  initialisation
        iopt(1) = 0
c  a large value for s for computing the least-squares polynomial
        s = 60.
        go to 200
c  iopt(1) = 1 from the second call on
 120    s = 0.05
        iopt(1) = 1
        go to 200
c  an interpolating spline
 130    s = 0.
        go to 200
c  a second set of approximations with c1-continuity at the poles
 140    iopt(2) = 1
        iopt(3) = 1
c  exact values at the poles.
        ider(1) = 1
        ider(3) = 1
        r0 = exr0
        r1 = exr1
c reinitialization
        iopt(1) = 0
        s = 0.05
        go to 200
c  vanishing derivatives at the poles
 150    ider(2) = 1
        ider(4) = 1
c reinitialization
        iopt(1) = 0
        go to 200
c finally we calculate the least-squares spline according to the current
c  set of knots
 160    iopt(1) = -1
 200    call spgrid(iopt,ider,mu,u,mv,v,r,r0,r1,s,nuest,nvest,
     *    nu,tu,nv,tv,c,fp,wrk,lwrk,iwrk,kwrk,ier)
c printing of the fitting results.
        if(iopt(1).ge.0) go to 210
        write(6,945)
        go to 220
 210    write(6,950) s
 220    write(6,955) iopt(2),iopt(3)
        if(ider(2).eq.1) write(6,960)
        if(ider(4).eq.1) write(6,965)
        write(6,970) fp,ier
        write(6,975) nu
        write(6,980)
        write(6,985) (tu(i),i=1,nu)
        write(6,990) nv
        write(6,980)
        write(6,985) (tv(i),i=1,nv)
        nc = (nu-4)*(nv-4)
        write(6,995)
        write(6,985) (c(i),i=1,nc)
c  evaluation of the spline approximation
        call bispev(tu,nu,tv,nv,c,3,3,u,mu,v,mv,f,wk,100,iw,25,ier)
        write(6,1000)
        write(6,915) (j,j=1,mv,2)
        write(6,920)
        ermax = 0.
        sum = 0.
        do 240 i=1,mu
          k = i
          do 230 j=1,mv
            sp(j) = f(k)
            err(j) = abs(exact(k)-f(k))
            sum = sum+err(j)
            if(err(j).gt.ermax) ermax = err(j)
            k = k+1
 230      continue
          if( (i/2)*2 .ne.i ) go to 240
          write(6,925) i,(sp(j),j=1,mv,2)
          write(6,930) (err(j),j=1,mv,2)
 240    continue
        sum = sum/ai
        write(6,935) sum,ermax
        write(6,1005) c(1),c(nc)
 300  continue
      stop
 900  format(7f8.3)
 905  format(14i3)
 910  format(49h1data value (exact function value) at grid points)
 915  format(8h0v(j),j=,3x,7(i2,6x))
 920  format(8h u(i),i=)
 925  format(1x,i5,7(2x,f6.3))
 930  format(7x,7(2h (,f5.3,1h)))
 935  format(19h0mean abs. error = ,f9.3,5x,18hmax. abs. error = ,f9.3)
 940  format(30h function values at the poles ,f7.3,5x,f7.3)
 945  format(21h0least-squares spline)
 950  format(25h0smoothing spline with s=,f7.2)
 955  format(1x,34horder of continuity at the poles =,2i5)
 960  format(1x,37hvanishing derivatives at the pole u=0)
 965  format(1x,38hvanishing derivatives at the pole u=pi)
 970  format(1x,23hsum squared residuals =,e15.6,5x,11herror flag=,i3)
 975  format(1x,42htotal number of knots in the u-direction =,i3)
 980  format(1x,22hposition of the knots )
 985  format(5x,8f9.4)
 990  format(1x,42htotal number of knots in the v-direction =,i3)
 995  format(23h0b-spline coefficients )
1000  format(50h0spline value (approximation error) at grid points)
1005  format(28h spline values at the poles ,f7.3,5x,f7.3)
      end
c
      real function tesspg(u,v)
c function program tesspg calculates the value of the test function
c underlying the data.
      real u,v,sin,cos
      tesspg = 2./(4.1+cos(3.*u)+3.*cos(v+v+u*0.25)*sin(u)**2)
      return
      end