      PROGRAM DRAWVEGETA
C---------------------------------------------------------------
#include "all.inc"
C---------------------------------------------------------------
      PARAMETER (IMB=-WBD/DLMD+1+0.5,JMB=2*(-SBD/DPHD)+1+0.5
     &          ,IM=2*IMB-1,JM=JMB/2+1,IMJM=IMB*JMB-JMB/2)
C
      PARAMETER(NCOLS=14)
      PARAMETER (LRWK=50000,LIWK=50000,LAMA=6000000)

      DIMENSION IBATS(IMJM),ISIB(IMJM),ISIB2D(IM,JM)
      dimension isibij(imb,jmb)

      DIMENSION LHGT(IM,JM)
      DIMENSION ALMD(IM,JM),APHD(IM,JM)
      DIMENSION TALMD(IM,JM),TAPHD(IM,JM)
      DIMENSION XCRA(5),YCRA(5)
      DIMENSION TXCRA(5),TYCRA(5)
C
      DIMENSION XBN(5),YBN(5)
      DIMENSION XBNM(5),YBNM(5)
C
      CHARACTER*3 CHAR_IT
      DIMENSION LND(NCOLS)
      CHARACTER*2 LLB(NCOLS)
C
      DIMENSION IAMA(LAMA), RWRK(LRWK),IWRK(LIWK)
      DIMENSION IASF(120)
      DATA IASF / 120*1 /
C
C Define boundary coordinates in metacode
C
      DATA XBNM/0,32767,32767,0,0/,YBNM/0,0,32767,32767,0/
C
      OPEN(UNIT=23,FILE='../tmp/vegeta.dat',STATUS='UNKNOWN'
     1,    FORM='FORMATTED')
      DO K=1,IMJM
       KK=K
cioannis      READ (23,*)KDUM,IBATS(K),ISIB(K)
      READ (23,*)KK,ISIB(K) ! 30sec
c      if(isib(k).eq.13) print*,k,isib(k)  ! ioannis
      END DO
      CLOSE (23)

      CALL REOR(ISIB,ISIB2D,IM,JM,IMJM) ! transform (imjm) into (imt,jmt)

c ************ transform (imjm) into (im,jm) ********
c      call iconh12(ISIB,ISIBij)
c      print*,'********* ',imb,jmb
c      do j=1,jmb
c         do i=1,imb
c            if (isibij(i,j).eq.13) print*,i,j,isibij(i,j)
c         enddo
c      enddo
c ****************************************************
C
C-----------------------------------------------------------------
C------ G R A P H I C S ------------------------------------------
C-----------------------------------------------------------------

      CALL OPNGKS
      CALL DFCLRS_SSIB
      CALL GSFAIS (1)
      CALL ARINAM (IAMA,LAMA)
C
C Conversion boundary coordinates in fractional coordinates
C
      DO I=1,5
         XBN(I)=CMFX(XBNM(I))
         YBN(I)=CMFY(YBNM(I))
      END DO
C
C Fill whole plot area in white
C
      CALL GSFACI (1)
      CALL GFA (5,XBN,YBN)
C
      CALL MAPPOS (0.1,0.9,0.1,0.9)
      CALL MAPROJ ('CE',TPH0D,TLM0D,0.)
      CALL RTLLNEW (WBD,SBD,TLM0D,TPH0D,RLONMN,RLATMN)
      CALL RTLLNEW (-WBD,-SBD,TLM0D,TPH0D,RLONMX,RLATMX)
      CALL MAPSET ('CO',RLATMN,RLONMN,RLATMX,RLONMX)
      CALL MAPINT

      CALL GSPLCI (0)
      CALL GSTXCI (0)
      CALL SFSETI ('TYPE OF FILL',0)
      CALL GSPLCI(1)
C
C Draw label bar
C
      DO I=1,NCOLS
         LND(I)=I+2
         WRITE(LLB(I),'(I2)')I
      END DO
C
      call LBSETI('CLB',10)
      call LBSETI('CBL',10)
      call LBSETR('WLB',2.0)
      call LBSETR('WBL',2.0)

      CALL GETSET (VL,VR,VB,VT,dum,dum,dum,dum,idum)
      CALL LBLBAR
     &   (1,VR+0.02,1.0,VB,VT,NCOLS,.3333,1.,LND,1,LLB,NCOLS,1)
