C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& 
      PROGRAM DRAW
C-----------------------------------------------------------------------
#include "all.inc"

c      parameter(TLM0D=15.0,TPH0D=39.5)
c      parameter(WBD=-24.6,SBD=-11.5,DLMD=0.1,DPHD=0.1,DTB=30.,LM=38)
C-----------------------------------------------------------------------
       PARAMETER(LMP1=LM+1)
       PARAMETER
     &(IM=-1*WBD/DLMD+1.5,JM=-2*SBD/DPHD+1.5
     &,IM1=IM-1)
       PARAMETER
     &( IMM1=IM1,IMJM=IM*JM-JM/2
     &, IMT=2*IM-1,JMT=JM/2+1)
       PARAMETER
     & (KNE=IM,KNW=IM1,KSE=-KNW,KSW=-KNE
     &, NINC=2*IM-1,KHL00=1,KHH00=IM*JM-JM/2,JM2=JM-2
     &, KHL22=2*IM+1,KHH22=IM*(JM-2)-(JM-2)/2-1  )

C-----------------------------------------------------------------------
c      PARAMETER (IWEST=0,IEAST=0,JNORTH=0,JSOUTH=0)
c      parameter(ZWBD=WBD+1.2, ZSBD=SBD+0.6)
      parameter(ZWBD=WBD, ZSBD=SBD)
c      parameter(ZWBD=WBD+0.1, ZSBD=SBD+0.3)
      parameter (ZLM0D = TLM0D, ZPH0D = TPH0D)

      parameter (DLM0D = -TLM0D+ZLM0D, DPH0D = -ZPH0D+TPH0D)
      parameter (IWEST = (-WBD+ZWBD+DLM0D)/(2*DLMD) + 0.5)
      parameter (IEAST = (-WBD+ZWBD-DLM0D)/(2*DPHD) + 0.5)
      parameter (JNORTH = 2*(-SBD+ZSBD+DPH0D)/(2*DPHD)+0.5)
      parameter (JSOUTH = 2*(-SBD+ZSBD-DPH0D)/(2*DPHD)+0.5)
      PARAMETER (IME=IM-IWEST-IEAST)
      PARAMETER (JME=JM-JNORTH-JSOUTH)
      PARAMETER (IMA=IME*2-1,JMA=JME)
C-----------------------------------------------------------------------
      PARAMETER (NCONT=15)  ! 15
      PARAMETER (DTR=0.01745329)
      PARAMETER (rundays=17.0)
C-----------------------------------------------------------------------
      DIMENSION DUM2(IM,JM),DFL(LMP1),VALUE(NCONT+1)
      DIMENSION S(IM,JM,LM),SINT(IM,JM),LMH(IM,JM),TMP(IM,JM)
      DIMENSION DEPD(IM,JM),DEPW(IM,JM),DEPT(IM,JM)
      DIMENSION SFLX(IM,JM),USTAR(IM,JM),SND(IM,JM)
      DIMENSION SZ0(IM,JM),S10(IM,JM),AKHS(IM,JM)
C
      dimension as(im,jm,lm,3),adepd(im,jm,3),adepw(im,jm,3)
      
      dimension asb(im,jm,lm,3),asv(im,jm,lm,3),ass(im,jm,3)

      dimension adepdb(im,jm,2),adepwb(im,jm,2)
      dimension adepw2(im,jm,3),adepwb2(im,jm,2)
      dimension vvv(im,jm,3),gas(im,jm,lm,3)
      dimension seap(im,jm),sean(im,jm)
      dimension pd(im,jm),t(0:im+1,0:jm+1,lm)
      dimension amerc(im,jm,lm,3),bmerc(im,jm,lm,3)
      dimension vmerc(im,jm,lm,3)

      DIMENSION APHD(IM,JM),TAPHD(IM,JM),ALMD(IM,JM),TALMD(IM,JM) 

      DIMENSION IDAT(3),ITERMIN(9)

      DIMENSION VARA(IMA,JMA)
      DIMENSION VARE(IME,JME)

      DIMENSION DUM1D(LM+1),DUM1DD(10)
      DIMENSION DUM2D(IM,JM),SM(IM,JM),FIS(IM,JM),HGT(IM,JM)
      dimension hgtk(imjm),hgtt(imt,jmt)
      dimension smk(imjm),smt(imt,jmt)
      DIMENSION DUM3D(IM,JM,LM),DUM3DD(IM,JM,LM-1)
      DIMENSION U10(IM,JM),V10(IM,JM)

      dimension dum2d1(im,jm),dum2d2(im,jm),dum2d3(im,jm)
      dimension dum3d1(im,jm,LM),dum3d2(im,jm,LM),dum3d3(im,jm,LM)

      dimension t2(im,jm,lm),q2(im,jm,lm)
      dimension u2(im,jm,lm),v2(im,jm,lm)
      dimension acprec2(im,jm),cuprec2(im,jm),acprec(im,jm)
      dimension idat2(3),T2M(im,jm),TPOT2M(im,jm)
      DIMENSION 
     & RES(IM,JM)
     &,ETA(LMP1)
     &,HTM(IM,JM,LM)
     &,VTM(IM,JM,LM)
     &,LMV(IM,JM)
     &,SICE(IM,JM)
     &,SST(IM,JM),IVGTYP(IM,JM),ISLTYP(IM,JM)


      CHARACTER*140 INFN,FILEIN,INFN2

      character*6 vmin,vmax ! ioannis
      CHARACTER*1 CVAR
      CHARACTER*70 TITLE1, TITLE2, TITLE3
      DATA ITERMIN / 01, 06, 12, 18, 24, 30, 36, 42, 48/
