C------------------------------------------------------------
C
C     CODE TO PRODUCE ETA SOIL;                
C     PROGRAM WRITTEN USING THE STRUCTURE OF mountnew.f
C     WRITTEN BY F. MESINGER, W. COLLINS AND Z. JANJIC.
C     (MODIFIED BY S. NICKOVIC FOR SURFACE TYPES)
C 
C************************************************************
C
C     PROGRAM HISTORY LOG:
C      Jan 1996 S. NICKOVIC     - ORIGINATOR
C      Jul 1996 A. PAPADOPOULOS - CONVERTION TO DIFFERENT
C                                 SURFACE SCHEMES
C************************************************************
C                                                           *
C     DATA SOURCE: OLSON WORLD ECOSYSTEM                    *
C     WITH 74 CATHEGORIES; RESOLUTION: 10'X10'DEG           *
C     DATA ARRAY STARTS FROM NBD=-90, EBD=-180              *
C                                                           *
C************************************************************
C                                                           *
C     METHOD: THE LARGEST PERCENT OF A VEG  TYPE  THAT      *
C     ENTERS IN THE H GRID-BOX IS USED AS                   * 
C     THE REPRESENTATIVE ONE                                *    
C                                                           *
C************************************************************
C                                                           *
#include "mountnew.inc"
C---------------------
      PARAMETER(NVEG=74,DLD=1./6.
     &, IROW=360./DLD, JCOL=180./DLD )

      BYTE IV(IROW,JCOL)

      DIMENSION NDB(IMJM,NVEG)
      DIMENSION IVEGT(IROW)

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

      DO K=1,IMJM
      RVEG(K)=0.
      DO I=1,NVEG  
      NDB(K,I)=0
      ENDDO
      ENDDO
C
C------------------------------------------------------------
C     NOTE: DATA IS RECORDED AS INTEGER*1
C
      OPEN(10,FILE=FN_OWE,FORM='unformatted'
     &,            ACCESS='direct',RECL=IROW*JCOL)
      READ(10,REC=1)((IV(I,J),I=1,IROW),J=JCOL,1,-1)
      CLOSE(10)
C
C     HERE STARTS PROCESSING ROW BY ROW
C
      LAST=.FALSE.
      NRRD=0
  102 CONTINUE
C
       DO I=1,IROW
          IVEGT(I)=IV(I,NRRD+1) 
       END DO
C
 900   FORMAT(135I3)
