      SUBROUTINE KMVEGT(TLM0D,TPH0D,DLMD,DPHD
     +               ,WBD,SBD,SM,sice,vegmap)
      INCLUDE "parmeta"
      PARAMETER (IMT=IM*2-1,JMT=JM/2+1)
      PARAMETER  (IMJM=IM*JM-JM/2)
      PARAMETER  (IMA=7200,JMA=9000,IMJMA=IMA*JMA)
      PARAMETER  (NTYP = 12)
      DIMENSION NSIBV(NTYP),NEGRD(IMJM,NTYP),NMAX(IMJM)
      DIMENSION NSIBG(6,3) ,NEGRDG(IMJM,   6)
      DIMENSION SM(IM,JM),sice(IM,JM),SMK(IMJM)
      DIMENSION vegmap(im,jm)
      REAL      ALAT,ALON
      INTEGER   NTY(IMA,JMA)
C
      DATA NSIBV/1,2,3,4,5,6,7,8,9,10,11,12/
      DATA (NSIBG(1,I), I=1,3) / 1, 2,0/
      DATA (NSIBG(2,I), I=1,3) / 3, 4,0/
      DATA (NSIBG(3,I), I=1,3) / 5, 6,0/
      DATA (NSIBG(4,I), I=1,3) / 7, 8,0/
      DATA (NSIBG(5,I), I=1,3) / 9,10,0/
      DATA (NSIBG(6,I), I=1,3) /11,12,0/
C
c      DATA MISS/-99/,NZERO/0/
      DATA MISS/0/,NZERO/0/
C
C *** ALAT AND ALON ARE THE GEODETIC LATITUE 
C *** AND LONGITUDE OF THE 1KM LAT/LON GRID. 
C *** FOR EACH PAIR OF ALAT/ALON, FIND A K-INDEX ON ETA
C *** E GRID H POINT. THEN AT EVERY H POINT, 
C *** WE HAVE THE NUMBERS OF OB FOR THE 
C *** 20 VEG TYPES AND THE NUMBERS OF OB FOR THE 6 VEG 
C *** TYPE GROUPS. AT LAST, DETERMINE THE DOMINANT VEG 
C *** TYPE IN THE DOMAINT VEG GROUP ON EVERY H GRID POINT
C
       DO K = 1,IMJM
          NMAX (K)=MISS
       DO N = 1,NTYP
          NEGRD(K,N) =NZERO
          NEGRDG(K,N)=NZERO
       ENDDO
       ENDDO
      DO K=1,IMJM
        JX=(K-1)/IMT+1
        IX=K-(JX-1)*IMT
        IF(IX.LE.IM)THEN
          I=IX
          J=2*JX-1
        ELSE
          I=IX-IM
          J=2*JX
        ENDIF
        SMK(K)=SM(I,J)
      ENDDO
         write(65) SMK

C***
C     READ SEA-LAND MASK FOR 40KM ETA-GRID (1d)
C      READ(10) SM
C
C     READ MASK in 2-D WITH 1KM RESOLUTION 
C     READ VEGETATION TYPE MASK in 1-D WITH 1KM RESOLUTION 
C     (INTEGER VARIABLE)
C
	NNN=1
	NNOB=1
c
       dxy=0.008333333   !resolucao do mapa de veg. (1km)
c
CJLG        read(44,'(i2)')((NTY(i,j),i=1,ima),j=1,jma)
CJLG        write(63)NTY
        n=0
        do 210 nj=1,jma
            alat = -59.99583453 + (nj-1)*dxy
          do 210 ni=1,ima
              alon = -89.15898363 + (ni-1)*dxy
              n=n+1
           read(44,'(i2)')NTY(ni,nj)
C
c
	IF (NTY(ni,nj) .LE. 0 .OR. NTY(ni,nj) .GT. NTYP )then 
           NNN=NNN+1
	   GO TO 210
	ENDIF
C
C    BE WARE THE WEST IS POSITIVE.
C
C    The comment above is assuming the given longitudes
C     were positive. But since we have given the western
C     boundary as negative, we comment the line below 
C     in which the longitude was set to negative
C
         FLAT = ALAT
         FLON = -ALON
C
C        FIND KINDEX OF OBSERVATION
C
        ALN =   FLON
        APH =   FLAT
C
CJLG        CALL TOM( APH, ALN, KINDEX)
        CALL TOM(DLMD,DPHD,WBD,SBD,APH,ALN,KINDEX)
        IF(KINDEX .LE. 0 .OR. KINDEX .GT. IMJM) THEN
