      subroutine bill(wina, woua,faclama, facphia,
     *                winb, woub,faclamb, facphib, 
     *                  rlatitude,rlongitude, ja, jb)

      include "slgcm.slgprm"
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 bilinear scheme to interpolate the fields wina and winb to the point
c  (rlatitude,rlongitude) and puts the interpolated value in woua woub.
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  ja and jb are the beggining and ending latitude indexes
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
       real   wina( imp3 ,jmp3 )   !first input field
       real   woua( imp3 ,jm   )   !first output field
       real   winb( imp3 ,jmp3 )   !second input field
       real   woub( imp3 ,jm   )   !second output field

       real     rlatitude  ( imp3 ,jm   ) ! arrays with coordinates on
       real     rlongitude ( imp3 ,jm   ) ! which the subroutine interpolates
c					  ! the fields wina and winb
       integer dk, dj


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c Bilinear interpolation function:
c
       bli(delx,dely,w1,w2,w3,w4) = (1.-delx)*((1.-dely)*w1+dely*w2) +
     >                                  delx *((1.-dely)*w3+dely*w4)
c                                |                         |
c                       ---------+-------------------------+----------
c                                |                         |
c                                |                         |
c                                |                         |
c                                |                         |
c                                |                         |
c                           ^    |        X                |
c                           |    |rlatitudeij,rlongitudeij |
c                           |    |                         |
c                     dely  |    |                         |
c                           |    |                         |
c                       ----x----+-------------------------+-----------
c                                |<------->                |          
c                                   delx
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
C1              print *, '+WINA'
C1              do dj = 1, jmp3
C1                 do dk = 1, imp3
C1                    print 101, wina(dk,dj)
C1                 end do
C1              end do
C1              print 101, ((wina(dk, dj), dj = 1, jmp3), dk = 1, imp3)
C1              print *, '+WINA OVER'
C1              print *, '+WINB'
C1              do dj = 1, jmp3
C1                 do dk = 1, imp3
C1                    print 101, winb(dk,dj)
C1                 end do
C1              end do
C1              print 101, ((winb(dk, dj), dj = 1, jmp3), dk = 1, imp3)
C1              print *, '+WINB OVER'
C1              print *, '+RLAT'
C1              do dj = 1, jm
C1                 do dk = 1, imp3
C1                    print 102, rlatitude(dk,dj)
C1                 end do
C1              end do
C1              print 102, ((rlatitude(dk, dj), dj = 1, jm), dk = 1, imp3)
C1              print *, '+RLAT OVER'
C1              print *, '+RLON'
C1              do dj = 1, jm
C1                 do dk = 1, imp3
C1                    print 102, rlongitude(dk,dj)
C1                 end do
C1              end do
C1              print 102, ((rlongitude(dk, dj), dj = 1, jm), 
C1     &             dk = 1, imp3)
C1              print *, '+RLON OVER'
C1              print 103, faclama, faclamb, facphia, facphib
C1 101   format ( F16.8 )
C1 102   format ( F16.8 )
C1
C 101  format( 1X, 8(4F16.8 //), 3F16.8 //)
C 102  format( 1X, 8(4F16.8 //))
 103  format( 1X, '+', 'Faclama : ', F16.8, ' Faclamb : ', 
     &         F16.8, ' Facphia : ', F16.8, ' Facphib : ', F16.8 //)
      
c
c r transforms radians into "turns" 
	r=one/twopi
c dlm is 2*pi/im
c u and v quantities are offset differently from the equator.
        facdlma = (2.-faclama)*dlm 
        facdlmb = (2.-faclamb)*dlm
c main loop over all grid point in the horizontal range 
	do j=ja,jb    
	 do i=1,imp3

c add enough turns as to have a positive lon
c and compute the index number of the grid point to the
c east of the lon 
          tem1 = (  rlongitude(i,j)-facdlma)*r   
          ii =  tem1
C1          PRINT *, 'First time ii', ii
c if started with a negative give a turn more to get into positive 
c angles
	 if(tem1.lt.0.) ii=ii-1
C1          PRINT *, 'Second Time ii', ii
c get the actual gridpoint index in the lon direction
	   tem1=(tem1-ii)*im+two
           ii=tem1
C1          PRINT *, 'Third Time ii', ii
           tem1=tem1-ii
           if(ii.eq.1)ii=imp1
C1          PRINT *, 'Third Time ii', ii
c same thing for the latitudes 
c the computation for core1 can be used for the second field
         core1= (  rlatitude(i,j) + pio2)*rdph 
         tem2 = core1+ facphia
         jj = tem2 
         tem2= tem2   - jj

        woua(i,j) =  bli ( tem1 , tem2,
     >                wina(ii  ,jj+2  ), wina (ii    ,jj+3),
     >                wina(ii+1,jj+2  ), wina (ii+1  ,jj+3) )

c
c
c
c
c the comments are the same as above
          tem1 = (  rlongitude(i,j)-facdlmb)*r  
          ii =  tem1
	  if(tem1.lt.0.) ii=ii-1
	  tem1=(tem1-ii)*im+two
          ii=tem1
          tem1=tem1-ii
          if(ii.eq.1)ii=imp1
c same thing for the latitudes 
          tem2= core1+ facphib
          jj = tem2 
          tem2= tem2   - jj

        woub(i,j) =  bli ( tem1 , tem2,
     >                winb(ii  ,jj+2  ), winb (ii    ,jj+3),
     >                winb(ii+1,jj+2  ), winb (ii+1  ,jj+3) )


         enddo 
         enddo 
C2         stop

	return
	end
        
        