C
      ALT0D= -90+(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 IVEGT(ID)
C
      KVEG=IVEGT(ID)
C
C IF olson TYPE IS 0 REDEFINE IT AS 74
C
      if (KVEG .eq. 0) KVEG = 74
C
      ALM=(-180.+ DLD/2. +(ID-1)*DLD)*DTR
      APH=(- 90.+ DLD/2. +(NRRD-1)*DLD)*DTR
      CALL PQK(ALM,APH,INSIDE,P,Q,K)
      IF (.NOT.INSIDE) GO TO 106
      ILAST=ILAST+1
C
C       THE VEGETATION TYPE (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 73 
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
C
  106 CONTINUE
      LAST=LAST.OR.(ILAST.EQ.0.AND.ALT0D.GT.TPH0D)
      IF (.NOT.LAST) GO TO 102
C---------------------------------------------------------
C     ALL ROWS OF VEGETATION DATA HAVE BEEN PROCESSED.
C     NOW FIND VEGETATION TYPE WHICH HAS MAX CONTRIBUTION IN A 
C     H-GRID BOX; MULTIPLY RVEG 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
         DO I=1,NVEG
            IF (NDB(K,I).GE.INIMAX) THEN
              INIMAX=NDB(K,I)
              IINDEX=I
            ENDIF
         ENDDO
         RVEG(K)=IINDEX
      ENDDO 
C-----------------------------------------------------------------------
      OPEN(UNIT=23,FILE='../tmp/vegeta.dat',STATUS='UNKNOWN')
      DO K=1,IMJM
         II=INT(RVEG(K))
         if (II.le.0) then
            write(*,*)k,(NDB(K,I),I=1,NVEG)
            print*, '!!!!', k,RVEG(K),ii
         end if
C
      CALL CONVEG2BATS(INT(RVEG(K)),IBATS(K))
      CALL CONBATS2SSIB(IBATS(K),ISSIB(K))
C
      WRITE(23,*)K,INT(RVEG(K)),ISSIB(K)
!      print*,K,RVEG(K),ISSIB(K)
      END DO
      CLOSE (23)

      OPEN(UNIT=23,FILE=FN_VEG
     &    ,FORM='UNFORMATTED',STATUS='UNKNOWN')
      WRITE(23) ISSIB
      CLOSE(23)
C
      STOP
      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?------------
      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---------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
      SUBROUTINE CONVEG2BATS(IRVEG,IBATS)
      DIMENSION IIBATS(74)

      IIBATS(1)=  2
      IIBATS(2)=  1
      IIBATS(3)=  99
      IIBATS(4)=  99
      IIBATS(5)=  99
      IIBATS(6)=  6
      IIBATS(7)=  99
      IIBATS(8)=  8
      IIBATS(9)=  99
      IIBATS(10)= 99
      IIBATS(11)= 99
      IIBATS(12)= 99
      IIBATS(13)= 99
      IIBATS(14)= 99
      IIBATS(15)= 99
      IIBATS(16)= 16
      IIBATS(17)= 12
      IIBATS(18)= 99
      IIBATS(19)= 99
      IIBATS(20)= 3
      IIBATS(21)= 3
      IIBATS(22)= 3
      IIBATS(23)= 4
      IIBATS(24)= 4
      IIBATS(25)= 5
      IIBATS(26)= 5
      IIBATS(27)= 3
      IIBATS(28)= 6
      IIBATS(29)= 5
      IIBATS(30)= 1
      IIBATS(31)= 1
      IIBATS(32)= 5
      IIBATS(33)= 6
      IIBATS(34)= 99
      IIBATS(35)= 99
      IIBATS(36)= 13
      IIBATS(37)= 10
      IIBATS(38)= 10
      IIBATS(39)= 10
      IIBATS(40)= 9
      IIBATS(41)= 2
      IIBATS(42)= 9
      IIBATS(43)= 7
      IIBATS(44)= 13
      IIBATS(45)= 13
      IIBATS(46)= 16
      IIBATS(47)= 17
      IIBATS(48)= 16
      IIBATS(49)= 6
      IIBATS(50)= 8
      IIBATS(51)= 11
      IIBATS(52)= 11
      IIBATS(53)= 9
      IIBATS(54)= 18
      IIBATS(55)= 9
      IIBATS(56)= 17
      IIBATS(57)= 2
      IIBATS(58)= 2
      IIBATS(59)= 17
      IIBATS(60)= 18
      IIBATS(61)= 18
      IIBATS(62)= 18
      IIBATS(63)= 2
      IIBATS(64)= 13
      IIBATS(65)= 2
      IIBATS(66)= 2
      IIBATS(67)= 2
      IIBATS(68)= 2
      IIBATS(69)= 9
      IIBATS(70)= 12
      IIBATS(71)= 8
      IIBATS(72)= 6
      IIBATS(73)= 14
      IIBATS(74)= 15

C
      IBATS=IIBATS(IRVEG)
C
      RETURN
      END

C


      SUBROUTINE CONBATS2SSIB(IBATS,ISSIB)
      DIMENSION IISSIB(18)

      IISSIB(1)=  12
      IISSIB(2)=  9
      IISSIB(3)=  4
      IISSIB(4)=  5
      IISSIB(5)=  2
      IISSIB(6)=  1
      IISSIB(7)=  10
      IISSIB(8)=  11
      IISSIB(9)=  10
      IISSIB(10)= 12
      IISSIB(11)= 11
      IISSIB(12)= 13
      IISSIB(13)= 10
      IISSIB(14)= 14
      IISSIB(15)= 14
      IISSIB(16)= 6
      IISSIB(17)= 6
      IISSIB(18)= 3

C
      ISSIB=IISSIB(IBATS)
C
      RETURN
      END
C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