C---------------------------------------------------------------------
      dimension month(12)
      data month/31,28,31,30,31,30,31,31,30,31,30,31/
C---------------------------------------------------------------------
      COMMON/MAPCOOR/XSTA,XSTO,YSTA,YSTO,YLBL

      if(MOD(IWEST,2).ne.0.or.MOD(IEAST,2).ne.0) then
         print*,'IWEST - IEAST : ',IWEST,IEAST
         STOP'Change ZWBD'
      endif
      if(MOD(JNORTH,2).ne.0.or.MOD(JSOUTH,2).ne.0) then
         print*,'JSOUTH - JNORTH : ',JSOUTH,JNORTH
         STOP'Change ZSBD'
      endif
      XSTA = 0.05
      XSTO = 0.95
      YSTA = 0.10 ! 0.05
      YSTO = 0.95
      YLBL = 0.15

      CVAR='1'          
C-----------------------------------------------------------------------
        print*,infn
        print*,'IM=',im,' JM=',jm
      INFN=
     $ '../../../../../../output/PERM/etatopo.dat'

        OPEN(UNIT=31,FILE=INFN
     &   ,STATUS='UNKNOWN',ACCESS='SEQUENTIAL',FORM='UNFORMATTED')
c        READ(31) HGTT,SMT
        READ(31) HGT,SM
        CLOSE(31)
C-----------------------------------------------------------------------
C********************* TEST *********************************  
      print*,'iwest=',iwest,' ieast=',ieast
      print*,'jnorth=',jnorth,' jsouth=',jsouth
	write(*,*) im,jm

c      open(53,file='./fis_orog-real_000.dat')
c      do j=jm,1,-1
c      write(53,'(i3,39f5.0))') j,(fis(i,j)/9.8,i=1,im)
c      write(53,'(i3,39f3.0))') j,(sm(i,j),i=1,im)
c      enddo
c      close (53)

c       do i=1,im
c      write(55,1)(acprec2(i,j)*1000., j=1,jm)
c1     format(1x,100F3.0)
c       enddo

c************************************************************
C
C-----------------------------------------------------------------------
      CALL OPNGKS
      CALL GSCLIP (0)
C-----------------------------------------------------------------------
C ********** Calculate the lat-long of each grid point **************
C ** the lat-longs are stored in the arrays TALMD(i,j), TAPHD(i,j) **
C *******************************************************************

      TPH0=TPH0D*DTR
      TLM0=TLM0D*DTR
      CTPH0=COS(TPH0)
      STPH0=SIN(TPH0)

         DO J=1,JM
         DO I=1,IM
            IF((I.EQ.IM).AND.(MOD(J,2).EQ.0)) go to 300
            APHD(I,J)=SBD+(J-1)*DPHD
            IF(MOD(J,2).EQ.0) THEN
               ALMD(I,J)=WBD+(I-1)*2.*DLMD + DLMD
            ELSE
               ALMD(I,J)=WBD+(I-1)*2.*DLMD
            END IF
            CALL RTLLNEW
     & (ALMD(I,J),APHD(I,J),TLM0D,DTR,CTPH0,STPH0,TALMD(I,J),TAPHD(I,J))
  300    ENDDO
         ENDDO
c **********************************************************************
C-----------------------------------------------------------------------
      IF(CVAR.EQ.'1') THEN
c       CALL CONH21(HGTT,HGTK)
c	CALL CONH12(HGTK,HGT)
c	CALL CONH21(SMT,SMK)
c	CALL CONH12(SMK,SM)
        AMX=0.0
        DO J=1,JM
        DO I=1,IM
c         if (((hgt(i,j)*9.81).gt.29000.0).and.
c     $       ((hgt(i,j)*9.81).lt.29500.0)) print*,i,j,hgt(i,j)*9.81
           sint(i,j)=0.0
           sint(i,j)=HGT(i,j)
c           sint(i,j)=SM(i,j)

c          if (sm(i,j).eq.1.0) sint(i,j)=0.0

           IF (SINT(I,J).GT.AMX) AMX=SINT(I,J)
        END DO
        END DO
        
c        ixval=119
c        iyval=228
c        print*,ixval,iyval
c        print*,'Height=',hgt(ixval,iyval)
c        print*,'Height=',hgt(ixval,iyval)*9.81
c        print*,'Mask=',sm(ixval,iyval)
c        print*,TAPHD(ixval,iyval),TALMD(ixval,iyval)
c        sint(ixval,iyval)=5000.0
        is=im
	js=jm
        call ddd(is,js,sint,ppmin,ppmax)
         
        print*,im,jm,ppmin,ppmax
c        print*,'Give MIN and MAX values for the plot'
c        read*, vlmin,vlmmax,app
        vlmin = -1.
        vlmmax = 4500.
        app = 1.
