1      subroutine fpsuev(idim,tu,nu,tv,nv,c,u,mu,v,mv,f,wu,wv,lu,lv)
2c  ..scalar arguments..
3      integer idim,nu,nv,mu,mv
4c  ..array arguments..
5      integer lu(mu),lv(mv)
6      real*8 tu(nu),tv(nv),c((nu-4)*(nv-4)*idim),u(mu),v(mv),
7     * f(mu*mv*idim),wu(mu,4),wv(mv,4)
8c  ..local scalars..
9      integer i,i1,j,j1,k,l,l1,l2,l3,m,nuv,nu4,nv4
10      real*8 arg,sp,tb,te
11c  ..local arrays..
12      real*8 h(4)
13c  ..subroutine references..
14c    fpbspl
15c  ..
16      nu4 = nu-4
17      tb = tu(4)
18      te = tu(nu4+1)
19      l = 4
20      l1 = l+1
21      do 40 i=1,mu
22        arg = u(i)
23        if(arg.lt.tb) arg = tb
24        if(arg.gt.te) arg = te
25  10    if(arg.lt.tu(l1) .or. l.eq.nu4) go to 20
26        l = l1
27        l1 = l+1
28        go to 10
29  20    call fpbspl(tu,nu,3,arg,l,h)
30        lu(i) = l-4
31        do 30 j=1,4
32          wu(i,j) = h(j)
33  30    continue
34  40  continue
35      nv4 = nv-4
36      tb = tv(4)
37      te = tv(nv4+1)
38      l = 4
39      l1 = l+1
40      do 80 i=1,mv
41        arg = v(i)
42        if(arg.lt.tb) arg = tb
43        if(arg.gt.te) arg = te
44  50    if(arg.lt.tv(l1) .or. l.eq.nv4) go to 60
45        l = l1
46        l1 = l+1
47        go to 50
48  60    call fpbspl(tv,nv,3,arg,l,h)
49        lv(i) = l-4
50        do 70 j=1,4
51          wv(i,j) = h(j)
52  70    continue
53  80  continue
54      m = 0
55      nuv = nu4*nv4
56      do 140 k=1,idim
57        l3 = (k-1)*nuv
58        do 130 i=1,mu
59          l = lu(i)*nv4+l3
60          do 90 i1=1,4
61            h(i1) = wu(i,i1)
62  90      continue
63          do 120 j=1,mv
64            l1 = l+lv(j)
65            sp = 0.
66            do 110 i1=1,4
67              l2 = l1
68              do 100 j1=1,4
69                l2 = l2+1
70                sp = sp+c(l2)*h(i1)*wv(j,j1)
71 100          continue
72              l1 = l1+nv4
73 110        continue
74            m = m+1
75            f(m) = sp
76 120      continue
77 130    continue
78 140  continue
79      return
80      end
81