      PROGRAM VEGETANEW
C------------------------------------------------------------
C
C     CODE TO PRODUCE ETA VEGETATION;
C     PROGRAM WRITTEN USING THE STRUCTURE OF mountnew.F
C     WRITTEN BY F. MESINGER, W. COLLINS AND Z. JANJIC.
C
C************************************************************
C
C     PROGRAM HISTORY LOG:
C     Jul 1998 A. PAPADOPOULOS - ADJUST THE STRUCTURE OF vegeta.F
C                                IN ORDER TO USE THE 30'' DATA SET
C************************************************************

#include "mountnew.inc"
C---------------------
      INCLUDE "vegparm.inc"
C---------------------
      PARAMETER(NVEG=20,DLD=1./120.
     &, IROW=(VEGEBD-VEGWBD)/DLD+1.5 - 1
     &, JCOL=(VEGNBD-VEGSBD)/DLD+1.5 - 1)
      PARAMETER(IVEGSEA=14, PERSEA=0.75)
      DIMENSION NDB(IMJM,NVEG), NVP(IMJM), PER(NVEG)
      DIMENSION IVT(IROW)

      DIMENSION IBATS(IMJM), ISIB(IMJM)
      LOGICAL LAST, INSIDE
      DATA DTR/.01745329/
C
      INCLUDE "../../mount/src/mount_nl.nml"
      INCLUDE "soil_nl.nml"
C
C------------------------------------------------------------
      OPEN (UNIT=11,FILE='../../namelists/name.list',FORM='FORMATTED')
      READ(11,SOIL_NL)
      CLOSE(11)
C------------------------------------------------------------
      WRITE(*,*)'JCOL=',JCOL,'IROW=',IROW

      DO K=1,IMJM
      IBATS(K)=0
      NVP(K)=0
      DO I=1,NVEG  
      NDB(K,I)=0
      ENDDO
      ENDDO
C
C     HERE STARTS PROCESSING ROW BY ROW
C
      LAST=.FALSE.
      NRRD=0
C
      OPEN (UNIT=99,FILE=FN_VEG4ETA,FORM='formatted')
  102 CONTINUE
C
       DO I=1,IROW
       READ(99,'(i2)')IVT(I)
       END DO
C
      ALT0D=VEGNBD-(NRRD-1)*DLD   - DLD/2  
      NRRD=NRRD+1
      ILAST=0
C
      IF (NRRD.EQ.JCOL) LAST=.TRUE.
C
      DO 106 ID=1,IROW
C
C       Vegetation type at this point is IVT(ID)
C
      KVEG=IVT(ID)
C**********************************************************
C ASSUME INLAND WATER TO OCEAN WATER
C
      IF (KVEG .eq. 15) KVEG = 14
C**********************************************************
C
      ALM=(VEGWBD+(ID-1)*DLD)*DTR
      APH=(VEGNBD-(NRRD-1)*DLD)*DTR
      CALL PQK(ALM,APH,INSIDE,P,Q,K)
      IF (.NOT.INSIDE) GO TO 106
      ILAST=ILAST+1
C
C       THE LNDCVR (KVEG) DATA POINT IS WITHIN THE BOX 
C       FORMED BY FOUR NEIGHBORING H POINTS WITH H(K) POINT 
C       AS SOUTH POINT 
C
C                          H(K+NINC)
C
C                    H(K+KNW)      H(K+KNE)
C
C                            H(K)
C
C       NOW CALCULATE THE NUMBER OF POINTS FOR EACH OF 20 
C       VEGETATION TYPES THAT ARE INSIDE THE H GRID-BOX:
C
C
           KK=K
           IF (P.GT.0.5 .AND. Q.LE.0.5) KK=K+KNE
           IF (P.GT.0.5 .AND. Q.GT.0.5) KK=K+NINC
           IF (P.LE.0.5 .AND. Q.GT.0.5) KK=K+KNW

           NDB(KK,KVEG)=NDB(KK,KVEG)+1
           NVP(KK) = NVP(KK) + 1
