      subroutine bic1nu (wina, woua,faclama,facphia,
     *                   winb, woub,faclamb,facphib,
     *                            rlatitude, rlongitude,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
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  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
c
c
      include "slgcm.slgprm"
c
c wina and win b are the two fields on a regular grid
c
      real   wina(imp3,jmp3)
      real   winb(imp3,jmp3)
c
c Rlatitude and rlongitude are the coordinates of the points
c to which the fields wina and winb are to be interpolated
c
      real   rlatitude(imp3,jm  )
      real   rlongitude(imp3,jm  )
c
c woua and woub are the interpolated values
c
c
      real   woua(imp3,jm  )
      real   woub(imp3,jm  )
CD
      real core1
      integer icount
      common /serial/ icount
      integer ny, ny0, nysth
      integer lo1, up1, width
      integer m1jj(imp3,jmp3), m2jj(imp3, jmp3), max_iter(2)
CD
c
c facdlma facdlmv are the distance in radias from the dateline to the
c first row of points in the longitude direction, dlm is 2*pi/im
c
      facdlma = (2.-faclama)*dlm
      facdlmb = (2.-faclamb)*dlm

c
c
c the quntities ijb and ije are still needed because this subroutine 
c wraps the fields and the wrap subroutine is still in one dimension.
c in the future wrapping should be taken out of here and put into slgcst2d
c for clarity. j0 
c
C1      print *, '!= start of jj'
C1      print *, '!= start of dely'
C      print *, '!= start of ii'

      r=one/twopi
c the procedure is the same as in the routine bill
      do j=ja, jb
      do i=1,imp3
CD
CD          if ((j .gt. 1) .and. (j .le. 4)) then
CD            print '(" ITER : J ", I4, " I ", I4)', j, i
CD            print '(" rlon", F16.8, 
CD     &       " facdlma", F4.2,
CD     &       " r", F16.8)', rlongitude(i,j), facdlma, r
CD         endif
CD
           temx = (rlongitude(i,j)-facdlma)*r   
           ii =  temx
CD
CD           if ((j .gt. 1) .and. (j .le. 4)) print '(" tem1", F16.8, 
CD     &     " II", I4)',temx, ii
CD
           if(temx.lt.0.) ii=ii-1
CD
CD          if ((j .gt. 1) .and. (j .le. 4)) PRINT *, 'Second Time ii ', 
CD     &     ii
CD
	   temx=(temx-ii)*im+two
           ii=temx
CD
CD          if ((j .gt. 1) .and. (j .le. 4)) print '(" tem1", F16.8, 
CD     &     " II ", I4)', temx, ii
CD

           temx=temx-ii
           if(ii.eq.1)ii=imp1
CD
CD          if ((j .gt. 1) .and. (j .le. 4)) print '(" tem1", F16.8, 
CD     &     " II ", I4)', temx, ii
CD
CD

           temy   = (rlatitude(i,j) + pio2)*rdph + facphia
             jj   = temy
CD       print *, '!=', jj
CD 101     format(E14.4E2)
CD 109     format(I4,":", I4, ":",F16.8,":", F16.8)
CD         print 109, ii, jj, temx, temy
         
           temy   = temy - jj
C1         print 109, ii, jj, temx, temy
CD        m1jj(i,j) = jj
      	call  cubit( temx,temy,
     .                wina(ii  ,jj+2  ),wina(ii,jj+3  ),
     .                wina(ii  ,jj+1  ),wina(ii,jj+4  ),
     .                wina(ii+1,jj+2  ),wina(ii+1,jj+3),
     .                wina(ii+1,jj+1  ),wina(ii+1,jj+4),
     .                wina(ii-1,jj+2  ),wina(ii-1,jj+3),
     .                wina(ii-1,jj+1  ),wina(ii-1,jj+4),
     .                wina(ii+2,jj+2  ),wina(ii+2,jj+3),
     .                wina(ii+2,jj+1  ),wina(ii+2,jj+4), woua(i,j))



CD
CD          if ((j .gt. 1) .and. (j .le. 4)) then
CD            print '(" ITER : J ", I4, " I ", I4)', j, i
CD            print '(" rlon", F16.8, 
CD     &       " facdlmb", F4.2,
CD     &       " r", F16.8)', rlongitude(i,j), facdlmb, r
CD         endif
CD
           temx = (rlongitude(i,j)-facdlmb)*r   
           ii =  temx
CD
CD           if ((j .gt. 1) .and. (j .le. 4)) print '(" tem1", F16.8, 
CD     &     " II", I4)',temx, ii
CD
	   if(temx.lt.0.) ii=ii-1
CD
CD          if ((j .gt. 1) .and. (j .le. 4)) PRINT *, 'Second Time ii ', 
CD     &     ii
CD
	   temx=(temx-ii)*im+two
           ii=temx
CD
CD          if ((j .gt. 1) .and. (j .le. 4)) print '(" tem1", F16.8, 
CD     &     " II ", I4)', temx, ii
CD
           temx=temx-ii
           if(ii.eq.1)ii=imp1
