      program flexpartplot_latlon


      parameter(ierrf=6,iworkmeta=1,iworkx=8,iworkps=20)
      parameter(lunitx=80,lunitps=81,lunitmeta=82)
      parameter(idmeta=1,idx=2,idps=3,maxcon=100)
      parameter(lrwk=29000,liwk=29000,lmap=9000000)
      parameter(nwrk=29000,isiz=7)
      real xwrk(nwrk),ywrk(nwrk),rwrk(lrwk),iwrk(liwk)
      integer map(lmap),iarea(isiz),igrp(isiz)

      real zlev(0:maxcon)
      character lbls(0:maxcon)*4,label*64
      character concchar*100,latchar1*6,agechar*20,specchar*10
      character concchar2*100,latchar2*6,latchar3*6,latchar4*6
      character headlonlat*9,unitlabel*20,xlabel*9
      character headlonlat2*9
      character alspec*10,aiage*6

      parameter(nxmax=420,nymax=180,nzmax=29)
      parameter(maxpoint=800000,maxspec=4,maxageclass=15)


      real outheight(nzmax),heightnn(nxmax,nymax,0:nzmax)
      integer ireleasestart(maxpoint),ireleaseend(maxpoint)
      integer itime,iiday
      integer npart(maxpoint),kind(maxpoint),lage(0:maxageclass)
      real xpoint(maxpoint),ypoint(maxpoint),zpoint1(maxpoint)
      real zpoint2(maxpoint),xmass(maxpoint,maxspec)
      real drygrid(nxmax,nymax,maxspec,maxageclass)
      real wetgrid(nxmax,nymax,maxspec,maxageclass)
      real grid(nxmax,nymax,nzmax,maxspec,maxageclass)
      real gridtotal(nxmax,nymax,nzmax,maxspec)
      parameter(nregrid=10,maxpoints=150)
      real zgrid(maxpoints,nregrid*nzmax)
      real oro(nxmax,nymax),ihour,iday
      character dirname*110,datesname*110,maxvalue*10,minvalue*10
      character meanvalue*10,compoint(maxpoint)*10,species(maxspec)*10
      character filename*120,aday*8,atime*6,charlocation*7
      character charkindsec*6

      real yra(5),xra(5)


      external cpdrpl,sfill_grau,fill



C Read CONTROL file
*******************

      open(10,file='CONTROL')
      read(10,'(a110)') dirname
      read(10,'(a110)') datesname
      read(10,'(a7)') charlocation
      read(10,'(a6)') charkindsec
      read(10,*) lconcormixratio
      read(10,*) decayconstant
      read(10,*) longflag            ! if negative, zonal mean is plotted
      read(10,*) xll,yll
      read(10,*) xur,yur
c     read(10,*) kindflag            ! if 1, against lat., if not, against long.
      read(10,'(a10)') alspec               ! number of species to be plotted
      read(10,*) lageplot            ! age classes in s to be plotted, 0 if all are to be plotted
      read(10,*) lcolor              ! color or B&W
      read(10,*) loglin              ! linear or logarithmic scale?
      read(10,*) scaleconc           ! scaling factor for conc./mixing ratios
      read(10,*) overrulemax         ! maximum value for conc./mixing ratios
      read(10,'(a20)') unitlabel      ! label for concentration unit
      read(10,*) numcon
      read(10,*) ifact               ! number of values reported on labelbar
      read(10,*) ltickx              ! number of ticks on x axis
      read(10,*) lticky              ! number of ticks on y axis
      read(10,*) heightmax
      read(10,*) iswitch
      close(10)
      if (numcon.gt.maxcon) numcon=maxcon
      ifact=numcon/(ifact-1)

      len=index(dirname,' ')-1


c plot against latitude or longitude
************************************

      if (xll.eq.xur) then
       kindflag = 1
      else
       kindflag = 0
      endif



C Open GKS, open and activate workstations
******************************************

      call gopks(ierrf,0)
      call gopwk(idmeta,lunitmeta,iworkmeta)
      call gacwk(idmeta)

      call pcseti('FN',22)
      call gaseti('LTY',1)
      call setusv('LW',1000)
      CALL PCSETI('CC',1)



