      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=99000,liwk=99000,lmap=1900000)
      parameter(nwrk=99000,isiz=5)
      real xwrk(nwrk),ywrk(nwrk),rwrk(lrwk),iwrk(liwk)
      integer map(lmap),iarea(isiz),igrp(isiz)
      real zlev(0:maxcon)
      real plim1(2),plim2(2),plim3(2),plim4(2)
      character lbls(0:maxcon)*4,label*80
      character concchar*80,levchar*5,agechar*20,specchar*10
      character asloragl*3,unitlabel*20,unitdepolabel*7
      character alspec*10,aiage*6

      parameter(nxmax=420,nymax=180,nzmax=29)
      parameter(maxpoint=800000,maxspec=4,maxageclass=15)
c     parameter(nxmax=720,nymax=180,nzmax=2)
c     parameter(maxpoint=5000,maxspec=268,maxageclass=1)
c     parameter(nxmax=120,nymax=80,nzmax=8)
c     parameter(maxpoint=50000,maxspec=50,maxageclass=12)


      real outheight(0: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 drygridtotal(nxmax,nymax,maxspec)
      real wetgridtotal(nxmax,nymax,maxspec)
      real gridtotal(nxmax,nymax,nzmax,maxspec)
      real z(nxmax,nymax),zw(nxmax,nymax),zd(nxmax,nymax)
      real zgrid(nxmax,nymax)
      real zt(nxmax,nymax),oro(nxmax,nymax),ihour,iday
      character dirname*110,datesname*110,maxvalue*10,minvalue*10
      character gridvalue*4
      character meanvalue*10,compoint(maxpoint)*10,species(maxspec)*10
      character filename*120,aday*8,atime*6

      real yra(5),xra(5)


      external cpdrpl,sfill_grau
      logical gridflag,dryflag,wetflag,totflag



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

      open(10,file='CONTROL')
      read(10,'(a110)') dirname
      read(10,'(a110)') datesname
      read(10,*) lconcormixratio
      read(10,*) decayconstant
      read(10,*) altitude,lev1,lev2  ! altitude level of plot, level range for column plots
      read(10,*) ltopo               ! asl or agl?
      read(10,'(a10)') alspec               ! kind of species to be plotted
c     read(10,*) nxstart,nxend       ! index of leftmost/rightmost column of data to plot
c     read(10,*) nystart,nyend       ! index of lowermost/uppermost row of data to plot
      read(10,*) rmin_col,rmax_col      ! degree leftmost/rightmost column of data to plot
      read(10,*) rmin_row,rmax_row      ! degree lowermost/uppermost row of data to plot
      read(10,*) lmark               ! mark release locations?
      read(10,*) lageplot            ! age classes in s to be plotted, 0 if all are to be plotted
      read(10,*) ldepoflag           ! plot deposition?
      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,*) scaledepo           ! scaling factor for deposition
      read(10,*) overrulemax         ! maximum value for conc./mixing ratios
      read(10,*) overruledmax        ! maximum value for deposition
      read(10,'(a20)') unitlabel      ! label for concentration unit
      read(10,'(a7)') unitdepolabel  ! label for deposition unit
      read(10,*) numcon
      read(10,*) ifact               ! number of values reported on labelbar
      read(10,*) xllon,yllat
      read(10,*) xulon,yulat
      read(10,*) iswitch
      close(10)
      if (numcon.gt.maxcon) numcon=maxcon
      ifact=numcon/(ifact-1)


      if (ltopo.eq.1) then  ! altitude in m asl
        asloragl='asl'
      else
        asloragl='agl'
      endif

      if ((altitude.lt.0.).and.(lconcormixratio.eq.2))
     +stop 'Cannot use mixing ratio output to calculate total columns,
     +use file grid_conc for that'
      len=index(dirname,' ')-1



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

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

      call pcseti('FN',22)
      call setusv('LW',1000)

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
c          write(*,*) ibdate,ibtime
      read(10) loutstep,loutaver,loutsample
      read(10) outlon0,outlat0,numxgrid,numygrid,dxout,dyout