C            WRITE(6,6002) KINDEX, APH, ALN, ALAT
            NNOB=NNOB+1
            GOTO 210
        ENDIF
 6002  FORMAT(1x,'REJECTED:    KINDEX, APH, ALN, ALAT',I5,1x,3F11.6) 
       K= KINDEX
       IF(SMK(K) .gt. 0.5) GO TO 210
C
C     COUNT THE NUMBERS OF EVERY VEG TYPE IN ONE GRID BOX
C
        DO 500 NV = 1, NTYP
          IF(NTY(ni,nj) .eq. NSIBV(NV) ) THEN
            NEGRD(K,NV)=NEGRD(K,NV) + 1
        ENDIF
 500    CONTINUE
C
C    GROUP THE 20 VEG TYPES INTO 6 AND COUNT THE NUMBER OF EACH GROUP
C
      DO 620 NG  = 1, 6
      DO 600 NGG = 1, 3
         IF (NTY(ni,nj) .EQ. NSIBG(NG,NGG)) THEN
            NEGRDG(K,NG)= NEGRDG(K,NG)+1
         ENDIF
 600   CONTINUE
 620   CONTINUE
 210   CONTINUE
c
c  fim do primeiro loop   
c
       DO 220 K=1,IMJM
        IF(SMK(K) .EQ. 1.0) GO TO 220
	NTYY=0
C     DETERMINE THE DOMINANT GROUP IN THE GRID BOX(exclude water)
C
      MMAX = -99999
      DO 650 NG=1,6     
        IF (NEGRDG(K,NG) .LE.MMAX .OR. NEGRDG(K,NG) .EQ. NZERO) THEN
            GO TO 650
        ENDIF
	MMAX=NEGRDG(K,NG)
        NTYY=NG
 650  CONTINUE
      IF(NTYY .EQ. 0) THEN
C	  write(47,*) 'No available data at the grid point ',K
          GO TO 220
       ENDIF
C
C      NOW DETERMINE THE DOMINANT VEG TYPE IN THE DOMINANT GROUP
C
      MMAX = -99999
      NGG=1
       DO 302 NV = 1, NTYP
	 IF (NSIBV(NV) .NE. NSIBG(NTYY,NGG))  GO TO 302
	 IF (NGG .GT. 3) THEN
	    GO TO 220
	 ELSEIF (NEGRD(K,NV) .LT. MMAX .OR. NEGRD(K,NV) .EQ. NZERO) THEN
            NGG=NGG+1
	    GO TO 302
	 ELSE
	    NGGG=NGG
            NVV=NV
	    NMAX(K)=NSIBG(NTYY,NGGG)
	    NGG=NGG+1
	    MMAX=NEGRD(K,NV)
	    GO TO 302
	 ENDIF
  302  CONTINUE
c	write(48,*) 'K=',K,' NTYY= ',NTYY,' NVV= ',NVV,' NGG= ',NGGG
 220   CONTINUE
C
CJLG      WRITE(30) NMAX
          write(62) NMAX
C
c
c Convert veg_map from k-index into 2-dim index im,jm
c
      DO J=1,JM
      DO I=1,IM
        vegmap(I,J)=0.
      ENDDO
      ENDDO
C
      DO K=1,IMJM
        JX=(K-1)/IMT+1
        IX=K-(JX-1)*IMT
        IF(IX.LE.IM)THEN
          I=IX
          J=2*JX-1
        ELSE
          I=IX-IM
          J=2*JX
        ENDIF
        vegmap(I,J)=NMAX(K)
        IF(IX.EQ.IMT)vegmap(I+1,J)=vegmap(I,J)
      ENDDO
       write(66)vegmap
      WRITE(61,1061) (NMAX(K),K=1,IMJM)
 1061 FORMAT(1x,16(1x,I4))
C
      DO 720 NN = 1,IMJM
	if (NMAX(NN) .eq. MISS) GO TO 720
c        WRITE(63,1062) NN, (NEGRDG(NN,N1),N1=1,6), NMAX(NN)
c        WRITE(62,1063) NN, (NEGRD(NN,MM),MM=1,20)
 720  CONTINUE
 1062 FORMAT(1x,I5,6(1x,I4), 6x,I4)
 1063 FORMAT(1x,I5,20I4)
C 
 9999 WRITE(6,*)  'Data length; ',N-1
      WRITE(6,*)  'Missing data number: ', NNN-1
      WRITE(6,*)  'Rejected data number: ',NNOB-1
C
      
      RETURN
      END
C
      SUBROUTINE TOM(DLMD,DPHD,WBD,SBD,APHD,ALMD,KMIN)
