      subroutine spjbfg(nx,ny,x,f,fgrad,task,ecc,b)
      character*(*) task
      integer nx, ny
      real f, ecc, b
      real x(nx*ny), fgrad(nx*ny)
c     **********
c
c     Subroutine spjbfg
c
c     This subroutine computes the function and gradient of the
c     pressure distribution in a journal bearing problem.
c
c     The subroutine statement is
c
c       subroutine spjbfg(nx,ny,x,f,fgrad,task,ecc,b)
c
c     where
c
c       nx is an integer variable.
c         On entry nx is the number of grid points in the first
c            coordinate direction.
c         On exit nx is unchanged.
c
c       ny is an integer variable.
c         On entry ny is the number of grid points in the second
c            coordinate direction.
c         On exit ny is unchanged.
c
c       x is a real array of dimension nx*ny.
c         On entry x specifies the vector x if task = 'F', 'G', or 'FG'.
c            Otherwise x need not be specified.
c         On exit x is unchanged if task = 'F', 'G', or 'FG'. Otherwise
c            x is set according to task.
c
c       f is a real variable.
c         On entry f need not be specified.
c         On exit f is set to the function evaluated at x if task = 'F'
c            or 'FG'.
c
c       fgrad is a real array of dimension nx*ny.
c         On entry fgrad need not be specified.
c         On exit fgrad contains the gradient evaluated at x if
c            task = 'G' or 'FG'.
c
c       task is a character variable.
c         On entry task specifies the action of the subroutine:
c
c            task               action
c            ----               ------
c             'F'     Evaluate the function at x.
c             'G'     Evaluate the gradient at x.
c             'FG'    Evaluate the function and the gradient at x.
c             'XS'    Set x to the standard starting point xs.
c             'XL'    Set x to the lower bound xl.
c
c         On exit task is unchanged.
c
c       ecc is a real variable
c         On entry ecc is the eccentricity in (0,1).
c         On exit ecc is unchanged
c
c       b is a real variable
c         On entry b defines the domain as D = (0,2*pi) X (0,2*b).
c         On exit b is unchanged.
c
c     MINPACK-2 Project. November 1993.
c     Argonne National Laboratory and University of Minnesota.
c     Brett M. Averick and Jorge J. More'.
c
c     **********
      real four, one, p5, six, two, zero
      parameter (zero=0.0,p5=0.5,one=1.0,two=2.0,four=4.0,six=6.0)

      logical feval, geval
      integer i, j, k
      real dvdx, dvdy, ehxhy, flin, fquad, hx, hxhy, hy, pi, temp,
     +     trule, v, vb, vl, vr, vt, xi

      real p

      p(xi) = (1+ecc*cos(xi))**3

c     Initialization.

      pi = four*atan(one)
      hx = two*pi/real(nx+1)
      hy = two*b/real(ny+1)
      hxhy = hx*hy
      ehxhy = ecc*hxhy

c     Compute the lower bound xl for x if task = 'XL'.

      if (task .eq. 'XL') then
         do 10 k = 1, nx*ny
            x(k) = zero
   10    continue

         return

      end if

c     Compute the standard starting point if task = 'XS'.

      if (task .eq. 'XS') then
         do 30 i = 1, nx
            temp = max(sin(real(i)*hx),zero)
            do 20 j = 1, ny
               k = nx*(j-1) + i
               x(k) = temp
   20       continue
   30    continue

         return

      end if

      if (task .eq. 'F' .or. task .eq. 'FG') then
         feval = .true.
      else
         feval = .false.
      end if
      if (task .eq. 'G' .or. task .eq. 'FG') then
         geval = .true.
      else
         geval = .false.
      end if

c     Compute the function if task = 'F', the gradient if task = 'G', or
c     both if task = 'FG'.

      if (feval) then
         fquad = zero
         flin = zero
      end if
      if (geval) then
         do 40 k = 1, nx*ny
            fgrad(k) = zero
   40    continue
      end if

c     Computation of the quadratic part of the function and
c     corresponding components of the gradient over the
c     lower triangular elements.

      do 60 i = 0, nx
         xi = real(i)*hx
         trule = hxhy*(p(xi)+p(xi+hx)+p(xi))/six
         do 50 j = 0, ny
            k = nx*(j-1) + i
            v = zero
            vr = zero
            vt = zero
            if (i .ne. 0 .and. j .ne. 0) v = x(k)
            if (i .ne. nx .and. j .ne. 0) vr = x(k+1)
            if (i .ne. 0 .and. j .ne. ny) vt = x(k+nx)
            dvdx = (vr-v)/hx
            dvdy = (vt-v)/hy
            if (feval) fquad = fquad + trule*(dvdx**2+dvdy**2)
            if (geval) then
               if (i .ne. 0 .and. j .ne. 0) fgrad(k) = fgrad(k) -
     +             trule*(dvdx/hx+dvdy/hy)
               if (i .ne. nx .and. j .ne. 0) fgrad(k+1) = fgrad(k+1) +
     +             trule*dvdx/hx
               if (i .ne. 0 .and. j .ne. ny) fgrad(k+nx) = fgrad(k+nx) +
     +             trule*dvdy/hy
            end if
   50    continue
   60 continue

c     Computation of the quadratic part of the function and
c     corresponding components of the gradient over the upper
c     triangular elements.

      do 80 i = 1, nx + 1
         xi = real(i)*hx
         trule = hxhy*(p(xi)+p(xi-hx)+p(xi))/six
         do 70 j = 1, ny + 1
            k = nx*(j-1) + i
            vb = zero
            vl = zero
            v = zero
            if (i .ne. nx+1 .and. j .ne. 1) vb = x(k-nx)
            if (i .ne. 1 .and. j .ne. ny+1) vl = x(k-1)
            if (i .ne. nx+1 .and. j .ne. ny+1) v = x(k)
            dvdx = (v-vl)/hx
            dvdy = (v-vb)/hy
            if (feval) fquad = fquad + trule*(dvdx**2+dvdy**2)
            if (geval) then
               if (i .le. nx .and. j .gt. 1) fgrad(k-nx) = fgrad(k-nx) -
     +             trule*dvdy/hy
               if (i .gt. 1 .and. j .le. ny) fgrad(k-1) = fgrad(k-1) -
     +             trule*dvdx/hx
               if (i .le. nx .and. j .le. ny) fgrad(k) = fgrad(k) +
     +             trule*(dvdx/hx+dvdy/hy)
            end if
   70    continue
   80 continue

c     Computation of the linear part of the function and
c     corresponding components of the gradient.

      do 110 i = 1, nx
         temp = sin(real(i)*hx)
         if (feval) then
            do 90 j = 1, ny
               k = nx*(j-1) + i
               flin = flin + temp*x(k)
   90       continue
         end if
         if (geval) then
            do 100 j = 1, ny
               k = nx*(j-1) + i
               fgrad(k) = fgrad(k) - ehxhy*temp
  100       continue
         end if
  110 continue

c     Finish off the function.

      if (feval) f = p5*fquad - ehxhy*flin

      end