c           write(*,*) outlon0,outlat0,numxgrid,numygrid,dxout,dyout
      read(10) numzgrid,(outheight(i),i=1,numzgrid)
      read(10) jjjjmmdd,ihmmss
c          write(*,*) jjjjmmdd,ihmss
      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)

      outheight(0)=0.
      close(10)
c         write(*,*) (outheight(i),i=1,numzgrid)

      if (loutstep.lt.0) nspec=numpoint
c          write(*,*) nspec

      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 if (lconcormixratio.eq.3) then
          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 if (lconcormixratio.eq.3) then
          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 If negative altitude is selected, columns from lev1 to lev2 are calculated
C If total column shall be plotted, either lev1 or lev2 can be set to zero
****************************************************************************

      if (lev2.gt.numzgrid) lev2=numzgrid
      if ((altitude.lt.0.).and.((lev1.le.0).or.(lev2.le.0))) then
        lev1=1
        lev2=numzgrid
      endif


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

      do 150 ix=1,numxgrid
        do 150 jy=1,numygrid
          if (ltopo.eq.1) then
            heightnn (ix,jy,0) = oro(ix,jy)
          else
            heightnn (ix,jy,0) = 0.
          endif
          do 150 i=1,numzgrid
            if (ltopo.eq.1) then
              heightnn (ix,jy,i) = outheight(i) + oro(ix,jy)
            else
              heightnn (ix,jy,i) = outheight(i)
            endif
150         continue         


C Set up sub-domain to be plotted
*********************************

      nxstart = int((rmin_col - outlon0)/dxout + 1)
      nxend = int((rmax_col - outlon0)/dxout + 1)
      nystart = int((rmin_row - outlat0)/dyout + 1)
      nyend = int((rmax_row - outlat0)/dyout + 1)

c     write(*,*) nxstart, nxend, numxgrid
c     write(*,*) nystart, nyend, numygrid

      if (nxstart.lt.1) nxstart = 1
      if (nystart.lt.1) nystart = 1
      if (nxend.le.1) nxend = 2 
      if (nyend.le.1) nyend = 2

      if (nxstart.gt.numxgrid) nxstart = numxgrid-1
      if (nystart.gt.numygrid) nystart = numygrid-1
      if (nxend.gt.numxgrid) nxend = numxgrid
      if (nyend.gt.numygrid) nyend = numygrid

c     write(*,*) nxstart, nxend, numxgrid
c     write(*,*) nystart, nyend, numygrid

      if (nxstart.gt.nxend) stop 'domain settings incorrect'
      if (nystart.gt.nyend) stop 'domain settings incorrect'

!!! Fix: NCAR doesn't plot global fields for unknown reason

      if (float(nxend-nxstart+1)*dxout.gt.359.99) nxend=nxend-1
      nx=nxend-nxstart+1       ! number of columns to plot
      ny=nyend-nystart+1       ! number of rows to plot

C Set map outline
*****************


      plim1(1)=outlat0+dyout/2.+float(nystart-1)*dyout
      plim2(1)=outlon0+dxout/2.+float(nxstart-1)*dxout
      plim3(1)=float(nyend-1)*dyout+dyout/2.+outlat0
      plim4(1)=float(nxend-1)*dxout+dxout/2.+outlon0
      rlatmid=0.
      rlonmid=0.

      rgridist=1.
      if ((plim4(1)-plim2(1)).ge.180.) rgridist=20.
      if ((plim4(1)-plim2(1)).ge.90..and.(plim4(1)-plim2(1)).lt.180.) 
     +                                  rgridist=10.
      if ((plim4(1)-plim2(1)).ge.20..and.(plim4(1)-plim2(1)).lt.90.) 
     +                                  rgridist=5.

      write(gridvalue,'(f4.1)') rgridist

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

C Loop about all times, given in file dates
*******************************************

      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,inday,intime

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
            wetgridtotal(ix,jy,k)=0.
            drygridtotal(ix,jy,k)=0.
            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
                read(10) (grid(ix,jy,kz,k,nage),jy=1,numygrid)