c         write(vmin(1:6),'(f7.4)') ppmin-0.0001 ! ioannis
c	 read(vmin,*) vlmin              ! ioannis
c         write(vmax(1:6),'(f7.4)') ppmax+0.001 ! ioannis
c         read(vmax,*) vlmmax                   ! ioannis
c         app=1.0                              ! ioannis

c	vlmin=0.0 ! ioannis
c	vlmmax=0.0 ! ioannis

	VLINC=(VLMMAX-VLMIN)/NCONT

c        DO N=1,NCONT+1
c	   VALUE(N)=VLMIN+(N-1)*VLINC
c        ENDDO
        VALUE(1) = -1.
        VALUE(2) = 1
        DO N=3,NCONT+1
	   VALUE(N)=VALUE(N-1)+VLINC
        ENDDO
        print*,VALUE

        DO J=1,JME
           DO I=1,IME
             VARE(I,J)=SINT(I+IWEST,J+JSOUTH)*app
           END DO
        END DO
        CALL ETOAH  (VARE,IME,JME,VARA,IMA,JMA)
        CALL DRAWDL (VARA,IMA,JMA,VALUE,NCONT,CVAR,
     &               IWEST,IEAST,JNORTH,JSOUTH)
C-----------------------------------------------------------------------
      ELSE
         print*,'Wrong CVAR'
         stop 30
C-----------------------------------------------------------------------
      ENDIF
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------

      CALL FRAME
      CALL CLSGKS
      STOP
      END
C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
      subroutine ddd(is,js,sss,pmin,pmax)
      dimension sss(is,js)
      qmax=-100000000000.0
      qmin=10000000000000000.0
      imax=0
      jmax=0
      do i=1,is
      do j=1,js
c      print*,i,j,sss(i,j)
      if(sss(i,j).lt.qmin) then
      qmin=sss(i,j)
      imin=i
      jmin=j
      endif
      if(sss(i,j).gt.qmax) then 
      qmax=sss(i,j)
c      imax=i
c      jmax=j
      endif
      enddo
      enddo
      pmin=qmin
      pmax=qmax
      print*,is,js,pmin,pmax
     $     ,imin,jmin
      return
      end
C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
      SUBROUTINE DRAWDL(ZDAT,IMT,JMT,VALUE,NCONT,CVAR,
     &                  IWEST,IEAST,JNORTH,JSOUTH)
      DIMENSION ZDAT(IMT,JMT)
      DIMENSION VALUE(NCONT+1)
      DIMENSION LND(NCONT+1)
      DIMENSION LLB(NCONT+1)
      CHARACTER*8 CLLB(NCONT+1)
      CHARACTER*1 CVAR
C
C Creating LLB i CLLB...
C
      DO I=1,NCONT+1
        LND(I)=I+2
        VALUEI=VALUE(I)
c        WRITE(CLLB(I),'(F6.4)')VALUEI
c        WRITE(CLLB(I),'(F5.2)')VALUEI
        WRITE(CLLB(I),'(F5.0)')VALUEI
c        WRITE(CLLB(I),'(F6.1)')VALUEI
c        WRITE(CLLB(I),'(F7.5)')VALUEI
c        WRITE(CLLB(I),'(F4.0)')VALUEI
        IF(CVAR.EQ.'4')WRITE(CLLB(I),'(F4.1)')VALUEI*1E09
      END DO
      CALL SCALARMONO
     &    (ZDAT,IMT,JMT,NCONT,VALUE,LND,CLLB,
     &     IWEST,IEAST,JNORTH,JSOUTH)
      RETURN
      END
C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
      SUBROUTINE RTLL(TLM,TPH,ALMD,APHD,tph0d,tlm0d)

      DTR=.01745329

      TPH0=TPH0D*DTR
      CTPH0=COS(TPH0)
      STPH0=SIN(TPH0)

      TLMr=TLM*DTR
      TPHr=TPH*DTR
      STPH=SIN(TPHr)
      CTPH=COS(TPHr)
      CTLM=COS(TLMr)
      STLM=SIN(TLMr)

      APH=ASIN(STPH0*CTPH*CTLM+CTPH0*STPH)
      CPH=COS(APH)

      ALMD=TLM0D+ASIN(STLM*CTPH/CPH)/DTR
      APHD=APH/DTR

      RETURN
      END

C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
      SUBROUTINE RLTLW(ALMD,APHD,TPUS,TPVS,TLM0D,TPH0D,PUS,PVS)                                         
      DTR=.01745329

      TPH0=TPH0D*DTR
      CTPH0=COS(TPH0)
      STPH0=SIN(TPH0)

      RELM=(ALMD-TLM0D)*DTR                                             
      SRLM=SIN(RELM)                                                    
      CRLM=COS(RELM)                                                    

      APH=APHD*DTR                                                      
      SPH=SIN(APH)                                                      
      CPH=COS(APH)                                                      

      CC=CPH*CRLM                                                       
      TPH=ASIN(CTPH0*SPH-STPH0*CC)

      RCTPH=1./COS(TPH)                                                 
      CRAY=STPH0*SRLM*RCTPH                                             
      DRAY=(CTPH0*CPH+STPH0*SPH*CRLM)*RCTPH                             
      DC=DRAY*DRAY+CRAY*CRAY                                            
      PUS=(DRAY*TPUS+CRAY*TPVS)/DC                                      
      PVS=(DRAY*TPVS-CRAY*TPUS)/DC                                      

      RETURN                                                            
      END                                                               