C Create the color table
************************

      xc=1.0
      xb=0.0
      xg=1.0
      xl=0.0
      xh=1.0
      jr=0
      jn=0
      do 10 i=1,numcon+1
        call colspec(lcolor,i-1,numcon+1,xl,xh,jr,jn,xc,xb,xg,0,i+1)
10      continue


C Read model results
********************

      if (iswitch.eq.0) then
       filename=dirname(1:len)//'header_nest'
      else
       filename=dirname(1:len)//'header'
      endif
      open(10,file=filename,form='unformatted',status='old')
      read(10) ibdate,ibtime
      read(10) loutstep,loutaver,loutsample
      read(10) outlon0,outlat0,numxgrid,numygrid,dxout,dyout
c           write(*,*) numxgrid,numygrid,dxout,dyout
      read(10) numzgrid,(outheight(i),i=1,numzgrid)
      read(10) jjjjmmdd,ihmmss
      read(10) nspec
      nspec=nspec/3
      do 8 n=1,nspec
        read(10) numzgrid,species(n)
        read(10) numzgrid,species(n)
8       read(10) numzgrid,species(n)


      read(10) numpoint
      do 13 i=1,numpoint
        read(10) ireleasestart(i),ireleaseend(i)
        read(10) xpoint(i),ypoint(i),xp2,yp2,zpoint1(i),zpoint2(i)
        read(10) npart(i),kind(i)
        read(10) compoint(i)
        do 13 j=1,nspec
          read(10)
          read(10)
13        read(10) xmass(i,j)
      read(10) method
      read(10) nageclass,(lage(i),i=1,nageclass)
      lage(0)=0
c         write(*,*) (lage(i),i=0,nageclass)


      do 130 ix=1,numxgrid
 130     read(10) (oro(ix,jy),jy=1,numygrid)

c     do ix=1,numxgrid
c     do jy=1,numygrid
c       if (oro(ix,jy).gt.0.) write(*,*) oro(ix,jy)
c     end do
c     end do


c         write(*,*) (outheight(i),i=1,numzgrid)

      if (iswitch.eq.0) then
        if (lconcormixratio.eq.1) then
          filename=dirname(1:len)//'grid_conc_nest'
        else if (lconcormixratio.eq.2) then
          filename=dirname(1:len)//'grid_pptv_nest'
        else
          filename=dirname(1:len)//'grid_time_nest'
        endif
      else
        if (lconcormixratio.eq.1) then
          filename=dirname(1:len)//'grid_conc'
        else if (lconcormixratio.eq.2) then
          filename=dirname(1:len)//'grid_pptv'
        else
          filename=dirname(1:len)//'grid_time'
        endif
      endif

c determine up to which ageclass to be summed up and plotted
************************************************************

      if (lageplot.eq.0) then
        iageclass = nageclass
        goto 888
      endif

      do nage=1,nageclass
       if (lageplot.eq.lage(nage)) then
        iageclass=nage
        goto 888
       endif
      end do
      stop 'selected ageclass does not exist'
 888  continue


c determine which species to be plotted
***************************************

      do n=1,nspec
       if (alspec.eq.species(n)) then
        lspec=n
        goto 777
       endif
      end do
      stop 'selected species alspec does not exist'
 777  continue

c check whether xll, yll, xur, and yur are within domain boundaries
*******************************************************************

       if (yll.lt.outlat0) yll=outlat0
       if (yur.lt.outlat0) yur=outlat0 + dyout
       if (yll.gt.outlat0 + float(numygrid)*dyout) 
     +     yll=outlat0 + float(numygrid)*dyout - dyout
       if (yur.gt.outlat0 + float(numygrid)*dyout) 
     +     yur=outlat0 + float(numygrid)*dyout
       if (xll.lt.outlon0) xll=outlon0
       if (xur.lt.outlon0) xur=outlon0 + dxout
       if (xll.gt.outlon0 + float(numxgrid)*dxout) 
     +     xll=outlon0 + float(numxgrid)*dxout - dxout
       if (xur.gt.outlon0 + float(numxgrid)*dxout) 
     +     xur=outlon0 + float(numxgrid)*dxout