39              continue
          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
              drygrid(ix,jy,k,nage)=drygrid(ix,jy,k,nage)*scaledepo
              wetgrid(ix,jy,k,nage)=wetgrid(ix,jy,k,nage)*scaledepo
              drygridtotal(ix,jy,k)=drygridtotal(ix,jy,k)+
     +        drygrid(ix,jy,k,nage)
              wetgridtotal(ix,jy,k)=wetgridtotal(ix,jy,k)+
     +        wetgrid(ix,jy,k,nage)
              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)
41              continue

 49      continue

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

      do 200 k=1,nspec

      specchar=alspec
      write(levchar,'(i5)') nint(altitude)

      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

        dryflag=.false.
        wetflag=.false.
        totflag=.false.
        gridflag=.false.


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

        zmax=0.
        zmin=1.e30
        zmean=0.
        zwmax=0.
        zwmin=1.e30
        zwmean=0.
        zdmax=0.
        zdmin=1.e30
        zdmean=0.
        ztmax=0.
        ztmin=1.e30
        ztmean=0.
        numbmean=0
	do 201 ix=1,numxgrid
	  do 201 jy=1,numygrid
              zw(ix,jy)=wetgridtotal(ix,jy,k)
              zd(ix,jy)=drygridtotal(ix,jy,k)
              zt(ix,jy)=wetgridtotal(ix,jy,k)+drygridtotal(ix,jy,k)
              if (altitude.gt.0.) then
                do 220 kz=1,numzgrid
                  if (heightnn(ix,jy,kz).ge.altitude) goto 221
220               continue
221             kz=min(kz,numzgrid)
                z(ix,jy)=gridtotal(ix,jy,kz,k)
                if (heightnn(ix,jy,0).ge.altitude) z(ix,jy)=1.e30  ! below topo
              else       ! calculate total column (mass/area)
                z(ix,jy)=0.
                do 203 kz=lev1,lev2
203               z(ix,jy)=z(ix,jy)+gridtotal(ix,jy,kz,k)*
     +            (outheight(kz)-outheight(kz-1))            ! ng/m2
              endif
           if ((z(ix,jy).lt.9.e29).and.(z(ix,jy).gt.zmax)) zmax=z(ix,jy)
            if (z(ix,jy).lt.zmin) zmin=z(ix,jy)
            if (z(ix,jy).lt.9.e29) then
              zmean=zmean+z(ix,jy)
              numbmean=numbmean+1
            endif
            if (zw(ix,jy).gt.zwmax) zwmax=zw(ix,jy)
            if (zw(ix,jy).lt.zwmin) zwmin=zw(ix,jy)
            zwmean=zwmean+zw(ix,jy)
            if (zd(ix,jy).gt.zdmax) zdmax=zd(ix,jy)
            if (zd(ix,jy).lt.zdmin) zdmin=zw(ix,jy)
            zdmean=zdmean+zd(ix,jy)
            if (zt(ix,jy).gt.ztmax) ztmax=zt(ix,jy)
            if (zt(ix,jy).lt.ztmin) ztmin=zw(ix,jy)
            ztmean=ztmean+zt(ix,jy)
201         continue
        zmean=zmean/float(numbmean)
        zwmean=zwmean/float(numxgrid)/float(numygrid)
        zdmean=zdmean/float(numxgrid)/float(numygrid)
        ztmean=ztmean/float(numxgrid)/float(numygrid)


        if (zdmax.gt.0.) dryflag=.true.
        if (zwmax.gt.0.) wetflag=.true.
        if (ztmax.gt.0.) totflag=.true.
        if (zmax.gt.0.) gridflag=.true.


C Draw map
**********


      if (gridflag) then
      call mapset('CO',plim1,plim2,plim3,plim4)   ! select area of plot
      call mpseti('LA',0)      ! no labels are drawn at GM
      call maproj('CE',rlatmid,rlonmid,0.)


      call mappos(0.02,0.98,0.15,0.84) ! set map position
      call mapstc('OU','CO')           ! choose shorelines
      call mapint    ! initialize map

      call arinam(map,lmap)         ! initialize areas
      call mapbla(map)              ! add areas to map
      call mpsetr('GR',rgridist)  ! set resolution of lat/long grid
      call mpsetr('GD',1.)
      call mapgrd                   ! plot lat/long grid
      call maplot    ! plot the map
      call maplbl    ! draw a perimeter

      include 'include.contour'

      do 44 i=nystart,nyend
        do 44 j=nxstart,nxend
          zgrid(j-nxstart+1,i-nystart+1)=z(j,i)
