      PROGRAM DRAWTEXTETA
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=11)
      PARAMETER (LRWK=50000,LIWK=50000,LAMA=6000000)

      DIMENSION SOIL(IMJM),SOIL2D(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
      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/texteta.dat',STATUS='UNKNOWN'
     1,    FORM='FORMATTED')
      DO K=1,IMJM
      READ (23,*)KDUM,SOIL(K)
      END DO
      CLOSE (23)
      CALL REOR(SOIL,SOIL2D,IM,JM,IMJM)
C
C-----------------------------------------------------------------
C------ G R A P H I C S ------------------------------------------
C-----------------------------------------------------------------

      CALL OPNGKS
      CALL DFCLRS
      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)
C
C Draw label bar
C
      DO I=1,NCOLS
         LND(I)=I+1
         WRITE(LLB(I),'(I2)')I
      END DO
C
      CALL SFSETI ('TYPE OF FILL',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,0,LLB,NCOLS,1)
C
      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 (INT(SOIL2D(I,J))+1)
           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
      PARAMETER (NCLRS = 13)
      DIMENSION RGBV(3,NCLRS)
      DATA RGBV /
     1      1.00 , 1.00 , 1.00 , !  White
     2       .00 ,  .39 ,  .00 , !  DarkGreen
     3       .24 ,  .70 ,  .44 , !  MediumSeaGreen
     4       .50 , 1.00 ,  .83 , !  aquamarine
     5       .42 ,  .56 ,  .14 , !  OliveDrab
     6       .55 ,  .27 ,  .07 , !  SaddleBrown
     7       .80 ,  .52 ,  .25 , !  peru
     8       .68 , 1.00 ,  .18 , !  GreenYellow
     9       .93 ,  .91 ,  .67 , !  PaleGoldenrod
     1       .93 ,  .87 ,  .51 , !  LightGoldenrod
     1       .00 ,  .00 , 1.00 , !  Blue
     2       .83 ,  .83 ,  .83 , !  GRAY_4
     3      1.00 , 1.00 ,  .00 / !  yellow

      CALL GSCR (1,0,0.,0.,0.)

      DO I=1,NCLRS
         CALL GSCR (1,I,RGBV(1,I),RGBV(2,I),RGBV(3,I))
      ENDDO

      RETURN
      END

C-----------------------------------------------------------------------

      SUBROUTINE DFCLRS_AUTO(N)
C
      CWCOMP = 1.4
      CWROT = 1.0
      RFACT = 0.5 *CWROT
C
        CALL GSCR (1,0,0.,0.,0.)
        CALL GSCR (1,1,1.,1.,1.)
C
        DO 101 I=2,N
        S=FLOAT(I-1)/FLOAT(N-1)
        T= CWCOMP*(S-RFACT)
         R=128.0 + 127.0 * ATAN( 7.0*T ) / 1.57
         G=128.0 + 127.0 * (2 * EXP(-7*T*T) - 1)
         B=128.0 + 127.0 * ATAN( -7.0*T ) / 1.57
         R=R/255
         G=G/255
         B=B/255
          CALL GSCR (1,I,R,G,B)
  101   CONTINUE
        RETURN
C
      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 
C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
        SUBROUTINE FILL (XWRK,YWRK,NWRK,IAREA,IGRP,NGRPS)
        DIMENSION XWRK(*),YWRK(*),IAREA(*),IGRP(*)

        DO 10, I=1,NGRPS
          IF (IGRP(I).EQ.3) IAREA3=IAREA(I)
   10   CONTINUE

        IF (IAREA3 .GT. 0) THEN
C If the area is defined by 3 or more points, fill it
           CALL GSFACI(IAREA3+1)
           CALL GFA(NWRK,XWRK,YWRK)
        ENDIF

C Otherwise, do nothing

        RETURN
        END
C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
      SUBROUTINE DRAWCL (XCS,YCS,NCS,IAI,IAG,NAI)
C
C This version of DRAWCL draws the polyline defined by the points
C ((XCS(I),YCS(I)),I=1,NCS) if and only if none of the area identifiers
C for the area containing the polyline are negative.  The dash package
C routine CURVED is called to do the drawing.
C
        DIMENSION XCS(*),YCS(*),IAI(*),IAG(*)
C
C Turn on drawing.
C
        IDR=1
C
C If any area identifier is negative, turn off drawing.
C
        DO 101 I=1,NAI
          IF (IAI(I).LT.0) IDR=0
  101   CONTINUE
C
C If drawing is turned on, draw the polyline.
C
        IF (IDR.NE.0) CALL CURVED (XCS,YCS,NCS)
C
C Done.
C
      RETURN
      END
C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