C Determine left and right coordinates of plot
**********************************************

      if (longflag.ne.1) then
        yll=outlat0
        yur=outlat0+float(numygrid)*dyout
        kindflag=1
      endif

      if (kindflag.eq.1) then
        xmin=yll
        xmax=yur
        xlabel='Latitude'
ccc     write(latchar,'(i5)') nint((xll+xur)/2.)
        write(latchar1,'(f6.1)') xll
        write(latchar2,'(f6.1)') xur
        write(latchar3,'(f6.1)') yll
        write(latchar4,'(f6.1)') yur
        headlonlat='LONGITUDE'
        headlonlat2='LATITUDE'
      else
        xmin=xll
        xmax=xur
        xlabel='longitude'
ccc     write(latchar,'(i5)') nint((yll+yur)/2.)
        write(latchar1,'(f6.1)') yll
        write(latchar2,'(f6.1)') yur
        write(latchar3,'(f6.1)') xll
        write(latchar4,'(f6.1)') xur
        headlonlat='LATITUDE'
        headlonlat2='LONGITUDE'
      endif


c Calculate height, which is outheight plus topography
******************************************************

      do 150 ix=1,numxgrid
        do 150 jy=1,numygrid
          heightnn (ix,jy,0) = oro(ix,jy)
          do 150 i=1,numzgrid
            heightnn (ix,jy,i) = outheight(i) + oro(ix,jy)
150         continue         



C Read data in sparse matrix format
***********************************

C Loop about all times
**********************

      open(20,file=datesname,form='formatted',status='old')
      it=0
100   read(20,'(i8,i6)',end=99) inday,intime
      write(aday,'(i8.8)') inday
      write(atime,'(i6.6)') intime
      if (iswitch.eq.0) then
         open(10,file=filename(1:len+14)//'_'//aday//atime,
     +form='unformatted',status='old')
      else
         open(10,file=filename(1:len+9)//'_'//aday//atime,
     +form='unformatted',status='old')
      endif
      it=it+1
      read(10,end=99) itime
c     write(*,*) it,itime

C Loop about all species
************************

      do 38 k=1,nspec
        if ((lspec.gt.0).and.(k.gt.lspec)) goto 48


C Initialize total fields
*************************

        do 62 ix=1,numxgrid
          do 62 jy=1,numygrid
            do 62 kz=1,numzgrid
62            gridtotal(ix,jy,kz,k)=0.

C Loop about all age classes
****************************

        do 38 nage=1,nageclass

C Initialize age class fields
*****************************

        do 63 ix=1,numxgrid
          do 63 jy=1,numygrid
            wetgrid(ix,jy,k,nage)=0.
            drygrid(ix,jy,k,nage)=0.
            do 63 kz=1,numzgrid
63            grid(ix,jy,kz,k,nage)=0.


C Read wet deposition
*********************

          read(10) lval
          if (lval.eq.1) then           !sparse matrix
33          read(10,end=99) n,xx
            if (n.lt.0) goto 32
            jy=n/numxgrid
            ix=n-numxgrid*jy
            ix=ix+1
            jy=jy+1
            wetgrid(ix,jy,k,nage)=xx
            goto 33
          else
            do 30 ix=1,numxgrid
30            read(10) (wetgrid(ix,jy,k,nage),jy=1,numygrid)
          endif

32        continue

C Read dry deposition
*********************

          read(10) lval
          if (lval.eq.1) then           !sparse matrix
36          read(10,end=99) n,xx
            if (n.lt.0) goto 35
            jy=n/numxgrid
            ix=n-numxgrid*jy
            ix=ix+1
            jy=jy+1
            drygrid(ix,jy,k,nage)=xx
            goto 36
          else
            do 31 ix=1,numxgrid
31            read(10) (drygrid(ix,jy,k,nage),jy=1,numygrid)
          endif

35        continue

C Read concentrations
*********************

          read(10) lval
          if (lval.eq.1) then           !sparse matrix
