!   &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
       SUBROUTINE PUTVEG(HLAT,HLON,  &
                         FVEG0,SM,SICE,FVEG1,im,jm,nm)
!   **************************************************************
!   * Interpolate NESDIS vegetation fraction product (five-year  *
!   * climatology with 0.144 degree resolution from 89.928S,180W * 
!   * to 89.928N, 180E)                                          *
!   * F. Chen 07/96                                              *
!   **************************************************************

!
      PARAMETER  (L0=2500*1250)
!
      integer::nk,im,jm,nm
      DIMENSION FVEG0(2500,1250), FVEG1(0:im+1,0:jm+1,nm)  
      DIMENSION HLAT(0:im+1,0:jm+1,nm),HLON(0:im+1,0:jm+1,nm)
      DIMENSION SM(0:im+1,0:jm+1,nm),SICE(0:im+1,0:jm+1,nm)
      real::dlmd,dphd,x1,y1,z1,x2,y2,z2
      integer::i1,j1,i2,j2,k
!      INTEGER IPOPT(20), JPDS(25), JGDS(22), KGDS0(22), KGDS1(22)
!     +      , KPDS0(25), KPDS1(25), BITFLG1(1)
!      CHARACTER GDS1(42)
!      LOGICAL BIT0(2500,1250), BIT1(0:im+1,0:jm+1,nm)
!FEI      DATA IG/96/

      if(nm.eq.6)then
        nk=1
      else if(nm.eq.14)then
        nk=3
      endif

      i1=(im+1)/2
      j1=(jm+1)/2
      i2=i1+1
      j2=j1+1

      x1=cos(hlat(i1,j1,nk))*cos(hlon(i1,j1,nk))
      y1=cos(hlat(i1,j1,nk))*sin(hlon(i1,j1,nk))
      z1=sin(hlat(i1,j1,nk))
      x2=cos(hlat(i2,j2,nk))*cos(hlon(i2,j2,nk))
      y2=cos(hlat(i2,j2,nk))*sin(hlon(i2,j2,nk))
      z2=sin(hlat(i2,j2,nk))              

      dlmd=sqrt((x1-x2)**2+(y1-y2)**2+(z1-z2)**2)*180/3.1415926
      dphd=dlmd

	dtr=acos(-1.)/180.
        hlat=hlat/dtr
        hlon=hlon/dtr

!   This subroutine AVERAGES the value stored in FVEG0 TO the FVEG1 array.
!   Special effort is made to ensure that islands are taken care of when
!   not resolved on the 0.144 x 0.144 deg. grid.
!  *** Define the search limit for isolated island point 
        NLIM=30*(1+AINT(0.15/DPHD))
!  *** Define the Eta averaging box size, NBOX=0 gives nearest neighbor
        NBOX1=NINT(MAX(DLMD,DPHD)/0.144)
        NBOX=NBOX1
!        print*,' 0:im+1,0:jm+1=',0:im+1,0:jm+1,'NBOX=',NBOX,'DPHD',DPHD,'NLIM',NLIM,
!     &         'DLMD',DLMD

        do k=1,nm
        DO J = 1,JM
        DO I = 1,IM

!mp 	added for workstation
	HLON(i,j,k) = 360.0 - HLON(i,j,k)
	IF(HLON(i,j,k) .GT. 360.) HLON(i,j,k) = HLON(i,j,k) - 360.
!mp

          IF((SM(i,j,k).GT.0.5).OR.(SICE(i,j,k).EQ.1.0)) THEN
            FVEG1(i,j,k) = 0.0
          ELSE
            DX=(HLON(i,j,k) - 179.928)/0.144
            INDX = NINT(DX)
!  ***  Here, 179.928 is the starting longitude (180.072) + one grid cell 
!       width (0.144) so that the first index is one NOT 0.
            IF(INDX.LT.1) THEN 
              DX=(HLON(i,j,k) + 180.072)/0.144
              INDX = NINT(DX)
            ENDIF
            DY=(HLAT(i,j,k) + 90.072)/0.144
            INDY = NINT(DY)
!  *** Get area-average value of FVEG1 from input grid FVEG0, for 
!      finer Eta grid, this becomes nearest-neighbour 
            ICON1=0
  100       CONTINUE
            IF(MOD(NBOX,2).EQ.0) THEN
              IF(DX.GT.REAL(INDX)) THEN
                NLO=MAX(1,INDX-NBOX/2+1)
                NHI=MIN(2500,INDX+NBOX/2)
              ELSE
                NLO=MAX(1,INDX-NBOX/2)
                NHI=MIN(2500,INDX+NBOX/2-1)
              ENDIF
              IF(DY.GT.REAL(INDY)) THEN
                MLO=MAX(1,INDY-NBOX/2+1)
                MHI=MIN(1250,INDY+NBOX/2) 
              ELSE
                MLO=MAX(1,INDY-NBOX/2)
                MHI=MIN(1250,INDY+NBOX/2-1)
              ENDIF
            ELSE
              NLO=MAX(1,INDX-NBOX/2)
              NHI=MIN(2500,INDX+NBOX/2)
              MLO=MAX(1,INDY-NBOX/2)
              MHI=MIN(1250,INDY+NBOX/2)
            ENDIF
!            NLO=MAX(1,INDX-NBOX)
!            NHI=MIN(2500,INDX+NBOX)
!            MLO=MAX(1,INDY-NBOX)
!            MHI=MIN(1250,INDY+NBOX)
            ICON=0
            VEGSUM=0.0
            DO N=NLO, NHI
              DO M=MLO, MHI
                IF(FVEG0(N,M).NE.0.0) THEN
                 ICON=ICON+1 
                 VEGSUM=VEGSUM+FVEG0(N,M)
                ENDIF
              ENDDO
            ENDDO
            IF(ICON.NE.0) THEN
              FVEG1(i,j,k)=VEGSUM/ICON
              NBOX=NBOX1
            ELSE
! *** Search for an isolated point 
              ICON1=ICON1+1
              NBOX=NBOX+1 
              IF(ICON1.LE.NLIM) THEN
                GO TO 100
              ELSE
                FVEG1(i,j,k)=0.55
                NBOX=NBOX1
!                print *,"OH OH ",i,j,k,HLAT(i,j,k),HLON(i,j,k)
              ENDIF
            ENDIF
! ** End of land case
!          if(NBOX.NE.NBOX1) print*,' NBOX is NOT NBOX1,=',NBOX
          END IF  
          IF(FVEG1(i,j,k).LT.0.0.OR.FVEG1(i,j,k).GT.1.0) THEN 
           print*,' FVEG1 out of range, FVEG=',FVEG1(i,j,k)
          ENDIF
! ** End of im,jm,nm loops
        END DO  
        END DO  
        END DO  
        RETURN
        END