      SUBROUTINE TERPV ( DEGLAT, DEGLON, FINT, LEN, F )
      include 'processors.h'
 
C  BILINEAR INTERPOLATION OF A FIELD, F, TO A SET OF LAT/LON LOCATIONS
C  DEGLAT(I) IS ASSUMED TO BE BETWEEN  -90.0 AND  +90.0
C  DEGLON(I) IS ASSUMED TO BE BETWEEN -180.0 AND +180.0


C**********************************************************************
C  Data parallel implementation (CM Fortran) : M. Makivic (NPAC)      *
C  Date:       4/94                                                   *
C**********************************************************************

 
      real,dimension(len):: deglat,yj,deglon,xi,fint,wtj1,wti1
      real,dimension(len):: a,b,c,d,d1,d2,d3,d4
      integer,dimension(len):: ii,i0,i1,j0,j1
      integer,dimension(len):: ifl
      real,dimension(im,jnp):: f
      real              fill
      REAL              QDEGDX
      REAL              QDEGDY
!hpf$ distribute deglat(block) onto PROC_ARRAY_1D
!hpf$ distribute deglon(block) onto PROC_ARRAY_1D
!hpf$ distribute fint(block) onto PROC_ARRAY_1D
!hpf$ distribute f(*,block) onto PROC_ARRAY_1D

!hpf$ align with deglat:: yj,xi,wtj1,wti1
!hpf$ align with deglat:: a,b,c,d,d1,d2,d3,d4
!hpf$ align with deglat:: ii,i0,i1,j0,j1,ifl
 
C  IM AND JNP ARE NEEDED FROM COMMON BLOCK /CONTRL/
      LOGICAL           STRAT,   SLP,     FDCOS,   FEXP,    SAVEAE,
     $                  SIGDAT
      COMMON /CONTRL/   IM,      JNP,     NB,      NDALT,   NHMS,
     $                  NLAY,    NLAYM1,  NLAYP1,  NYMD,    STRAT,
     $                  IDIMEN,  SLP,     FDCOS,   FEXP,    SAVEAE,
     $                  SIGDAT,  MAXREP
 
      data fill /1.e10/
 
      QDEGDX = IM / 360.0
      QDEGDY = (JNP-1) / 180.0
 
 
      YJ    = ( DEGLAT +  90.0 ) * QDEGDY + 1.0
      J0    = YJ
      J0    = MAX0( J0, 1 )
      J0    = MIN0( J0, JNP-1 )
      J1    = J0 + 1

      XI    = ( DEGLON + 180.0 ) * QDEGDX + 1.0
      WHERE ( XI.GE.(IM+1) ) 
           XI = XI - IM
      end where
      WHERE( XI.LT.1 )
           XI = XI + IM
      ENDWHERE
      I0    = XI
      I0    = MAX0( I0, 1 )
      I0    = MIN0( I0, IM )
      I1    = I0 + 1
      WHERE( I0.EQ.IM ) 
          I1 = 1
      ENDWHERE

      WTJ1  = YJ - J0
      WTI1  = XI - I0

c  Screen for fill values on level - replace fill values according 
c  to predetermined criteria. 

C
C For brevity, the four corner points will be identified as follows:
C
C         X X   8               X X   4
C         X X                   X X
C         C  (D3)               D  (D4)
C
C
C
C
C         X X   1               X X   2
C         X X                   X X
C         A   (D1)              B   (D2)
C
C
C
C Assign corner values
C

c      FORALL(L=1:LEN) A(L)=F(I0(L),J0(L))
c      FORALL(L=1:LEN) B(L)=F(I1(L),J0(L))
c      FORALL(L=1:LEN) C(L)=F(I0(L),J1(L))
c      FORALL(L=1:LEN) D(L)=F(I1(L),J1(L))
 
      do L=1,LEN 
          A(L)=F(I0(L),J0(L))
          B(L)=F(I1(L),J0(L))
          C(L)=F(I0(L),J1(L))
          D(L)=F(I1(L),J1(L))
      end do
C Identify pattern of filled points by using a 4 bit binary representation,
C e.g., 2 (decimal) = 0010 (binary) = corner B only is filled. 
 
      IFL  = 0
 
      WHERE( A .EQ. FILL)
              IFL = 1
      ENDWHERE
      WHERE ( B .EQ. FILL) 
              IFL = IFL + 2
      ENDWHERE
      WHERE( D .EQ. FILL) 
              IFL = IFL + 4
      ENDWHERE
      WHERE( C .EQ. FILL)
              IFL = IFL + 8
      ENDWHERE
 
C
C Filled values are repalced with values obtained from the other corners
C so that the horizontal interpolation can be done.  Note that there are 
C 16 special cases of filled patterns. 
C
      WHERE( IFL .EQ. 0 )
              D1 = A
              D2 = B
              D3 = C
              D4 = D
      end where
      WHERE( IFL .EQ. 1 )
              D1 = .5 * ( B + C )
              D2 = B
              D3 = C
              D4 = D
      end where
      WHERE( IFL .EQ. 2 )
              D1 = A
              D2 = .5 * ( A + D )
              D3 = C
              D4 = D
      end where
      WHERE( IFL .EQ. 3 )
              D1 = .5 * ( C + D )
              D2 = D1
              D3 = C
              D4 = D
      end where
      WHERE( IFL .EQ. 4 )
              D1 = A
              D2 = B
              D3 = C
              D4 = .5 * ( B + C )
      end where
      WHERE( IFL .EQ. 5 )
              D1 = .5 * ( B + C )
              D2 = B
              D3 = C
              D4 = D1
      end where
      WHERE( IFL .EQ. 6 )
              D1 = A
              D2 = .5 * ( A + C )
              D3 = C
              D4 = D2
      end where
      WHERE( IFL .EQ. 7 )
              D1 = C
              D2 = C
              D3 = C
              D4 = C
      end where
      WHERE( IFL .EQ. 8 )
              D1 = A
              D2 = B
              D3 = .5 * ( A + D )
              D4 = D
      end where
      WHERE( IFL .EQ. 9 )
              D1 = .5 * ( B + D )
              D2 = B
              D3 = D1
              D4 = D
      end where
      WHERE( IFL .EQ. 10)
              D1 = A
              D2 = .5 * ( A + D )
              D3 = D2
              D4 = D
      end where
      WHERE( IFL .EQ. 11)
              D1 = D
              D2 = D
              D3 = D
              D4 = D
      end where
      WHERE( IFL .EQ. 12)
              D1 = A
              D2 = B
              D3 = .5 * ( A + B )
              D4 = D3
      end where
      WHERE( IFL .EQ. 13)
              D1 = B
              D2 = B
              D3 = B
              D4 = B
      end where
      WHERE( IFL .EQ. 14)
              D1 = A
              D2 = A
              D3 = A
              D4 = A
      end where
      WHERE( IFL .EQ. 15)
              D1 = FILL
              D2 = FILL
              D3 = FILL
              D4 = FILL
      ENDWHERE
 
C  HORIZONTAL INTERPOLATION
 
      WHERE (D1.EQ.FILL) 
         FINT =  FILL 
      ELSEWHERE 
         A = D1 
         B = D2 
         C = D3 
         D = D4 
         FINT =  WTJ1  * ( WTI1*D+(1.0-WTI1)*C )
     $        + (1.0-WTJ1) * ( WTI1*B+(1.0-WTI1)*A )
      ENDWHERE
 
 
      RETURN
      END