40          read(10,end=99) n,xx
            if (n.lt.0) goto 78
            kz=n/(numxgrid*numygrid)
            jy=(n-kz*numxgrid*numygrid)/numxgrid
            ix=n-numxgrid*numygrid*kz-numxgrid*jy
            ix=ix+1
            jy=jy+1
            grid(ix,jy,kz,k,nage)=xx
            goto 40
          else
            do 39 kz=1,numzgrid
              do 39 ix=1,numxgrid
39              read(10) (grid(ix,jy,kz,k,nage),jy=1,numygrid)
          endif

78        continue

C End species loop, end age class loop
**************************************

38      continue

48    continue


C Sum up all age classes to total fields and scale values as defined by user
****************************************************************************

          do 41 k=1,nspec
           if ((lspec.gt.0).and.(k.gt.lspec)) goto 49
          do 41 nage=1,iageclass
          do 41 ix=1,numxgrid
            do 41 jy=1,numygrid
              do 41 kz=1,numzgrid
                grid(ix,jy,kz,k,nage)=grid(ix,jy,kz,k,nage)*
     +          exp(-1.*float(lage(nage)+lage(nage-1))/2./
     +          decayconstant)
                grid(ix,jy,kz,k,nage)=grid(ix,jy,kz,k,nage)*scaleconc
                gridtotal(ix,jy,kz,k)=gridtotal(ix,jy,kz,k)+
     +          grid(ix,jy,kz,k,nage)
                grid(ix,jy,kz,k,nage)=grid(ix,jy,kz,k,nage)/
     +          float(lage(nage)-lage(nage-1))*86400.
41              continue

 49      continue



c     if (it.gt.8) then
c       if (mod(it,2).eq.1) goto 3333
c     endif


C Begin species loop, begin age class loop
******************************************

      do 200 k=1,nspec
        specchar=alspec

      if ((lspec.gt.0).and.(k.ne.lspec)) goto 200

         if (lage(iageclass).lt.86400) then
          rage=float(lage(iageclass))/3600.
          write(aiage,'(f6.2)') rage
          agechar=' 0 - '//aiage//' HRS'
         else
          rage=float(lage(iageclass))/86400.
          write(aiage,'(f6.2)') rage
          agechar=' 0 - '//aiage//' DAYS'
         endif

         if (lage(iageclass).eq.999999999) then
          agechar=' all'
         endif


          do 87 ix=1,numxgrid
            do 87 jy=1,numygrid
              do 87 kz=1,numzgrid
87              grid(ix,jy,kz,k,1)=gridtotal(ix,jy,kz,k)


C Define contour field for that time step
*****************************************

        fact=outheight(numzgrid)/float(nregrid*numzgrid)
        zmax=0.
        zmin=1.e30
        zmean=0.
        numbmean=0
        do 202 kzn=1,nregrid*numzgrid
          altitude=(float(kzn)-0.5)*fact
          if (longflag.gt.0) then
            npoints=maxpoints
            do 230 i=1,npoints
              xx=xll+float(i-1)/float(npoints-1)*(xur-xll)
              yy=yll+float(i-1)/float(npoints-1)*(yur-yll)
              indx=(xx-outlon0)/dxout+1
              indy=(yy-outlat0)/dyout+1
              do 220 kz=1,numzgrid
                if (heightnn(indx,indy,kz).gt.altitude) goto 221
220             continue
221           kz=min(kz,numzgrid)
              zgrid(i,kzn)=grid(indx,indy,kz,k,1)
              if (heightnn(indx,indy,0).gt.altitude)
     +        zgrid(i,kzn)=1.e30  ! below topo
230           continue
          else       ! calculate zonal average
            npoints=numygrid
	    do 201 jy=1,numygrid
              zgrid(jy,kzn)=0.
              numb=0
              do 204 ix=1,numxgrid
                do 222 kz=1,numzgrid
                  if (heightnn(ix,jy,kz).gt.altitude) goto 223
222               continue
223             kz=min(kz,numzgrid)
                if (heightnn(ix,jy,0).lt.altitude) then ! above topo
                  numb=numb+1
                  zgrid(jy,kzn)=zgrid(jy,kzn)+grid(ix,jy,kz,k,1)
                endif
204             continue
              if (numb.gt.0) then
                zgrid(jy,kzn)=zgrid(jy,kzn)/float(numb)
              else
                zgrid(jy,kzn)=1.e30  ! all below topo
              endif