C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
      SUBROUTINE RTLLNEW (TLMD,TPHD,TLM0D,DTR,CTPH0,STPH0,ALMD,APHD)
!     ****************************************************************
!     *                                                              *
!     *  programer: z. janjic, shmz, feb. 1981                       *
!     *  ammended:  z. janjic, ncep, jan. 1996                       *
!     *                                                              *
!     *  transformation from rotated lat-lon to lat-lon coordinates  *
!     ****************************************************************
!     ****************************************************************
!     *  tlm   - transformed longitude, rad.                         *
!     *  tph   - transformed latitude, rad.                          *
!     *  tlm0d - the angle of rotation of the transformed lat-lon    *
!     *          system in the longitudinal direction, degs          *
!     *  ctph0 - cos(tph0), tph0 is the angle of rotation of the     *
!     *          transformed lat-lon systemn in the latitudinal      *
!     *          direction, precomputed                              *
!     *  stph0 - sin(tph0), tph0 is the angle of rotation of the     *
!     *          transformed lat-lon systemn in the latitudinal      *
!     *          direction, precomputed                              *
!     *  almd  - geographical longitude, degs, range -180.,180       *
!     *  aphd  - geographical latitude,  degs, range - 90., 90.,     *
!     *          poles are singular                                  *
!     ****************************************************************
!-----------------------------------------------------------------------
!
      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 COLRLN (XCS,YCS,NCS,IAI,IAG,NAI)

      DIMENSION XCS(*),YCS(*),IAI(*),IAG(*)
C
C For each line segment, one gets a set of points (using normalized
C device coordinates), two group identifiers and their associated
C area identifiers.  If both of the area identifiers are zero or
C negative, the segment is not drawn; otherwise, we use MAPACI to
C see if the segment is over water and, if so, we draw the segment.
C If the segment is defined by more than 150 points, we would like to
C know about it.
C
        IF (IAI(1).GE.0.AND.IAI(2).GE.0) THEN
          ITM=MAX0(IAI(1),IAI(2))
          IF (MAPACI(ITM).EQ.1) THEN
            IF (NCS.GT.150) PRINT * , 'COLRLN - NCS TOO BIG - ',NCS
            CALL GPL (NCS,XCS,YCS)
          END IF
        END IF
C
C Done.
C
      RETURN
      END

C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
      SUBROUTINE SCALARMONO(ZDAT,IMT,JMT,NCONT,VALUE,LND,CLLB
     &,IWEST,IEAST,JNORTH,JSOUTH)
#include "all.inc"

c      parameter(TLM0D=15.0,TPH0D=39.5)
c      parameter(WBD=-24.6,SBD=-11.5,DLMD=0.1,DPHD=0.1,DTB=30.,LM=38)
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
      PARAMETER (LRWK=50000,LIWK=50000,LAMA=6000000)
      PARAMETER (NCRA=20000)
C
c      DIMENSION XBN(5),YBN(5)
c      DIMENSION XBNM(5),YBNM(5)
C
      DIMENSION ZDAT(IMT,JMT)
      DIMENSION VALUE(NCONT+1) 
      DIMENSION LND(NCONT+1)
      CHARACTER*(*) CLLB(NCONT+1)
      COMMON/MAPCOOR/XSTA,XSTO,YSTA,YSTO,YLBL
C
C Declare the area map array.
C
      DIMENSION IAMA(LAMA)
C
C Declare the real and integer workspace arrays for CONPACK.
C
      DIMENSION RWRK(LRWK),IWRK(LIWK)
C
C Declare arrays for ARSCAM and MAPGRM to use in calls to COLRAM and
C COLRLL, respectively.  XCRA and YCRA hold X and Y coordinates; IAIA
C and IGIA hold area identifiers and group identifiers.
C
      DIMENSION XCRA(NCRA),YCRA(NCRA),IAIA(240),IGIA(240)
C
C Declare an array in which to retrieve aspect source flags.
C
      DIMENSION IASF(240)
      DATA IASF / 240*1 /
C
C Define boundary coordinates in metacode
C
c      DATA XBNM/0,32767,32767,0,0/,YBNM/0,0,32767,32767,0/
C
      EXTERNAL COLRLN
      EXTERNAL DRAWCL
      EXTERNAL COLRAM
      EXTERNAL COLRCL
      EXTERNAL COLRLL
      EXTERNAL FILL
      DO J = 1, JMT
      DO I = 1, IMT
	 IF (ZDAT(I,J) .LE. VALUE(1)-0.02) ZDAT(I,J)=0.
      ENDDO
      ENDDO
C
C Set all aspect source flags to "individual".
C
      CALL GSCLIP (0)
      CALL GSASF (IASF)
      CALL GSFAIS (1)
C
C Force solid fill.
C
      CALL DFCLRS_INC(NCONT+1)
c     CALL DFCLRS(NCONT+2)
      CALL ARINAM (IAMA,LAMA)
      CALL MAPPOS (XSTA,XSTO,YSTA,YSTO)
      CALL RTLL(WBD+IWEST*2*DLMD ,SBD+JSOUTH*DLMD
     &         ,RLONMN         ,RLATMN         ,TPH0D,TLM0D)
      CALL RTLL(-WBD-IEAST*2*DLMD,-SBD-JNORTH*DLMD
     &         ,RLONMX         ,RLATMX          ,TPH0D,TLM0D)