C
      CALL GSPLCI(0)
      DO I=2,IM-2
      DO J=2,JM-2
         IF(J.EQ.JM.AND.(MOD(I,2).EQ.0)) go to 300

         ALMD(I,J)=WBD+(I-1)*DLMD
         IF(MOD(I,2).EQ.0) THEN
            APHD(I,J)=SBD+(J-1)*2*DPHD + DPHD
         ELSE
            APHD(I,J)=SBD+(J-1)*2*DPHD
         END IF

         XCRA(1)= ALMD(I,J)-DLMD
         XCRA(2)= ALMD(I,J)
         XCRA(3)= ALMD(I,J)+DLMD
         XCRA(4)= ALMD(I,J)
         XCRA(5)= ALMD(I,J)-DLMD
 
         YCRA(1)= APHD(I,J)
         YCRA(2)= APHD(I,J)-DPHD
         YCRA(3)= APHD(I,J)
         YCRA(4)= APHD(I,J)+DPHD
         YCRA(5)= APHD(I,J)

        CALL RTLLNEW (ALMD(I,J),APHD(I,J),TLM0D,TPH0D,
     &                TALMD(I,J),TAPHD(I,J))

        DO N=1,5
         CALL RTLLNEW (XCRA(N),YCRA(N),TLM0D,TPH0D,TXCRA(N),TYCRA(N))
        CALL MAPTRA(TYCRA(N),TXCRA(N),TXCRA(N),TYCRA(N))
        END DO

        CALL MAPTRA(TAPHD(I,J),TALMD(I,J),X,Y)
        IF (X.NE.1.E12) THEN
           CALL POINTS (X, Y, 1, -1, 0)
           CALL GSFACI (ISIB2D(I,J)+2)
           CALL GFA (5,TXCRA,TYCRA)
        END IF

  300 END DO
      END DO

      CALL GSPLCI(0)
        CALL MPLNAM ('Earth..1',4,IAMA)
        CALL MAPSTI ('LA',0)
        CALL MAPSTI ('MV',1)
        CALL GSPLCI (1)
        CALL MAPLBL
        CALL MPLNDR ('Earth..1',4)
        CALL MAPGRD
      CALL FRAME
      CALL CLSGKS

      STOP
      END
C-----------------------------------------------------------------------
      SUBROUTINE TLL (ALMD,APHD,TLMD,TPHD,TPH0D,TLM0D)
C-----------------------------------------------------------------------
      DTR=.01745329
C
      TPH0=TPH0D*DTR
      CTPH0=COS(TPH0)
      STPH0=SIN(TPH0)
C
      RELM=(ALMD-TLM0D)*DTR
      SRLM=SIN(RELM)
      CRLM=COS(RELM)
      APH=APHD*DTR
      SPH=SIN(APH)
      CPH=COS(APH)
      CC=CPH*CRLM
      ANUM=CPH*SRLM
      DENOM=CTPH0*CC+STPH0*SPH
C
      TLMD=ATAN2(ANUM,DENOM)/DTR
      TPHD=ASIN(CTPH0*SPH-STPH0*CC)/DTR
C
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE RTLL(TLMD,TPHD,ALMD,APHD,TPH0D,TLM0D)
C-----------------------------------------------------------------------
      DTR=.01745329
C
      TPH0=TPH0D*DTR
      CTPH0=COS(TPH0)
      STPH0=SIN(TPH0)
C
      TLMR=TLMD*DTR
      TPHR=TPHD*DTR
      STPH=SIN(TPHR)
      CTPH=COS(TPHR)
      CTLM=COS(TLMR)
      STLM=SIN(TLMR)
C
      APH=ASIN(STPH0*CTPH*CTLM+CTPH0*STPH)
      CPH=COS(APH)
C
      ALMD=TLM0D+ASIN(STLM*CTPH/CPH)/DTR
      APHD=APH/DTR
C
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE DFCLRS_BATS
      PARAMETER (NCLRS=23)
      DIMENSION RGBV(3,NCLRS)
      DATA RGBV /
     &     1.00 , 1.00 , 1.00 ,
     1     0.00 , 0.00 , 0.00 ,

     1       .00 , 1.00 ,  .00 , !  green                       35
     2       .13 ,  .55 ,  .13 , !  ForestGreen                 41
     3       .82 ,  .41 ,  .12 , !  chocolate                   60
     4       .96 ,  .64 ,  .38 , !  SandyBrown                  58
     5       .80 ,  .52 ,  .25 , !  peru                        55
     6       .55 ,  .27 ,  .07 , !  SaddleBrown                 53
     7       .00 ,  .39 ,  .00 , !  DarkGreen                   27
     8      1.00 , 1.00 ,  .00 , !  yellow                      48
     9      1.00 ,  .00 , 1.00 , !  magenta                     81
     1       .60 ,  .80 ,  .20 , !  YellowGreen                 40
     1       .93 ,  .87 ,  .51 , !  LightGoldenrod              50
     2       .88 ,  .88 ,  .88 , !  GRAY_3                      97
     3       .58 ,  .44 ,  .86 , !  MediumPurple                87
     4       .00 ,  .75 , 1.00 , !  DeepSkyBlue                 14
     5       .00 ,  .00 , 1.00 , !  Blue                         1
     6       .74 ,  .72 ,  .42 , !  DarkKhaki                   43
     7       .93 ,  .91 ,  .67 , !  PaleGoldenrod               45
     8       .65 ,  .16 ,  .16 , !  brown                       62
     9       .91 ,  .59 ,  .48 , !  DarkSalmon                  63
     1       .00 , 1.00 , 1.00 , !  cyan                        22
     1       .00 ,  .00 ,  .00 / !  black                        22