C     ******************************************************************
C     *                                                                *
C     *     ROUTINE FOR CALCULATING THE K INDEX OF THE NEAREST         *
C     *       HEIGHT POINT ON THE E GRID TO THAT OF A POINT            *
C     *                GIVEN IN GEODETIC COORDINATES                   *
C     *  LAT POS NORTH    LON POS WEST                                 *
C     ******************************************************************
      INCLUDE "parmeta"
      INTEGER IM
      REAL DLMD,DPHD,WBD,SBD
                             PARAMETER
     1 (IM1=IM-1,NINC=2*IM-1,DTR=1.745329E-2,NMAX=4)
       DIMENSION AH(NMAX),KH(NMAX)
       WB=WBD*DTR
       SB=SBD*DTR
       DLM=DLMD*DTR
       DPH=DPHD*DTR
C----------------------------------------------------------------------
      AIM1=REAL(IM1)
      RDLM=1./DLM
      RDPH=1./DPH
C
C***  CONVERT LOCATION FROM GEODETIC TO TRANSFORMED LAT/LON
C
CJLG      CALL TLLRAD(APHD,ALMD,TLM,TPH)
      CALL TLLRAD(TLM0D,TPH0D,APHD,ALMD,TLM,TPH)
C
C***  X1,Y1 ARE THE COORDINATES WITH RESPECT TO THE X,Y TRANSFORMED
C***  SYSTEM
C
      X1=(TLM-WB)*RDLM
      Y1=(TPH-SB)*RDPH
C
C***  X2,Y2 ARE THE COORDINATES WITH RESPECT TO THE X',Y' TRANSFORMED
C***  SYSTEM (TRANSLATED BY AIM1 ALONG Y')
C
      X2=.5*( X1+Y1)
      Y2=.5*(-X1+Y1)+AIM1
C
C***  I2,J2 ARE THE COORDINATES OF THE SOUTHERN HEIGHT POINT IN
C***  A CLUSTER OF FOUR
C
      I2=INT(X2)
      J2=INT(Y2)
C
      P=X2-I2
      Q=Y2-J2
      PQ=P*Q
C
      JR=J2-IM1
      I3=I2-JR
      J3=I2+JR
C***
C***  INDEX K CORRESPONDING TO POINT (I2,J2) (SOUTHERNMOST OF THE
C***  FOUR SURROUNDING HEIGHT POINTS)
C***
      K=J3*IM-J3/2+(I3+2)/2
C***
C***  FIND WHICH OF THE FOUR H POINTS IS NEAREST
C***
      AH(1)=P*Q
      AH(2)=(1.-P)*Q
      AH(3)=P*(1.-Q)
      AH(4)=(1.-P)*(1.-Q)
      KH(1)=K
      KH(2)=K+IM
      KH(3)=K+IM1
      KH(4)=K+NINC
C
      KMIN=KH(1)
      NM=1
  100 NL=NMAX1-NM
      DO 200 NN=1,NL
      NX=NM+NN
      IF(AH(NM).GT.AH(NX))THEN
        KMIN=KH(NX)
        NM=NX
        IF(NX.LT.NMAX1)GO TO 100
      ENDIF
  200 CONTINUE
C----------------------------------------------------------------------
      RETURN
      END
C
C----------------------------------------------------------------------
      SUBROUTINE TLLRAD(TLM0D,TPH0D,APHD,ALMD,TLM,TPH)
C     ******************************************************************
C     *                                                                *
C     *  ROUTINE FOR CONVERTING FROM GEODETIC COORDINATES (RADIANS)    *
C     *            TO TRANSFORMED COORDINATES (RADIANS)                *
C     *                                                                *
C     ******************************************************************
                           P A R A M E T E R
C
     1 (DTR=1.745329E-2)
       TLM0=TLM0D*DTR
       TPH0=TPH0D*DTR
C---------------------------------------------------------------------
      STPH0=SIN(TPH0)
      CTPH0=COS(TPH0)
      ALM=ALMD*DTR
      APH=APHD*DTR
      RELM=ALM-TLM0
      SRLM=SIN(RELM)
      CRLM=COS(RELM)
C
      SPH=SIN(APH)
      CPH=COS(APH)
      CC=CPH*CRLM
C
      ANUM=CPH*SRLM
      DENOM=CTPH0*CC+STPH0*SPH
C----------------------------------------------------------------------
      TLM=-ATAN2(ANUM,DENOM)
      TPH=ASIN(CTPH0*SPH-STPH0*CC)
C----------------------------------------------------------------------
      RETURN
      END