C
      CALL MAPROJ ('CE',TPH0D,TLM0D,0.)
      CALL MAPSET ('CO',RLATMN,RLONMN,RLATMX,RLONMX)
      CALL MPSETC ('OU','PO')
      CALL MAPINT
      CALL MAPBLA (IAMA)

      call MPLNAM ('Earth..2',4,IAMA) !ioannis, filling

C
C Pre-process the area map.
C
        CALL ARPRAM (IAMA,0,0,0)
C
C
C Tell CONPACK not to do the SET call (since it is already been done)
C and to use mapping function 1 (EZMAP background).
C
      CALL CPSETI ('SET - DO-SET-CALL FLAG',0)
      CALL CPSETI ('MAP - MAPPING FLAG',0)
C
C Tell CONPACK what ranges of X and Y coordinates to send into the
C mapping function.
C
      CALL MAPTRN (RLATMN,RLONMN,XMIN,YMIN)
      CALL MAPTRN (RLATMX,RLONMX,XMAX,YMAX)
C
      CALL CPSETR ('XC1 - X COORDINATE AT I=1',XMIN)
      CALL CPSETR ('XCM - X COORDINATE AT I=M',XMAX)
      CALL CPSETR ('YC1 - Y COORDINATE AT J=1',YMIN)
      CALL CPSETR ('YCN - Y COORDINATE AT J=N',YMAX)
C
C Tell CONPACK exactly what contour levels to use.
C
      CALL CPSETI ('LLP - LINE LEBEL POSITIONING',3)
      CALL CPSETC ('ILT - INFORMATIONAL LABEL TEXT',' ')
      CALL CPSETI ('CLS - CONTOUR LEVEL SELECTOR',0)
      CALL CPSETI ('NCL - NUMBER OF CONTOUR LEVELS',NCONT)
      DO I=1,NCONT
         VALUEI=VALUE(I)
         CALL CPSETI ('PAI - PARAMETER ARRAY IDENTIFIER',I)
         CALL CPSETR ('CLV - CONTOUR LEVEL VALUE',VALUEI)
         CALL CPSETI ('CLU - CONTOUR LEVEL USE',0)  !3
         CALL CPSETI ('LLC - LINE LABEL COLOR INDEX',0) !1
         CALL CPSETI ('LLB - LINE LABEL BOX FLAG',0)
         CALL CPSETI ('LBC - LABEL BOX COLOR INDEX',0)
         CALL CPSETR ('CLL - CONTOUR LEVEL LINE WIDTH',1.)
      END DO
      
      CALL CPSETI('NSD',4)  ! standard value is 4 
      CALL CPSETI ('LIS - LABEL INTERVAL SPECIFIER',1)
      CALL CPSETI ('LLO - LINE LEBEL ORIENTATION',1)
      CALL CPSETR ('LLS - LINE LEBEL SIZE',0.01)
      CALL CPSETI ('LBC - Label Box Color Index ',1)
!      CALL CPSETI ('HLB - High/Low Label Box Flag',0)
!      CALL CPSETI ('HLO - High/Low Label Overlap Flag',1)
!      CALL CPSETI ('HLL - High/Low Line Width Flag',5)
!      CALL CPSETC ('HLT - HIGH/LOW TEXT','H:B:$ZDV$:E:''L:B:$ZDV$:E:')
      CALL CPSETC ('HLT - HIGH/LOW TEXT',' ')
cioannis      CALL CPSETC ('CFT - Constant-field label',' ')
C
      CALL CPRECT (ZDAT,IMT,IMT,JMT,RWRK,LRWK,IWRK,LIWK)
      CALL CPCLAM (ZDAT,RWRK,IWRK,IAMA)
C
      CALL GSFAIS (1)
      CALL ARSCAM (IAMA,XCRA,YCRA,NCRA,IAIA,IGIA,240,COLRAM)


        CALL GSPLCI (1)     
        CALL MAPLBL
        call mapgrd

      call MPLNDR('Earth..2',4)   ! ioannis, coastline

      CALL CPCLDR(ZDAT,RWRK,IWRK)
      CALL CPLBDR (ZDAT,RWRK,IWRK)

cc      CALL MAPDRW
 
      CALL PCSETR ('CS - CONSTANT SPACING',1.25)
      CALL SFSETI ('TYPE OF FILL',0)
      CALL LBLBAR 
CIOANNIS  &(0,XSTA,XSTO,0.02,0.1,NCONT,1.,.3333,LND,0,CLLB,NCONT+1,1)
cioannis  &(0,0.10,0.90,0.08,0.16,NCONT,1.,.3333,LND,0,CLLB,NCONT+1,1)
cc     &(0,0.10,0.90,0.001,0.081,NCONT,1.,.3333,LND,0,CLLB,NCONT+1,1)
cpk     &(0,0.05,0.95,0.021,0.101,NCONT,1.,.3333,LND,0,CLLB,NCONT+1,1)
     &(0,0.05,0.95,0.141,0.241,NCONT,1.,.3333,LND,0,CLLB,NCONT+1,1)

      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
C
      END