201           continue
          endif
          do 203 jy=1,npoints
            if ((zgrid(jy,kzn).lt.9.e29).and.(zgrid(jy,kzn).gt.zmax))
     +      zmax=zgrid(jy,kzn)
            if (zgrid(jy,kzn).lt.zmin) zmin=zgrid(jy,kzn)
            if (zgrid(jy,kzn).lt.9.e29) then
              zmean=zmean+zgrid(jy,kzn)
              numbmean=numbmean+1
            endif
203         continue
202       continue
        zmean=zmean/float(numbmean)

      write(maxvalue,'(e10.3)') zmax
      write(minvalue,'(e10.3)') zmin
      write(meanvalue,'(e10.3)') zmean



      call gsplci(1)
      call gsln(1)
      call gslwsc(1.)


C Map the plotter frame to user coordinates
*******************************************

      call set(0.16,0.91,0.28,0.83,xmin,xmax,0.,
     +heightmax,1)
      call cpseti('SET',0)

C Initialize conpack
********************

      call cpsetr('XC1',xmin)
      call cpsetr('XCM',xmax)
      call cpsetr('YC1',0.)
      call cpsetr('YCN',outheight(numzgrid)/1000.)


C Make logarithmic or linear scale
**********************************

      if (overrulemax.gt.0.) zmax=overrulemax

      call cpseti('CLS',0)
      call cpseti('NCL',numcon+1)
      if (loglin.eq.1) then
        zmax=float(int(2.*log10(zmax)))/2.
        zlev(0)=10**(zmax-float(numcon)/2.)
        do 28 i=1,numcon
          call cpseti('PAI - Parameter Array Index',i)
          zlev(i)=10**(zmax-float(numcon-i)/2.)
          call cpsetr('CLV - CONTOUR LEVEL VALUES',zlev(i))
28        call cpseti('CLU - Contour Level Use Flag',2)
      else
        zlev(0)=zmin
        do 27 i=1,numcon
          call cpseti('PAI - Parameter Array Index',i)
          zlev(i)=zmax-float(numcon-i)/float(numcon)*zmax
          call cpsetr('CLV - CONTOUR LEVEL VALUES',zlev(i))
27        call cpseti('CLU - Contour Level Use Flag',2)
      endif


C Avoid nasty colors if maximum pre-set value is exceeded
*********************************************************

      do 37 i=1,npoints
        do 37 kz=1,nregrid*numzgrid
37        if ((zgrid(i,kz).gt.zmax).and.
     +    (zgrid(i,kz).lt.9.e29)) zgrid(i,kz)=zmax


      call cpsetr('SPV',1.e30)
      call cprect(zgrid,maxpoints,npoints,nregrid*numzgrid,rwrk,lrwk,
     +iwrk,liwk)
      call CPPKCL(zgrid,RWRK,IWRK)
      call CPPKLB(zgrid,RWRK,IWRK)


      CALL ARINAM(MAP, LMAP)
      CALL CPCLAM(zgrid, RWRK, IWRK, MAP)
      CALL ARSCAM(MAP, XWRK, YWRK, NWRK, IAREA, IGRP, ISIZ, FILL)

      call labmod('(F6.1)','(I2)',0,0,15,15,20,20,0)
      call gridal(ltickx,1,lticky,1,1,1,5,0.,0.)


      do 34 i=1,npoints
        if (kindflag.eq.1) then
          yy=yll+float(i-1)/float(npoints-1)*(yur-yll)
        else
          yy=xll+float(i-1)/float(npoints-1)*(xur-xll)
        endif
        do 34 kz=1,nregrid*numzgrid
          if (zgrid(i,kz).gt.9.e29) then     ! mark regions below topography
            call points(yy,(float(kz)-.5)/float(nregrid*numzgrid)*
     +      outheight(numzgrid)/1000.,1,-2,0)
          endif
34        continue


C Draw the label bar
********************

      if (loglin.eq.1) then
        do 26 i=0,numcon
26      write(lbls(i),'(i4)') nint(10**(zmax-float(numcon-1-i)/2.))
      else
        do 24 i=0,numcon