C
  106 CONTINUE
      LAST=LAST.OR.(ILAST.EQ.0.AND.ALT0D.LT.TPH0D)
      IF (.NOT.LAST) GO TO 102
C---------------------------------------------------------
C     ALL ROWS OF LNDCVR DATA HAVE BEEN PROCESSED.
C     NOW FIND LNDCVR WHICH HAS MAX CONTRIBUTION IN A 
C     H-GRID BOX; MULTIPLY IBATS WITH SEA-MASK FACTOR TO BE 
C     SURE THAT SEA POINTS ARE SEEN AS SEA POINTS IN 
C     THE VEGETATION TYPE DATA
C---------------------------------------------------------
      DO K=1,IMJM
         INIMAX=0
         IINDEX=0
         INIMAX2=0
         IINDEX2=0
         DO I=1,NVEG
            IF (NDB(K,I).GE.INIMAX) THEN
              INIMAX=NDB(K,I)
              IINDEX=I
            ENDIF
         ENDDO
         IBATS(K)=IINDEX
         DO I=1,NVEG
            PER(I) = FLOAT(NDB(K,I))/FLOAT(NVP(K))
         ENDDO
         IF(PER(IVEGSEA).LT.PERSEA .AND. IINDEX.EQ.IVEGSEA) THEN
           DO I=1,NVEG
             IF (NDB(K,I) .NE. NDB(K,IVEGSEA)) THEN
                IF (NDB(K,I).GE.INIMAX2) THEN
                   INIMAX2=NDB(K,I)
                   IINDEX2=I
                ENDIF
             ENDIF
           ENDDO
           IBATS(K)=IINDEX2
cpk           print*,'IINDEX, IINDEX2 :',IINDEX, IINDEX2,NVP(K)
         ENDIF
C
C --- Routine to convert BATS classes to SiB classification
C
         CALL CONBATS2SIB(IBATS(K),ISIB(K))
      ENDDO 
C
      OPEN(UNIT=23,FILE='../tmp/vegeta.dat',STATUS='UNKNOWN')
      DO K=1,IMJM
      WRITE(23,*)K,IBATS(K),ISIB(K)
      END DO
      CLOSE (23)
C
      OPEN(UNIT=23,FILE=FN_VEG
     &        ,STATUS='UNKNOWN',FORM='UNFORMATTED')
      WRITE(23)ISIB
      CLOSE (23)
C
      STOP
      END
C***********************************************************************
      SUBROUTINE CONBATS2SIB(IBATS,ISIB)
      DIMENSION IISIB(21)

      IISIB(1)=  12
      IISIB(2)=  9
      IISIB(3)=  4
      IISIB(4)=  5
      IISIB(5)=  2
      IISIB(6)=  1
      IISIB(7)=  10
      IISIB(8)=  11
      IISIB(9)=  10
      IISIB(10)= 12
      IISIB(11)= 11
      IISIB(12)= 13
      IISIB(13)= 10
      IISIB(14)= 14
      IISIB(15)= 14
      IISIB(16)= 6
      IISIB(17)= 6
      IISIB(18)= 3
      IISIB(19)= 3
      IISIB(20)= 10
      IISIB(21)= 15
C
      ISIB=IISIB(IBATS)
C
      RETURN
      END
C***********************************************************************
                             SUBROUTINE PQK
     &(                    ALM,APH,INSIDE,P,Q,K      )