C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

      SUBROUTINE COLRAM (XCRA,YCRA,NCRA,IAIA,IGIA,NAGI)
        DIMENSION XCRA(*),YCRA(*),IAIA(*),IGIA(*)
        IAI1= -1
        IAI3=-1
        DO 101 I=1,NAGI
          IF (IGIA(I).EQ. 1) IAI1=IAIA(I)
          IF (IGIA(I).EQ. 3) IAI3=IAIA(I)
  101   CONTINUE

      IF (IAI1.GT.0) THEN
       IF (MAPACI(IAI1) .EQ. 1)  THEN
         CALL GSFACI(94)
         CALL GFA(NCRA-1,XCRA,YCRA)
       ELSE
          CALL GSFACI(93)
          CALL GFA(NCRA-1,XCRA,YCRA)
       ENDIF
      ENDIF
        IF (IAI3.GT.1) THEN     !0
            CALL GSFACI (4+MOD(IAI3,2))
            CALL GFA (NCRA-1,XCRA,YCRA)
            IF (IAI3.GT.1.AND.IAI3.LE.240) THEN
              CALL GSFACI (IAI3+1)
              CALL GFA (NCRA-1,XCRA,YCRA)
            END IF
        END IF
      RETURN
      END
C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

      SUBROUTINE COLRAMOLD (XCRA,YCRA,NCRA,IAIA,IGIA,NAGI)
        DIMENSION XCRA(*),YCRA(*),IAIA(*),IGIA(*)
        IAI1=-1
        IAI3=-1
        DO 101 I=1,NAGI
          IF (IGIA(I).EQ.1) IAI1=IAIA(I)
          IF (IGIA(I).EQ.3) IAI3=IAIA(I)
  101   CONTINUE
        IF (IAI3.GT.0) THEN
            CALL GSFACI (4+MOD(IAI3,2))
            CALL GFA (NCRA-1,XCRA,YCRA)
            IF (IAI3.GE.1.AND.IAI3.LE.240) THEN
              CALL GSFACI (IAI3+1)
              CALL GFA (NCRA-1,XCRA,YCRA)
            END IF
        END IF
      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 COLRCL (XCRA,YCRA,NCRA,IAIA,IGIA,NAGI)
      DIMENSION XCRA(*),YCRA(*),IAIA(*),IGIA(*)
      RETURN
      END

      SUBROUTINE COLRLL (XCRA,YCRA,NCRA,IAIA,IGIA,NAGI)
      DIMENSION XCRA(*),YCRA(*),IAIA(*),IGIA(*)
      IAI1=-1
      DO 101 I=1,NAGI
        IF (IGIA(I).EQ.1.AND.IAIA(I).GT.0) IAI1=IAIA(I)
  101 CONTINUE
      IF (IAI1.GT.0.AND.MAPACI(IAI1).EQ.1) CALL GPL (NCRA,XCRA,YCRA)
      RETURN
      END

C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
      SUBROUTINE ETOAH(HDATE,IME,JME,HDATA,IMA,JMA)
C----------------------------------------------------------------------
C       PROGRAM PANTLJICARA U E-MREZU + E-MREZA U A-MREZU
C       A-MREZA U GEOGRAFSKU OBLAST
C----------------------------------------------------------------------

      DIMENSION HDATA(IMA,JMA)
      DIMENSION HDATE(IME,JME)

C---------------------------------------------------------------------
C     UCITAVANJE PANTLJICARA
C---------------------------------------------------------------------
      DO J=1,JMA
         IF(MOD(J,2).EQ.0) THEN
           DO I=2,IMA-1,2
              HDATA(I,J)=HDATE((I+1)/2 ,J)
           END DO
         ELSE
           DO I=1,IMA,2
              HDATA(I,J)=HDATE((I+1)/2 ,J)
           END DO
         END IF
      END DO

C---------------------------------------------------------------------
C     INTERPOLACIJA E-MREZE U A-MREZU
C---------------------------------------------------------------------

          DO 50  J=2,JMA-1
              MJ=MOD(J,2)
              IF(MJ.EQ.0) THEN
              DO 60 I=3,IMA-2,2
   60               HDATA(I,J)=(HDATA(I-1,J)+HDATA(I+1,J)+HDATA(I,J-1
     &            )+HDATA(I,J+1))*.25
           ELSE
              DO 65 I=2,IMA-1,2
   65               HDATA(I,J)=(HDATA(I-1,J)+HDATA(I+1,J)+HDATA(I,J-1
     &            )+HDATA(I,J+1))*.25
                END IF
   50       CONTINUE


          DO 70 I=2,IMA-1,2
   70               HDATA(I,1)=(HDATA(I-1,1)+HDATA(I+1,1)+
     &            HDATA(I,2))/3
          DO 75 I=2,IMA-1,2
   75               HDATA(I,JMA)=(HDATA(I-1,JMA)+HDATA(I+1,JMA)+
     &            HDATA(I,JMA-1))/3.
          DO 80 J=2,JMA-1,2
   80               HDATA(1,J)=(HDATA(1,J-1)+HDATA(1,J+1)+
     &            HDATA(2,J))/3.
          DO 85 J=2,JMA-1,2
   85               HDATA(IMA,J)=(HDATA(IMA,J-1)+HDATA(IMA,J+1)+
     &            HDATA(IMA-1,J))/3.


       RETURN
       END