24        write(lbls(i),'(i4)') nint(zlev(i))
      endif


C Number of text labels


      call set(0.,1.,0.,1.,0.,1.,0.,1.,1)

      xlic1=0.1
      xlic2=0.9
      ylic1=0.165
      ylic2=0.190
      do i=1,numcon
        xxo=xlic1+(xlic2-xlic1)*float(i-1)/float(numcon)
        xxpp=xlic1+(xlic2-xlic1)*float(i)/float(numcon)
        xra(1)=xxo
        xra(2)=xxo
        xra(3)=xxpp
        xra(4)=xxpp
        yra(1)=ylic1
        yra(2)=ylic2
        yra(3)=ylic2
        yra(4)=ylic1
        call sfsgfa(xra,yra,4,dst,nst,ind,nnd,i+2)
        if (jblack.eq.1) then
        call setusv('LW',1000)
        call curved(xra,yra,2)
        call setusv('LW',2000)
        end if
      end do

      do i=0,numcon
      xlp=xlic1+(xlic2-xlic1)*float(i)/float(numcon)
      if (i.eq.0) then
        call plchhq(xlp-0.020,0.15,lbls(i),.015,0.,0.)
      else if (i.eq.numcon) then
        call plchhq(xlp-0.010,0.15,lbls(i),.015,0.,0.)
      else if (mod(i,ifact).eq.0) then
        call plchhq(xlp,0.15,lbls(i),.015,0.,0.)
      endif
      end do
      call setusv('LW',2000)
      do i=1,numcon-1
      if (mod(i,ifact).eq.0) then
        xlp=xlic1+(xlic2-xlic1)*float(i)/float(numcon)
        call frstd(xlp,ylic1)
        call vectd(xlp,ylic2)
      endif
      end do
      xra(1)=xlic2
      xra(2)=xlic2
      xra(3)=xlic1
      xra(4)=xlic1
      xra(5)=xra(1)
      yra(1)=ylic1
      yra(2)=ylic2
      yra(3)=ylic2
      yra(4)=ylic1
      yra(5)=yra(1)
      call curved(xra,yra,5)

      call pcseti('FN',21)
      call plchhq(.5,.120,unitlabel,.016,0.,0.)

      call setusv('LW',1000)


C Draw the release date and time and current time
*************************************************
      iday=(itime/3600)/24
      iiday=nint(iday)
      ihour=(itime/3600)
      iihour=nint(ihour)
      write(label,900) ibdate,ibtime-30000,inday,intime
900   format('LATEST ANALYSIS   ',i8,' ',i6,' UTC   ACTUAL TIME  ',
     +i8,' ',i6,' UTC')
      call set(0.,1.,0.,1.,0.,1.,0.,1.,1)
      if (longflag.gt.0) then
ccc     concchar='MIXING RATIO OF '//specchar//' AT '
ccc  +  //latchar//' DEG '//headlonlat//' FOR AGE CLASS '//agechar
      if (lconcormixratio.eq.2) then
        concchar=specchar//' MIXING RATIO FOR AGE CLASS '//agechar
      else
        concchar=specchar//' CONCENTRATION FOR AGE CLASS '//agechar
      endif
ccc     concchar2='           '//charkindsec//' CROSS SECTION FOR '
ccc  +//charlocation//'        BASED ON ECMWF ANALYSIS'
      else
        concchar='ZONALLY AVERAGED MIXING RATIO OF '//specchar//
     +  ' FOR AGE CLASS '//agechar
      endif

      concchar2='CROSS SECTION FROM '//latchar1//' TO '//latchar2//' '
     +//headlonlat//' AND '//latchar3//' TO '//latchar4//
     +' '//headlonlat2 

ccc   call plchhq(.5,.95,concchar,.013,0.,0.)
ccc   call plchhq(.5,.91,concchar2,.013,0.,0.)
      call plchhq(.6,.92,concchar2,.013,0.,0.)
      call plchhq(.6,.89,concchar,.013,0.,0.)
      call plchhq(.5,.864,label,.013,0.,0.)

