      subroutine polinb(cob,inb,jnb,ww,wfbc,rank,nxin,nyin,im,jm,nm)
!c-----------------------------------------------------------------------c
!c       routine for bilinear interpolation of fields                    c
!c                                                                       c
!c      --> cob    relative position of b-grid points in the avn grid    c
!c                 box                                                   c
!c      -->inb,jnb  indices of the nearest four points to the b-grid     c
!c                  point                                                c
!c      -->nxin,nyin  dimension of the input dataset                     c
!c      -->ww       input initial dataset                                c
!c      <--wfbc     dataset interpolated to b-grid points                c
!c      -->rank     type of interpolations                               c
!c                  2: wind                                              c
!!      <-- im,jm,nm output dimension                                    c
!c                                                                       c
!c-----------------------------------------------------------------------c 

      implicit none

      integer::nxin,nyin,im,jm,nm,i,j,n,rank
      real,dimension(3,0:im+1,0:jm+1,nm):: cob
      integer,dimension(4,0:im+1,0:jm+1,nm)::inb,jnb
      real,dimension(nxin,nyin):: ww
      real,dimension(0:im+1,0:jm+1,nm)::wfbc
      integer::rand,imt,jmt,i00,i10,i01,i11,j00,j01,j10,j11,ii,jj
      real::pq,p,q,dy,x1a(4),x2a(4),ya(4,4)

      x1a=(/-1.,0.,1.,2./)
      x2a=(/-1.,0.,1.,2./)
      
      if(rank.eq.2)then
         imt=im-1
         jmt=jm-1
      else
         imt=im
         jmt=jm
      endif
!
      do n=1,nm
      do i=1,imt
      do j=1,jmt
        i00=inb(1,i,j,n)
        i10=inb(2,i,j,n)
        i01=inb(3,i,j,n)
        i11=inb(4,i,j,n)
        j00=jnb(1,i,j,n)
        j10=jnb(2,i,j,n)
        j01=jnb(3,i,j,n)
        j11=jnb(4,i,j,n)

        p=cob(2,i,j,n)
        q=cob(3,i,j,n)
        pq=cob(1,i,j,n)

        if(j00.eq.1.or.j00.eq.nyin-1)then
      wfbc(i,j,n)=ww(i00,j00)+p*(ww(i10,j10)-ww(i00,j00))  &
                          +q*(ww(i01,j01)-ww(i00,j00))     &
        +pq*(ww(i00,j00)-ww(i10,j10)-ww(i01,j01)+ww(i11,j11))
	else
	  if(i00.eq.1)then
	    ya(1,:)=ww(nxin,j00-1:j00+2)
	    ya(2:4,:)=ww(1:3,j00-1:j00+2)
          else if(i00.eq.nxin-1)then
	    ya(1:3,:)=ww(nxin-2:nxin,j00-1:j00+2)
	    ya(4,:)=ww(1,j00-1:j00+2)
          else
	    ya=ww(i00-1:i00+2,j00-1:j00+2)
          endif
	  call polin2(x1a,x2a,ya,4,4,p,q,wfbc(i,j,n),dy)
        endif
	  
       enddo
       enddo
       enddo
!-----------------------------------------------------------------------
      return
      end