C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
      SUBROUTINE ETOAV(VDATE,IME,JME,VDATA,IMA,JMA)
C----------------------------------------------------------------------
C       PROGRAM PANTLJICARA U E-MREZU + E-MREZA U A-MREZU
C       A-MREZA U GEOGRAFSKU OBLAST
C----------------------------------------------------------------------

      DIMENSION   VDATA(IMA,JMA)
      DIMENSION   VDATE(IME,JME)

C---------------------------------------------------------------------
C     UCITAVANJE PANTLJICARA
C---------------------------------------------------------------------
      DO J=1,JMA
         IF(MOD(J,2).EQ.0) THEN
           DO I=1,IMA,2
              VDATA(I,J)=VDATE((I+1)/2 ,J)
           END DO
         ELSE
           DO I=2,IMA-1,2
              VDATA(I,J)=VDATE((I+1)/2 ,J)
           END DO
         END IF
      END DO

C---------------------------------------------------------------------
C     INTERPOLACIJA E-MREZE U A-MREZU
C---------------------------------------------------------------------

          DO 50  J=2,JMA-1
              IF(MOD(J,2).EQ.0) THEN
              DO 60 I=2,IMA-1,2
   60               VDATA(I,J)=(VDATA(I-1,J)+VDATA(I+1,J)+VDATA(I,J-1
     &            )+VDATA(I,J+1))*.25
           ELSE
              DO 65 I=3,IMA-2,2
   65               VDATA(I,J)=(VDATA(I-1,J)+VDATA(I+1,J)+VDATA(I,J-1
     &            )+VDATA(I,J+1))*.25
                END IF
   50       CONTINUE


          DO 70 I=3,IMA-2,2
   70               VDATA(I,1)=(VDATA(I-1,1)+VDATA(I+1,1)+
     &            VDATA(I,2))/3
          DO 75 I=3,IMA-2,2
   75               VDATA(I,JMA)=(VDATA(I-1,JMA)+VDATA(I+1,JMA)+
     &            VDATA(I,JMA-1))/3.
          DO 80 J=3,JMA-2,2
   80               VDATA(1,J)=(VDATA(1,J-1)+VDATA(1,J+1)+
     &            VDATA(2,J))/3.
          DO 85 J=3,JMA-2,2
   85               VDATA(IMA,J)=(VDATA(IMA,J-1)+VDATA(IMA,J+1)+
     &            VDATA(IMA-1,J))/3.
          VDATA(1,1)=0.5*(VDATA(2,1)+VDATA(1,2))
          VDATA(IMA,1)=0.5*(VDATA(IMA-1,1)+VDATA(IMA,2))
          VDATA(IMA,JMA)=0.5*(VDATA(IMA-1,JMA)+VDATA(IMA,JMA-1))
          VDATA(1,JMA)=0.5*(VDATA(1,JMA-1)+VDATA(2,JMA))



       RETURN
       END
C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
      SUBROUTINE DFCLRS_INC(N)
      parameter (nc=52)
      dimension R(nc), G(nc), B(nc)
C
cc      OPEN(UNIT=10,FILE='../src/data1.dat',
cc     & ACCESS='SEQUENTIAL',STATUS='OLD',FORM='FORMATTED')
      CALL GSCR (1,0,1.,1.,1.)
      CALL GSCR (1,1,0.,0.,0.)

cc       DO I = 1, nc
cc           READ(10,*) R(i),G(i),B(i)
cc        ENDDO

c        incr = ifix(nc/float(N) + 0.5)
c        ii = 0
c        DO I = 2, nc, incr
c           ii = ii + 1
c           if(ii.eq.2)then
c              r(i)=1
c              g(i)=1
c              b(i)=1
c           endif
c           CALL GSCR (1,II,R(i),G(i),B(i))
c        ENDDO

      CALL GSCR (1,0,1.,1.,1.)
      CALL GSCR (1,1,0.,0.,0.)
      CALL GSCR (1,93,0.65,0.46,0.46)
      CALL GSCR (1,93,1.,1.,1.)
      CALL GSCR (1,94,1.,1.,1.)

      CALL GSCR (1,2,1.  ,1.  ,1.)

c ******* colours for concentrations *******
c      CALL GSCR (1,3,0.  ,0.  ,0.82)   !  1
c      CALL GSCR (1,4,0.  ,0.3 ,1. )    !  2
c      CALL GSCR (1,5,0.  ,0.75,1. )    !  3
c      CALL GSCR (1,6,0.  ,1.  ,1. )    !  4
c      CALL GSCR (1,7,0.  ,1.  ,0.85)   !  5
c      CALL GSCR (1,8,0.  ,1.  ,0.)     !  6
cc      CALL GSCR (1,9,0.  ,0.9  ,0.)    !  7c
cc      CALL GSCR (1,10,0.  ,0.8  ,0.)   !  8c
c      CALL GSCR (1,9,0.  ,0.5  ,0.)   !  9

c      CALL GSCR (1,10,1.  ,1.  ,0.)    !  10
c      CALL GSCR (1,11,1.  ,0.9 ,0.)    !  11
c      CALL GSCR (1,12,1.  ,0.7 ,0.)    !  12 
c      CALL GSCR (1,13,1.  ,0.  ,0.75)  !  13
c      CALL GSCR (1,14,0.9 ,0.  ,0.)    !  14
cc      CALL GSCR (1,17,1.  ,0.  ,1.)    !  15c