44        continue



C Settings of conpack
*********************

      call cpseti('SET',0)
      call cpsetc('ILT - Information Label Text String',' ')
      call cpsetc('HLT - High/Low label Text Strings', ' '' ')
      call cpseti('RWC - Real Workspace for Contours',500)
      call cpseti('CLS - Contour Level Selection',0)
      call cpseti('NCL - Number of Contour Levels',numcon+1)
      call cpsetr('SPV',1.e30)

      call arinam(map,lmap)
      call cprect(zgrid,nxmax,nx,ny,rwrk,lrwk,iwrk,liwk)
      call cppkcl(zgrid,rwrk,iwrk)

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


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

      if (overrulemax.gt.0.) then
        zmax=overrulemax
        zmin=0.
      endif

      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)=zlev(0)+float(i)/float(numcon)*(zmax-zmin)
          call cpsetr('CLV - CONTOUR LEVEL VALUES',zlev(i-1))
27        call cpseti('CLU - Contour Level Use Flag',2)
      endif
     
      call cpback(zgrid,rwrk,iwrk)      ! Draw perimeter
      call cpclam(zgrid,rwrk,iwrk,map)  ! Add contours to area map
      call cplbam(zgrid,rwrk,iwrk,map)  ! Add label boxes to area map

C Define the area shading

      call arscam(map,xwrk,ywrk,nwrk,iarea,igrp,isiz,sfill_grau)
      call cpcldm(zgrid,rwrk,iwrk,map,cpdrpl) ! draw contours
      call cplbdr(zgrid,rwrk,iwrk)            ! draw labels

      call mapgrd                   ! plot lat/long grid
      call maplot    ! plot the map
      call maplbl    ! draw a perimeter

C Mark release points
*********************

      if (lmark.eq.1) then
      do 37 i=1,numpoint
        call maptra(ypoint(i),xpoint(i),xval,yval)
        if (xval.ne.1.e12) call points(xval,yval,1,-3,0)
37      continue
      endif

      do 34 ix=1,numxgrid
        do 34 jy=1,numygrid
          if (z(ix,jy).gt.9.e29) then     ! mark regions below topography
            call maptra((float(jy)-0.5)*dyout+outlat0,
     +      (float(ix)-0.5)*dxout+outlon0,xval,yval)
            if (xval.ne.1.e12) call points(xval,yval,1,-2,0)
          endif
34        continue


C Draw line of cross-section
****************************
     
      if (xllon.lt.plim2(1)) xllon=plim2(1)
      if (xulon.gt.plim4(1)) xulon=plim4(1)
      if (yllat.lt.plim1(1)) yllat=plim1(1)
      if (yulat.gt.plim3(1)) yulat=plim3(1)

      call setusv('LW',3000)
      call mapfst (yllat,xllon)
      call mapvec (yulat,xulon)
      call setusv('LW',1000)


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.115
      ylic2=0.140
      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.10,lbls(i),.015,0.,0.)
      else if (i.eq.numcon) then
        call plchhq(xlp-0.010,0.10,lbls(i),.015,0.,0.)
      else if (mod(i,ifact).eq.0) then
        call plchhq(xlp,0.10,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,.070,unitlabel,.015,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)-(iday*24) 
      iihour=nint(ihour)
      ibtime_3=ibtime-30000
      if (inday.lt.ibdate) then
      write(label,901) inday,intime
      else
      write(label,900) ibdate,ibtime_3,inday,intime
      endif
900   format('Latest analysis time ',i8,'.',i6,3x,
     +'Actual time ',i8,'.',i6)
901   format('Actual time ',i8,'.',i6)
      call set(0.,1.,0.,1.,0.,1.,0.,1.,1)
      if (altitude.gt.0.) then
      if (lconcormixratio.eq.2) then
        concchar='Mixing ratio of species '//specchar//' at '
     +  //levchar//' m '//asloragl//' for age class '//agechar
      else
        concchar='Concentration of species '//specchar//' at '
     +  //levchar//' m '//asloragl//' for age class '//agechar
      endif
      else
        concchar='Total column of species '//specchar//
     +  ' for age class '//agechar
      endif

      call plchhq(.5,.95,label,.016,0.,0.)
