      subroutine bicpnu(vst, win, wou, rlatitude, rlongitude, 
     *    dtoa, faclam,facphi,ja,jb)

c
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c                       #     # ######   #####   #####
c                       #     # #     # #     # #     #
c                       #     # #     # #       #
c                       ####### ######  #       #
c                       #     # #       #       #
c                       #     # #       #     # #     #
c                       #     # #        #####   #####
c                       
c                          #######  #####   #####
c                          #       #     # #     #
c                          #       #       #
c                          #####    #####   #####
c                          #             #       #
c                          #       #     # #     #
c                          #######  #####   #####
c                       
c                          ######     #     #####
c                          #     #   # #   #     #
c                          #     #  #   #  #     #
c                          #     # #     # #     #
c                          #     # ####### #     #
c                          #     # #     # #     #
c                          ######  #     #  #####
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c	Jose Zero NASA GSFC DAO -- zero@talisman.gsfc.nasa.gov
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c Purpose:     
c  Uses a bicubic scheme to interpolate the field win to the point
c  (rlatitude,rlongitude) and puts the interpolated value in wout.
c  on a c-grid, faclam and facphi must be set as follows depending on
c  whether the points being used are the u v or fi points.
c  Note that this routine acts on a single field because its call 
c  in slgcst2d.f is inmediately preceded by a call to one of the 
c zangs routines that probably changest the values of the latitude and 
c longitude 
c  u points: faclam = 0.5; facphi = 1.0
c  v points: faclam = 1.0; facphi = 0.5
c  h points: faclam = 1.0; facphi = 1.0
      include "slgcm.slgprm"
      include "slgcm.slgclm"
c
      real       win(imp3,jmp3)
      real       wou(imp3,jm  )
c
      real   rlatitude (imp3,jm)
      real   rlongitude(imp3,jm)
      real          vst(imp3,jm)
      real          cvv(2      )
      real        glm00(2      )

c
      facdlm = (2.-faclam)*dlm
      r=one/twopi
c
c I feel that the values at the poles can be computed outside of this
c subroutine and then it would not be neccessay to have a special 
c subroutine for this .... may be ... the answer will be in a future 
c version of the model.
c
c
c     calculate pole values
c
      call vpole (im, dlm, twopi, twooim, pi
     *,       vst, clm(1,2), slm(1,2), cvv(1), glm00(1))

      call vpole (im, dlm, twopi, twooim, zero
     *,      vst(1,jmm1), clm(1,2), slm(1,2), cvv(2), glm00(2))

      do j1x = 1, 2
         cvv(j1x) = cvv(j1x)*dtoa
      end do

      call fipole(win, wou, cvv   , glm00   , .true.)

      call fipole(win, wou, cvv(2), glm00(2), .false.)
c
c
c  main loop over all the points in the grid
c

      do j=ja, jb
      do i=1,imp3
c
c  first make sure dlm < lam* < (twopi+dlm) for the v and fi points
c  and that 1.5*dlm < lam* < (twopi+1.5*dlm) for the u points
c
c      the procedure is identical to the one in bill add enough "turns
c      so that the longitude in the desired range find the pointers
c      compute the deltas and interpolate. the bicubi interpolation is
c      identical to the one in bic1nu so the same subroutine "cubit" is used
c compute number of turns
           temx = (rlongitude(i,j)-facdlm)*r   
           ii =  temx
	   if(temx.lt.0.) ii=ii-1
c compute the deltax 
	   temx=(temx-ii)*im+two
           ii=temx
           temx=temx-ii
c in case it is the first gridpoint make is its "wrapped up" so we can
c name "ii-1" easily
           if(ii.eq.1)ii=imp1
c compute delta y and its index.
           temy   = (rlatitude(i,j) + pio2)*rdph + facphi
             jj   = temy
           temy   = temy - jj
c interpolate
      	call  cubit( temx,temy,
     .                win(ii  ,jj+2  ),win(ii,jj+3  ),
     .                win(ii  ,jj+1  ),win(ii,jj+4  ),
     .                win(ii+1,jj+2  ),win(ii+1,jj+3),
     .                win(ii+1,jj+1  ),win(ii+1,jj+4),
     .                win(ii-1,jj+2  ),win(ii-1,jj+3),
     .                win(ii-1,jj+1  ),win(ii-1,jj+4),
     .                win(ii+2,jj+2  ),win(ii+2,jj+3),
     .                win(ii+2,jj+1  ),win(ii+2,jj+4), wou(i,j))
        enddo
        enddo


      return
      end