CD
CD          if ((j .gt. 1) .and. (j .le. 4)) print '(" tem1", F16.8, 
CD     &     " II ", I4)', temx, ii
CD
CD           if ((j .gt. 1) .and. (j .le. 4)) then
CD              print '(" rlat", F16.8, " facphib", F16.8)',
CD     &        rlatitude(i,j), facphib
CD              print '(" rdph", F16.8, " pio2", F16.8)', rdph, pio2
CD            core1 = (rlatitude(i,j) + pio2) * rdph
CD            print '(" core1", F16.8)', core1
CD         endif

CD
           temy   = (rlatitude(i,j) + pio2)*rdph + facphib
             jj   = temy
CD
CD          if ((j .gt. 1) .and. (j .le. 4)) print '(" tem2", F16.8, 
CD     &     " JJ ", I4)', temy, jj
CD
           temy   = temy - jj
CD
CD          if ((j .gt. 1) .and. (j .le. 4)) print '(" tem2", F16.8)',
CD     &    temy
CD
CD        m2jj(i,j) = jj
      	call  cubit( temx,temy,
     .                winb(ii  ,jj+2  ),winb(ii,jj+3  ),
     .                winb(ii  ,jj+1  ),winb(ii,jj+4  ),
     .                winb(ii+1,jj+2  ),winb(ii+1,jj+3),
     .                winb(ii+1,jj+1  ),winb(ii+1,jj+4),
     .                winb(ii-1,jj+2  ),winb(ii-1,jj+3),
     .                winb(ii-1,jj+1  ),winb(ii-1,jj+4),
     .                winb(ii+2,jj+2  ),winb(ii+2,jj+3),
     .                winb(ii+2,jj+1  ),winb(ii+2,jj+4), woub(i,j))

       enddo
       enddo

CD       call get_lsd_bc_rng(width, ja, jb, m1jj, max_iter,
CD     &                     ny, ny0, nysth)       
CD       call get_lsd_bc_rng(width, ja, jb, m2jj, max_iter,
CD     &                     ny, ny0, nysth)       
 100   format ( 4(F16.10))
 107   format ( I4, I4, F16.8, F16.8)
C1          open (unit = 10, file = 'dwina.dat', status = 'new' )
C1          write (10, 100) ((wina(i,j), i = 1, imp3), j = 1, jmp3)
C1          close(10)
C1          open (unit = 11, file = 'dwinb.dat', status = 'new' )
C1          write (11, 100) ((winb(i,j), i = 1, imp3), j = 1, jmp3)
C1          close(11)
C1          open (unit = 12, file = 'drlat.dat', status = 'new' )
C1          write (12, 100) ((rlatitude(i,j), i = 1, imp3), j = 1, jm)
C1          close(12)
C1          open (unit = 13, file = 'drlon.dat', status = 'new' )
C1          write (13, 100) ((rlongitude(i,j), i = 1, imp3), j = 1, jm)
C1          close(13)
C1          
CD          open (unit = 14, file = 'woua.dat', status = 'new' )
CD          do j = 1, jm
CD              do k = 1, imp3
CD                 write (14,107), k, j, woua(k,j), woub(k,j)
CD              end do
CD           end do
CD          close(14)
C1          open (unit = 15, file = 'woub.dat', status = 'new' )
C1           do j = 1, jm
C1              do k = 1, imp3
C1                 write (15,107), k, j, woub(k,j)
C1              end do
C1           end do
C1          close(15)


C1        print *, '!= end of jj'
C1          print *, '!= end of dely'
C1          print *, '!= end of ii'

      call wrap(ja,jb,woua)
      call wrap(ja,jb,woub)
CD
      icount = icount + 1
 109  format ( 4F18.13 )

C1      open (unit = 14, file = 'woua64_n.out.ser', status = 'scratch' )
C1      close(14)
C1      open (unit = 14, file = 'woua64_n.out.ser', status = 'new' )
C1      write (14,*), "ICOUNT : ", icount
C1      write (14,*), "WOUA"
C1      write (14,109), ((woua(k,j), k = 1, imp3), j = 1, jm)
C1      close(14)
C1      if (icount .eq. 1) then
C1      do j = 1, jm
C1         do k = 1, imp3
C1            print *, woua(k,j)
C1         enddo
C1      enddo
C1      endif
C1      print *, "Writing results 2"
C1      open (unit = 15, file = 'woub64_n.out.ser', status = 'scratch' )
C1      close(15)
C1      open (unit = 15, file = 'woub64_n.out.ser', status = 'new' )
C1      write (15,*), "ICOUNT : ", icount
C1      write (15,*), "WOUB"
C1      write (15, 109), ((woub(k,j), k = 1, imp3), j = 1, jm)
C1      close (15)
C1
C1      if (icount .eq. 4) then
C1         open (unit = 13, file = 'wou64_4.out.ser', status = 'scratch' )
C1         close(13)
C1         open (unit = 13, file = 'wou64_4.out.ser', status = 'new' )
C1         write(13,*), "WOUA"
C1         write (13, 109), ((woua(k,j), k = 1, imp3), j = 1, jm)
C1         write(13,*), "WOUB"
C1         write (13, 109), ((woub(k,j), k = 1, imp3), j = 1, jm)
C1         close (13)
C1      endif

      return
      end