C
      DO 101 I=1,NCLRS
         CALL GSCR (1,I,RGBV(1,I),RGBV(2,I),RGBV(3,I))
  101 CONTINUE
C
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE DFCLRS_SSIB
      PARAMETER (NCLRS=17)
      DIMENSION RGBV(3,NCLRS)
      DATA RGBV /
     &     1.00 , 1.00 , 1.00 ,
     1     0.00 , 0.00 , 0.00 ,

     1       .55 ,  .27 ,  .07 , !  SaddleBrown                 53
     2       .80 ,  .52 ,  .25 , !  peru                        55
     3       .65 ,  .16 ,  .16 , !  brown                       62
     4       .82 ,  .41 ,  .12 , !  chocolate                   60
     5       .96 ,  .64 ,  .38 , !  SandyBrown                  58
     6       .74 ,  .72 ,  .42 , !  DarkKhaki                   43
     7       .00 ,  .00 ,  .00 , !  black                       22
     8       .00 ,  .00 ,  .00 , !  black                       22
     9       .13 ,  .55 ,  .13 , !  ForestGreen                 41
     1       .00 ,  .39 ,  .00 , !  DarkGreen                   27
     1      1.00 , 1.00 ,  .00 , !  yellow                      48
     2       .60 ,  .80 ,  .20 , !  YellowGreen                 40
     3       .88 ,  .88 ,  .88 , !  GRAY_3                      97
     4       .00 ,  .00 , 1.00 , !  Blue                         1
     1       .00 ,  .00 ,  .00 / !  black                        22
C
      DO 101 I=1,NCLRS
         CALL GSCR (1,I,RGBV(1,I),RGBV(2,I),RGBV(3,I))
  101 CONTINUE
C
      RETURN
      END


C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
                     SUBROUTINE REOR(A1D,A2D,IMT,JMT,IMJM)
C***
C***  TRANSPOSE THE ONE-DIMENSIONAL LAYER DATA TO THE OLD
C***  TWO-DIMENSIONAL CONFIGURATION SINCE THE INTERPOLATION
C***  ROUTINES ARE DESIGNED FOR THAT.
C***
C-----------------------------------------------------------------------
                        D I M E N S I O N
     1 A1D(IMJM),A2D(IMT,JMT)
C-----------------------------------------------------------------------
      N=0
      IBEG=1
      DO 100 J=1,JMT
      DO 50 I=IBEG,IMT,2
      N=N+1
      A2D(I,J)=A1D(N)
   50 CONTINUE
      IF(J.EQ.JMT)GO TO 100
      I2=3-IBEG
      DO 75 I=I2,IMT,2
      N=N+1
      A2D(I,J)=A1D(N)
   75 CONTINUE
  100 CONTINUE
      RETURN
      END

      SUBROUTINE RTLLNEW (TLMD,TPHD,TLM0D,TPH0D,ALMD,APHD)
!-----------------------------------------------------------------------
      parameter (dtr=.01745329)
!-----------------------------------------------------------------------
      TPH0=TPH0D*DTR
!
      CTPH0=COS(TPH0)
      STPH0=SIN(TPH0)
!
      TLMR=TLMD*DTR
      TPHR=TPHD*DTR
!
      STLM=SIN(TLMR)
      CTLM=COS(TLMR)
      STPH=SIN(TPHR)
      CTPH=COS(TPHR)
!
      SPH=CTPH0*STPH+STPH0*CTPH*CTLM
      APH=ASIN(SPH)
      APHD=APH/DTR
      ANUM=CTPH*STLM
      DENOM=(CTLM*CTPH-STPH0*SPH)/CTPH0
      RELM=ATAN2(ANUM,DENOM)
      ALMD=RELM/DTR+TLM0D
!
      IF(ALMD.GT. 180.)    ALMD=ALMD-360.
      IF(ALMD.LT.-180.)    ALMD=ALMD+360.
!
      RETURN
      END 
                         SUBROUTINE ICONH12(IH1,IH2)
C     ******************************************************************
C     *                                                                *
C     *  ROUTINE FOR REORDERING THIBUE-HEIGHT-POINT-1-DIMENSIONAL      *
C     *  MATRICES FOR 2-DIMENSIONAL INDEXING                           *
C     *                                                                *
C     ******************************************************************
#include "all.inc"
                        P A R A M E T E R
     &(IM=-WBD/DLMD+1.5,JM=-2*SBD/DPHD+1.5,IMJM=IM*JM-JM/2)
C-----------------------------------------------------------------------
                             D I M E N S I O N
     & IH2    (IM,JM),IH1    (IMJM)
C-----------------------------------------------------------------------
      print*,'program in ICONH'
      K=0
              DO 100 J=1,JM
          DO 100 I=1,IM-1+MOD(J,2)
      K=K+1
  100 IH2(I,J)=IH1(K)
              DO 101 J=2,JM-1,2
  101 IH2(IM,J)=IH2(IM-1,J)
C-----------------------------------------------------------------------
                             RETURN
                             END