c     call plchhq(.5,.92,'Mean value '//meanvalue,.016,0.,0.)
c     call plchhq(.5,.892,'Maximum value '//maxvalue,.016,0.,0.)
c     call plchhq(.5,.864,'Minimum value '//minvalue,.016,0.,0.)
      call plchhq(.5,.92,'Mean value '//meanvalue//'  Maximum value '
     +//maxvalue//'  Minimum value '//minvalue,.016,0.,0.)
      call plchhq(.5,.892,'Distance of grid lines '//gridvalue//' deg'
     +,.016,0.,0.)
      call pcseti('FN',22)
      call plchhq(.5,.98,concchar,.017,0.,0.)



      call frame
      endif


C Plot depositions as well?
***************************

      if (ldepoflag.ne.1) goto 200


C DRAW WET DEPOSITION
*********************

      if (wetflag) then
      call mapset('CO',plim1,plim2,plim3,plim4)   ! select area of plot
      call mpseti('LA',0)      ! no labels are drawn at GM
      call maproj('CE',rlatmid,rlonmid,0.)


      call mappos(0.02,0.98,0.15,0.84) ! set map position
      call mapstc('OU','CO')           ! choose shorelines
      call mapint    ! initialize map

      call arinam(map,lmap)         ! initialize areas
      call mapbla(map)              ! add areas to map
      call mpsetr('GR',rgridist)  ! set resolution of lat/long grid
      call mpsetr('GD',1.)
      call mapgrd                   ! plot lat/long grid
      call maplot    ! plot the map
      call maplbl    ! draw a perimeter

      include 'include.contour'

      do 144 i=nystart,nyend
        do 144 j=nxstart,nxend
          zgrid(j-nxstart+1,i-nystart+1)=zw(j,i)
144       continue


C Settings of conpack
*********************

      call cpseti('SET',0)
      call cpsetc('ILT - Information Label Text String',' ')
      call cpseti('RWC - Real Workspace for Contours',500)
      call cpseti('CLS - Contour Level Selection',0)
      call cpseti('NCL - Number of Contour Levels',numcon+1)

      call arinam(map,lmap)
      call cprect(zgrid,nxmax,nx,ny,rwrk,lrwk,iwrk,liwk)
      call cppkcl(zgrid,rwrk,iwrk)

      write(maxvalue,'(e10.3)') zwmax
      write(minvalue,'(e10.3)') zwmin
      write(meanvalue,'(e10.3)') zwmean


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

      if (overruledmax.gt.0.) zwmax=overruledmax

      if (loglin.eq.1) then
        zwmax=float(int(2.*log10(zwmax)))/2.
        zlev(0)=10**(zwmax-float(numcon)/2.)
        do 128 i=1,numcon
          call cpseti('PAI - Parameter Array Index',i)
          zlev(i)=10**(zwmax-float(numcon-i)/2.)
          call cpsetr('CLV - CONTOUR LEVEL VALUES',zlev(i))
128       call cpseti('CLU - Contour Level Use Flag',2)
      else
        zlev(0)=zwmin
        do 127 i=1,numcon
          call cpseti('PAI - Parameter Array Index',i)
          zlev(i)=zlev(0)+float(i)/float(numcon)*zwmax
          call cpsetr('CLV - CONTOUR LEVEL VALUES',zlev(i-1))
127       call cpseti('CLU - Contour Level Use Flag',2)
      endif
     
      call cpback(zgrid,rwrk,iwrk)      ! Draw perimeter
      call cpclam(zgrid,rwrk,iwrk,map)  ! Add contours to area map
      call cplbam(zgrid,rwrk,iwrk,map)  ! Add label boxes to area map

C Define the area shading

      call arscam(map,xwrk,ywrk,nwrk,iarea,igrp,isiz,sfill_grau)
      call cpcldm(zgrid,rwrk,iwrk,map,cpdrpl) ! draw contours
      call cplbdr(zgrid,rwrk,iwrk)            ! draw labels

      call mapgrd                   ! plot lat/long grid
      call maplot    ! plot the map
      call maplbl    ! draw a perimeter


C Mark release points
*********************

      if (lmark.eq.1) then
      do 137 i=1,numpoint
        call maptra(ypoint(i),xpoint(i),xval,yval)
        if (xval.ne.1.e12) call points(xval,yval,1,-3,0)
137     continue
      endif


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

      if (loglin.eq.1) then
        do 126 i=0,numcon
126       write(lbls(i),'(i4)') nint(10**(zwmax-float(numcon-1-i)/2.))
      else
        do 124 i=0,numcon
124       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.115
      ylic2=0.140
      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.10,lbls(i),.015,0.,0.)
      else if (i.eq.numcon) then
        call plchhq(xlp,0.10,lbls(i),.015,0.,0.)
      else if (mod(i,ifact).eq.0) then
        call plchhq(xlp,0.10,lbls(i),.015,0.,0.)
      endif
      end do
      call setusv('LW',2000)
      do i=1,numcon-1
      xlp=xlic1+(xlic2-xlic1)*float(i)/float(numcon)
      call frstd(xlp,ylic1)
      call vectd(xlp,ylic2)
      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,.070,unitdepolabel,.015,0.,0.)

      call setusv('LW',1000)



C Draw the release date and time and current time
*************************************************

      write(label,900) ibdate,ibtime,itime/3600
      call set(0.,1.,0.,1.,0.,1.,0.,1.,1)
      concchar='Wet deposition of species '//specchar//
     +' for age class '//agechar

      call plchhq(.5,.95,label,.016,0.,0.)
      call plchhq(.5,.92,'Mean value '//meanvalue,.016,0.,0.)
      call plchhq(.5,.892,'Maximum value '//maxvalue,.016,0.,0.)
      call plchhq(.5,.864,'Minimum value '//minvalue,.016,0.,0.)
      call pcseti('FN',22)
      call plchhq(.5,.98,concchar,.017,0.,0.)


      call frame
      endif


C DRY DEPOSITION


      if (dryflag) then
      call mapset('CO',plim1,plim2,plim3,plim4)   ! select area of plot
      call mpseti('LA',0)      ! no labels are drawn at GM
      call maproj('CE',rlatmid,rlonmid,0.)


      call mappos(0.02,0.98,0.15,0.84) ! set map position
      call mapstc('OU','CO')           ! choose shorelines
      call mapint    ! initialize map

      call arinam(map,lmap)         ! initialize areas
      call mapbla(map)              ! add areas to map
      call mpsetr('GR',rgridist)  ! set resolution of lat/long grid
      call mpsetr('GD',1.)
      call mapgrd                   ! plot lat/long grid
      call maplot    ! plot the map
      call maplbl    ! draw a perimeter

      include 'include.contour'

      do 244 i=nystart,nyend
        do 244 j=nxstart,nxend
          zgrid(j-nxstart+1,i-nystart+1)=zd(j,i)
244       continue


C Settings of conpack
*********************

      call cpseti('SET',0)
      call cpsetc('ILT - Information Label Text String',' ')
      call cpseti('RWC - Real Workspace for Contours',500)
      call cpseti('CLS - Contour Level Selection',0)
      call cpseti('NCL - Number of Contour Levels',numcon+1)

      call arinam(map,lmap)
      call cprect(zgrid,nxmax,nx,ny,rwrk,lrwk,iwrk,liwk)
      call cppkcl(zgrid,rwrk,iwrk)

     
      write(maxvalue,'(e10.3)') zdmax
      write(minvalue,'(e10.3)') zdmin
      write(meanvalue,'(e10.3)') zdmean


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

      if (overruledmax.gt.0.) zdmax=overruledmax

      if (loglin.eq.1) then
        zdmax=float(int(2.*log10(zdmax)))/2.
        zlev(0)=10**(zdmax-float(numcon)/2.)
        do 228 i=1,numcon
          call cpseti('PAI - Parameter Array Index',i)
          zlev(i)=10**(zdmax-float(numcon-i)/2.)
          call cpsetr('CLV - CONTOUR LEVEL VALUES',zlev(i))
228       call cpseti('CLU - Contour Level Use Flag',2)
      else
        zlev(0)=zdmin
        do 227 i=1,numcon
          call cpseti('PAI - Parameter Array Index',i)
          zlev(i)=zlev(0)+float(i)/float(numcon)*zdmax
          call cpsetr('CLV - CONTOUR LEVEL VALUES',zlev(i-1))
227       call cpseti('CLU - Contour Level Use Flag',2)
      endif
     
      call cpback(zgrid,rwrk,iwrk)      ! Draw perimeter
      call cpclam(zgrid,rwrk,iwrk,map)  ! Add contours to area map
      call cplbam(zgrid,rwrk,iwrk,map)  ! Add label boxes to area map

C Define the area shading

      call arscam(map,xwrk,ywrk,nwrk,iarea,igrp,isiz,sfill_grau)
      call cpcldm(zgrid,rwrk,iwrk,map,cpdrpl) ! draw contours
      call cplbdr(zgrid,rwrk,iwrk)            ! draw labels

      call mapgrd                   ! plot lat/long grid
      call maplot    ! plot the map
      call maplbl    ! draw a perimeter


C Mark release points
*********************

      if (lmark.eq.1) then
      do 237 i=1,numpoint
        call maptra(ypoint(i),xpoint(i),xval,yval)
237     if (xval.ne.1.e12) call points(xval,yval,1,-3,0)
      endif


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

      if (loglin.eq.1) then
        do 226 i=0,numcon
226       write(lbls(i),'(i4)') nint(10**(zdmax-float(numcon-1-i)/2.))
      else
        do 224 i=0,numcon
224       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.115
      ylic2=0.140
      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.10,lbls(i),.015,0.,0.)
      else if (i.eq.numcon) then
        call plchhq(xlp,0.10,lbls(i),.015,0.,0.)
      else if (mod(i,ifact).eq.0) then
        call plchhq(xlp,0.10,lbls(i),.015,0.,0.)
      endif
      end do
      call setusv('LW',2000)
      do i=1,numcon-1
      xlp=xlic1+(xlic2-xlic1)*float(i)/float(numcon)
      call frstd(xlp,ylic1)
      call vectd(xlp,ylic2)
      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,.070,unitdepolabel,.016,0.,0.)

      call setusv('LW',1000)


C Draw the release date and time and current time
*************************************************

      write(label,900) ibdate,ibtime,itime/3600
      call set(0.,1.,0.,1.,0.,1.,0.,1.,1)
      concchar='Dry deposition of species '//specchar//
     +' for age class '//agechar

      call plchhq(.5,.95,label,.016,0.,0.)
      call plchhq(.5,.92,'Mean value '//meanvalue,.016,0.,0.)
      call plchhq(.5,.892,'Maximum value '//maxvalue,.016,0.,0.)
      call plchhq(.5,.864,'Minimum value '//minvalue,.016,0.,0.)
      call pcseti('FN',22)
      call plchhq(.5,.98,concchar,.017,0.,0.)


      call frame
      endif



C TOTAL DEPOSITION


      if (totflag) then
      call mapset('CO',plim1,plim2,plim3,plim4)   ! select area of plot
      call mpseti('LA',0)      ! no labels are drawn at GM
      call maproj('CE',rlatmid,rlonmid,0.)


      call mappos(0.02,0.98,0.15,0.84) ! set map position
      call mapstc('OU','CO')           ! choose shorelines
      call mapint    ! initialize map

      call arinam(map,lmap)         ! initialize areas
      call mapbla(map)              ! add areas to map
      call mpsetr('GR',rgridist)  ! set resolution of lat/long grid
      call mpsetr('GD',1.)
      call mapgrd                   ! plot lat/long grid
      call maplot    ! plot the map
      call maplbl    ! draw a perimeter

      include 'include.contour'

      do 344 i=nystart,nyend
        do 344 j=nxstart,nxend
          zgrid(j-nxstart+1,i-nystart+1)=zt(j,i)
344       continue


C Settings of conpack
*********************

      call cpseti('SET',0)
      call cpsetc('ILT - Information Label Text String',' ')
      call cpseti('RWC - Real Workspace for Contours',500)
      call cpseti('CLS - Contour Level Selection',0)
      call cpseti('NCL - Number of Contour Levels',numcon+1)

      call arinam(map,lmap)
      call cprect(zgrid,nxmax,nx,ny,rwrk,lrwk,iwrk,liwk)
      call cppkcl(zgrid,rwrk,iwrk)

      write(maxvalue,'(e10.3)') ztmax
      write(minvalue,'(e10.3)') ztmin
      write(meanvalue,'(e10.3)') ztmean


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

      if (overruledmax.gt.0.) ztmax=overruledmax

      if (loglin.eq.1) then
        ztmax=float(int(2.*log10(ztmax)))/2.
        zlev(0)=10**(ztmax-float(numcon)/2.)
        do 328 i=1,numcon
          call cpseti('PAI - Parameter Array Index',i)
          zlev(i)=10**(ztmax-float(numcon-i)/2.)
          call cpsetr('CLV - CONTOUR LEVEL VALUES',zlev(i))
328       call cpseti('CLU - Contour Level Use Flag',2)
      else
        zlev(0)=ztmin
        do 327 i=1,numcon
          call cpseti('PAI - Parameter Array Index',i)
          zlev(i)=zlev(0)+float(i)/float(numcon)*ztmax
          call cpsetr('CLV - CONTOUR LEVEL VALUES',zlev(i-1))
327       call cpseti('CLU - Contour Level Use Flag',2)
      endif
     
      call cpback(zgrid,rwrk,iwrk)      ! Draw perimeter
      call cpclam(zgrid,rwrk,iwrk,map)  ! Add contours to area map
      call cplbam(zgrid,rwrk,iwrk,map)  ! Add label boxes to area map

C Define the area shading

      call arscam(map,xwrk,ywrk,nwrk,iarea,igrp,isiz,sfill_grau)
      call cpcldm(zgrid,rwrk,iwrk,map,cpdrpl) ! draw contours
      call cplbdr(zgrid,rwrk,iwrk)            ! draw labels

      call mapgrd                   ! plot lat/long grid
      call maplot    ! plot the map
      call maplbl    ! draw a perimeter


C Mark release points
*********************

      if (lmark.eq.1) then
      do 337 i=1,numpoint
        call maptra(ypoint(i),xpoint(i),xval,yval)
337     if (xval.ne.1.e12) call points(xval,yval,1,-3,0)
      endif


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

      if (loglin.eq.1) then
        do 326 i=0,numcon
326       write(lbls(i),'(i4)') nint(10**(ztmax-float(numcon-1-i)/2.))
      else
        do 324 i=0,numcon
324       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.115
      ylic2=0.140
      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.10,lbls(i),.015,0.,0.)
      else if (i.eq.numcon) then
        call plchhq(xlp,0.10,lbls(i),.015,0.,0.)
      else if (mod(i,ifact).eq.0) then
        call plchhq(xlp,0.10,lbls(i),.015,0.,0.)
      endif
      end do
      call setusv('LW',2000)
      do i=1,numcon-1
      xlp=xlic1+(xlic2-xlic1)*float(i)/float(numcon)
      call frstd(xlp,ylic1)
      call vectd(xlp,ylic2)
      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,.070,unitdepolabel,.016,0.,0.)

      call setusv('LW',1000)


C Draw the release date and time and current time
*************************************************

      write(label,900) ibdate,ibtime,itime/3600
      call set(0.,1.,0.,1.,0.,1.,0.,1.,1)
      concchar='Tot deposition of species '//specchar//
     +' for age class '//agechar

      call plchhq(.5,.95,label,.016,0.,0.)
      call plchhq(.5,.92,'Mean value '//meanvalue,.016,0.,0.)
      call plchhq(.5,.892,'Maximum value '//maxvalue,.016,0.,0.)
      call plchhq(.5,.864,'Minimum value '//minvalue,.016,0.,0.)
      call pcseti('FN',22)
      call plchhq(.5,.98,concchar,.017,0.,0.)

      call frame
      endif


200   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      
      