c ****** NEW SCALE --- ANTIGONI --- *********
      CALL GSCR (1,3, 0.,1.,1.)        !  1
      CALL GSCR (1,4,0.  ,0.85  ,1. )  !  2
      CALL GSCR (1,5,0.  ,0.75,1. )    !  3
      CALL GSCR (1,6,0.  ,0.3 ,1. )    !  4
      CALL GSCR (1,7,0.  ,1.  ,0.)     !  5
      CALL GSCR (1,8,0.  ,0.7  ,0.)    !  6
      CALL GSCR (1,9,0.  ,0.5  ,0.)    !  7
      CALL GSCR (1,10,0., 0.3  ,0.)    !  8
      CALL GSCR (1,11,1.  ,1.  ,0.)    !  9
      CALL GSCR (1,12,1.  ,0.8 ,0.)    ! 10
      CALL GSCR (1,13,1.  ,0.6 ,0.)    ! 11
      CALL GSCR (1,14,1.  ,0.4 ,0.)    ! 12
      CALL GSCR (1,15,1. ,0.  ,0.)     ! 13
      CALL GSCR (1,16,0.7 ,0.  ,0.)    ! 14
      CALL GSCR (1,17,1.  ,0.  ,0.75)  ! 15


        RETURN
	END

C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
      SUBROUTINE DFCLRS(N)
C
      OPEN(UNIT=10,FILE='../src/data1.dat',
     & ACCESS='SEQUENTIAL',STATUS='OLD',FORM='FORMATTED')

      CWCOMP = 1.4
      CWROT = 1.0
      RFACT = 0.5 *CWROT
C
ccc        CALL GSCR (1,0,0.,0.,0.)
ccc        CALL GSCR (1,1,1.,1.,1.)
      CALL GSCR (1,0,1.,1.,1.)
      CALL GSCR (1,1,0.,0.,0.)

C
ccccc        DO 101 I=2,N
        DO 101 I=2,N
cccc        DO 101 I=2,N+1
      READ(10,*) R,G,B
cccc      print*,R,G,B
cccc        S=FLOAT(I-1)/FLOAT(N-1)
cccc        T= CWCOMP*(S-RFACT)
cccc         R=128.0 + 127.0 * ATAN( 7.0*T ) / 1.57
cccc         G=128.0 + 127.0 * (2 * EXP(-7*T*T) - 1)
cccc         B=128.0 + 127.0 * ATAN( -7.0*T ) / 1.57
cccc         R=R/255
cccc         G=G/255
cccc         B=B/255
        if(i.eq.2)then
        r=1
        g=1
        b=1
        endif
          CALL GSCR (1,I,R,G,B)
  101   CONTINUE
        RETURN
C
      END
C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
                         SUBROUTINE CONH12(HI1,HI2)
C     ******************************************************************
C     *                                                                *
C     *  ROUTINE FOR REORDERING THIBUE-HEIGHT-POINT-1-DIMENSIONAL      *
C     *  MATRICES FOR 2-DIMENSIONAL INDEXING                           *
C     *                                                                *
C     ******************************************************************
#include "all.inc"

c      parameter(TLM0D=15.0,TPH0D=39.5)
c      parameter(WBD=-24.6,SBD=-11.5,DLMD=0.1,DPHD=0.1,DTB=30.,LM=38)
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
                             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
     & HI2    (IM,JM),HI1    (IMJM)
C-----------------------------------------------------------------------
      K=0
              DO 100 J=1,JM
          DO 100 I=1,IM-1+MOD(J,2)
      K=K+1
  100  HI2(I,J)=HI1(K)
              DO 101 J=2,JM-1,2
  101  HI2(IM,J)=HI2(IM-1,J)
C-----------------------------------------------------------------------
                             RETURN
                             END 
C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
                         SUBROUTINE CONH21(H2,H1B)
C     ******************************************************************
C     *                                                                *
C     *  ROUTINE FOR REORDERING THIBUE-HEIGHT-POINT-TWO-DIMENSIONAL    *
C     *  MATRICES FOR ONE-DIMENIONAL INDEXING                          *
C     *                                                                *
C     ******************************************************************
#include "all.inc"

c      parameter(TLM0D=15.0,TPH0D=39.5)
c      parameter(WBD=-24.6,SBD=-11.5,DLMD=0.1,DPHD=0.1,DTB=30.,LM=38)
C-----------------------------------------------------------------------
                             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)

      parameter(IMT=2*IM-1,JMT=JM/2+1)
C-----------------------------------------------------------------------
                             D I M E N S I O N
     & H2    (IMT,JMT),H1B   (IMJM)
C-----------------------------------------------------------------------
      K=0
              DO 100 J=1,JMT
          DO 110 I=1,IMT,2
      K=K+1
  110  H1B(K)=H2(I,J)
C
      IF(J.EQ.JMT)           GO TO 100
C
          DO 120 I=2,IMT,2
      K=K+1
  120  H1B(K)=H2(I,J)
C
  100          CONTINUE
C-----------------------------------------------------------------------
      return
      end