C     ******************************************************************
C     *                                                                *
C     *  ROUTINE TO TRANSFORM ALAMBDA,APHI (LONGITUDE,LATITUDE) INTO   *
C     *  TRANSFORMED (ROTATED) LONGITUDE,LATITUDE, CHECK WHETHER THE   *
C     *  VALUES OBTAINED ARE WITHIN THE MODEL REGION, AND, IF THEY     *
C     *  ARE, CALCULATE VALUES OF P,Q (COORDINATES WITHIN A SQUARE     *
C     *  FORMED BY CONNECTING FOUR NEIGHBORING HEIGHT POINTS) AND K.   *
C     *                                                                *
C     ******************************************************************
C     *                                                                *
C     *  GRID CONSTANTS:                                               *
C     *  WBD,SBD - WESTERN AND SOUTHERN BOUNDARIES ON LL OR TLL GRID   *
C     *            IN DEGREES                                          *
C     *  TLM0D,TPH0D - ANGLES OF ROTATION OF THE LL COORDINATE SYSTEM  *
C     *            IN THE DIRECTION OF LAMBDA AND PHI RESPECTIVELY     *
C     *            IN ORDER TO OBTAIN THE TLL COORDINATE SYSTEM        *
C     *  DLMD,DPHD - MESH SIDES IN DEGREES                             *
C     *                                                                *
C     ******************************************************************
C
#include "mountnew.inc"
                          L O G I C A L
     &  INSIDE
C
                             D A T A
     &  DTR/.01745329/, D50/0.5/, H1/1./
C
C---------------------------CONSTANTS-----------------------------------
      RIM1=IM1
      TLM0=TLM0D*DTR
      TPH0=TPH0D*DTR
      DLM=DLMD*DTR
      DPH=DPHD*DTR
      RDLM=H1/DLM
      RDPH=H1/DPH
      ALMWB=WBD*DTR !+DLM
      APHSB=SBD*DTR !+DPH
      ALMEB=ALMWB   +2*(IM-1)*DLM
      APHNB=APHSB   +  (JM-1)*DPH
      WB=ALMWB !-DLM
      SB=APHSB !-DPH
          STPH0= SIN(TPH0)
          CTPH0= COS(TPH0)
C
C--------------ALM,APH IS LON,LAT IN RADIAN MEASURE;--------------------
C------------TRANSFORM INTO ROTATED LONGITUDE,LATITUDE------------------
C
          RLM=ALM-TLM0
          SRLM= SIN(RLM)
          CRLM= COS(RLM)
          SAPH= SIN(APH)
          CAPH= COS(APH)
          CC=CRLM*CAPH
C
          ANUM=SRLM*CAPH
          DENOM=CTPH0*CC+STPH0*SAPH
C
          ALM= ATAN2(ANUM,DENOM)
          APH= ASIN(CTPH0*SAPH-STPH0*CC)
C
C-----IS ALM,APH INSIDE THE (ONE LINE REDUCED) MODEL DOMAIN?------------
C
      IF (ALM.LE.ALMWB .OR. ALM.GE.ALMEB .OR.
     &    APH.LE.APHSB .OR. APH.GE.APHNB)  THEN
             INSIDE=.FALSE.
             RETURN
      END IF
             INSIDE=.TRUE.
C
C
C---------X1,Y1 IS A COORDINATE SYSTEM WITH DLM,DPH AS LENGTH UNITS-----
      X1=(ALM-WB)*RDLM
      Y1=(APH-SB)*RDPH
C---------X2,Y2 ROTATED FOR +45 DEG. & TRANSLATED FOR IM1---------------
      X2=D50*( X1+Y1)
      Y2=D50*(-X1+Y1)+RIM1
C---------I2,J2 ARE COORDINATES OF SOUTHERN POINTS OF GRID BOXES--------
      I2=  INT(X2)
      J2=  INT(Y2)
C---------REMAINING PARAMETERS NEEDED TO KNOW THE POSITION--------------
C         OF THE TRANSFORMED POINT WITHIN THE MODEL REGION
      P=X2-I2
      Q=Y2-J2
C-----------------INDEX K CORRESPONDS TO I2,J2--------------------------
      JR=J2-IM1
      I3=I2-JR
      J3=I2+JR
      K=J3*IM-J3/2+(I3+2)/2
C-----------------------------------------------------------------------
                             RETURN
                             END
C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