c     call plchhq(.5,.92,'Mean value: '//meanvalue,.013,0.,0.)
c     call plchhq(.5,.892,'Maximum value: '//maxvalue,.013,0.,0.)
c     call plchhq(.5,.864,'Minimum value: '//minvalue,.013,0.,0.)

      call plchhq(.09,.5,'Altitude [km]',.014,90.,-1.)
      call plchhq(.54,.215,xlabel,.014,0.,0.)


      call frame

200   continue

3333  continue

      close(10)
      goto 100

99    close(10)

C Deactivate and close workstations, close GKS
**********************************************

      call gdawk(idmeta)
      call gclwk(idmeta)
      call gclks


      end



       SUBROUTINE SFILL_grau (XWRK, YWRK, NWRK, IAREA, IGRP, NGRPS)

       REAL XWRK(*), YWRK(*), RSCR(5000)
       INTEGER IAREA(*), IGRP(*), ISCR(5000)
       integer i , iarea3

       DO 10, I=1, NGRPS
          IF (IGRP(I) .EQ. 3) IAREA3 = IAREA(I)
   10  CONTINUE
       IF (IAREA3 .GE. 2) THEN
c          if (iarea3 .ge. 6) write (*,*) iarea3
           CALL SFSGFA (XWRK, YWRK, NWRK, RSCR, 5000, ISCR, 
     +                  5000, iarea3+1)
       ENDIF
       RETURN
       END



      subroutine colspec(idp,kc,nc,xl,xh,jr,jn,xc,xb,xg,jl,nnn)
C Sets a colour for NCAR-Graphics applications
C based on a flexible smooth colour spectrum generator
C
C Used for producing colour spectra with almost continuous
C colour gradients
C
C Version 1.1, Written by Paul James, May 2001
C
C Options / Key
C
C IDP = 0 (greyscale), otherwise colourscale
C KC is current index (eg. current contour)
C      must be a fraction of NC and not less than 1
C NC is total number of contours / lines (for example)
C XL (0.,1.) is the desired lowest colour index within
C      the full spectrum range f=0-1
C XH (0.,1.) is the desired highest colour index within
C      the full spectrum range f=0-1
C JR allow a colour rotation to be made
C      0 = RGB (default), 1 = RBG, 2 = GBR
C      3 = GRB, 4 = BRG, 5 = BGR
C JN reverses the direction of the colour scale if < > 0
C XC reduces the colour contrast between 0. and 1. (default)
C XB increases (0.,1.) or decreases (0.,-1.) the brightness
C XG is a post-spectral gamma adjustment of the RGB curves
C     > 1. for a brightness curve skewed to darker shades
C          (lighter parts of spectrum more contrasty)
C     < 1. for a brightness curve skewed to lighter shades
C          (darker parts of spectrum more contrasty)
C      XGA must be within the range (0.2,5.0)
C JL applies a linear brightness gradient to the colour scale
C    so that one end fades to white or black
C JL=0 (default) means no gradient is applied
C JL=1 white low end    >=2 white high end
C JL=-1 black low end    <=-2 black high end
C 
C False parameter settings are corrected to the defaults:
C (2,kc,nc,0.,1.,0,0,1.,0.,1.,0)
C 
C The default colour spectrum has the following key points:
C
C f = 0.000 WHITE
C f = 0.125 LILAC-MAGENTA
C f = 0.250 BLUE
C f = 0.375 CYAN
C f = 0.500 GREEN
C f = 0.625 YELLOW
C f = 0.700 ORANGE
C f = 0.775 RED
C f = 0.875 PINK-MAGENTA
C f = 1.000 BLACK
C
      real r,g,b,f

      call gscr(1,0,1.,1.,1.)         ! Background = white
      call gscr(1,1,0.,0.,0.)         ! Text = black
      yf=float(kc-1)/float(nc-1)
      if (yf.gt.1.) yf=1.
      if (yf.lt.0.) yf=0.
C Check xl and xh range for errors
C On error re-apply default settings (0.,1.)
      if (xh.gt.1..or.xl.gt.1.) then
         xl=0.
         xh=1.
      end if
      if (xh.lt.0..or.xl.lt.0.) then
         xl=0.
         xh=1.
      end if
      if (xg.lt.0.2) xg=0.2
      if (xg.gt.5.) xg=5.
C Convert location in desired colour range
C to real colour fraction
      f=xl+yf*(xh-xl)
      if (f.gt.1.) f=1.
      if (f.lt.0.) f=0.
      if (jn.ne.0) then
         f=1.-f
      end if
      if (idp.eq.0) then           ! Greyscale
         f=1.-f
         r=f
         g=f
         b=f
      else                            ! Colourscale
         r=0.
         g=0.
         b=0.
C Curves for RED
         if (f.ge.0.625.and.f.lt.0.875) r=1.
         if (f.ge.0.5.and.f.lt.0.625) then
            r=sqrt((f-0.5)/0.125)
         end if
         if (f.lt.0.25) then
            r=(0.25-f)/0.25
            r=r**(1.-sqrt(r))
         end if
         if (f.ge.0.875) then
            r=sqrt((1.-f)/0.125)
         end if
C Curves for GREEN
         if (f.ge.0.625.and.f.lt.0.775) then
            g=sqrt((0.775-f)/0.15)
         end if
         if (f.ge.0.375.and.f.lt.0.625) g=1.
         if (f.ge.0.25.and.f.lt.0.375) then
            g=sqrt((f-0.25)/0.125)
         end if
         if (f.lt.0.125) then
            g=sqrt((0.125-f)/0.125)
         end if
C Curves for BLUE
         if (f.ge.0.375.and.f.lt.0.5) then
            b=sqrt((0.5-f)/0.125)
         end if
         if (f.lt.0.375) b=1.
         if (f.ge.0.775) then
            b=sqrt((f-0.775)/0.225)
            b=b*r
         end if
C Rotate Colours if necessary
         if (jr.gt.5.or.jr.lt.0) jr=0
         if (jr.eq.1) then
            g1=b
            b=g
            g=g1
         else if (jr.eq.2) then
            r1=g
            g=b
            b=r
            r=r1
         else if (jr.eq.3) then
            r1=g
            g=r
            r=r1
         else if (jr.eq.4) then
            b1=g
            g=r
            r=b
            b=b1
         else if (jr.eq.5) then
            r1=b
            b=r
            r=r1
         end if
      end if
C Adjust Contrast
      if (xc.gt.1.) xc=1.
      if (xc.lt.0.) xc=1.
      r=(r-0.5)*xc+0.5
      g=(g-0.5)*xc+0.5
      b=(b-0.5)*xc+0.5
C Adjust brightness
      if (xb.gt.1.) xb=0.
      if (xb.lt.-1.) xb=0.
      if (xb.gt.0.) then
         r=1.-(1.-r)*(1.-xb)
         g=1.-(1.-g)*(1.-xb)
         b=1.-(1.-b)*(1.-xb)
      else if (xb.lt.0.) then
         r=r*(1.+xb)
         g=g*(1.+xb)
         b=b*(1.+xb)
      end if
C Apply gradient (colour only)
      if (jl.ne.0.and.idp.ne.0) then
      if (jl.eq.1) then
         r=r+(1.-r)*(1.-yf)
         g=g+(1.-g)*(1.-yf)
         b=b+(1.-b)*(1.-yf)
      else if (jl.ge.2) then
         r=r+(1.-r)*yf
         g=g+(1.-g)*yf
         b=b+(1.-b)*yf
      else if (jl.eq.-1) then
         r=r*yf
         g=g*yf
         b=b*yf
      else
         r=r*(1.-yf)
         g=g*(1.-yf)
         b=b*(1.-yf)
      end if
      end if
C Adjust Gamma
      if (xg.ne.1.) then
      r=r**xg
      g=g**xg
      b=b**xg
      end if
C Set Colour
      call gscr(1,nnn,r,g,b)
      return
      end      
      

      SUBROUTINE FILL (XWRK,YWRK,NWRK,IAREA,IGRP,NGRPS)
C
      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
C If the area is defined by 3 or more points, fill it
C
         CALL GSFACI(IAREA3+2)
         CALL GFA(NWRK,XWRK,YWRK)
      ENDIF
C
C Otherwise, do nothing
C
      RETURN
      END
