      program flexplot_ht
C
C $Id: flexplot_h_t.f,v 1.16 2001/08/28 12:28:26 pauli Exp $
C
c change 10/10/2005 Caroline Forster
c plots time-height profiles of tracer along flighttracks
c 

C options (read from CONTROL - file)
      character*120 dirname,datesname,filename,datafilename
      character aday*8,atime*6

      integer   lspec                ! species to be plotted
      character*10 alspec              ! species to be plotted
      integer   tstart, tend       ! start and end (time)
      real      pktfi, pktlambda   ! coordinates of desired point   lon/lat
      integer   do_cont            ! draw contours (0 off, 1 on)
      integer   do_fill            ! fill with colors (0 off, 1 on)
      integer   do_color           ! 
      integer   do_linlab          ! draw line labels
      integer   do_gridpts         ! mark gridpoints
      integer   do_lblbar          ! draw labelbar (0 off, 1 on)
      character unitlabel*10        ! label for concentration unit
      real      minheight          ! min. height to plot (m ASL) -1. = orography
      real      maxheight          ! max. height to plot (m ASL) 0. = max modelh

      integer   do_surface         ! plot surfaceline (only if startlev is MSL)
      integer   fixcont            ! fixed conour-intervals specified


      integer pktx,pkty            ! grid-coorinates of desired point

      integer maxncon
      parameter (maxncon=220)

      integer ierrf, iworkmeta
      integer lunitmeta,idmeta
      integer lrwk,liwk,lmap,nwrk,nogrps,lzrg
      parameter(ierrf=6,iworkmeta=1)
      parameter(lunitmeta=82, idmeta=1)
      parameter(lrwk=299000,liwk=299000)
      parameter(lmap=2900000,nwrk=299000)
      parameter(nogrps=299000,lzrg=2*49900000)
      real xwrk(nwrk),ywrk(nwrk),rwrk(lrwk),iwrk(liwk)
      integer map(lmap)
      integer iarea(nogrps),igrp(nogrps)

      integer lfin(maxncon)
      real zlev(0:maxncon)
      real plim1,plim2,plim3,plim4           ! limits
      character lbls(0:maxncon)*6            ! labels for labelbar
      character latstr*5, lonstr*6
      character label*39 

      real ihour,iday,i2hour,i2day

C Variables for reading model results
      integer nxmax,nymax,nzmax,ntmax,maxpoint,maxspec,maxageclass
      parameter(nxmax=360,nymax=180,nzmax=30,ntmax=5000)
      parameter(maxpoint=200000,maxspec=3,maxageclass=6)
      integer ireleasestart(maxpoint),ireleaseend(maxpoint)
      integer npart(maxpoint),kind(maxpoint),lage(0:maxageclass)
      real outheight(0:nzmax)
      real oro(nxmax,nymax)         ! topography
      real xpoint(maxpoint),ypoint(maxpoint)
      real zpoint1(maxpoint), zpoint2(maxpoint)
      real 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)
      real zw(nxmax,nymax),zd(nxmax,nymax)
      character compoint(maxpoint)*45,species(maxspec)*10
      real heightnn(nxmax,nymax,0:nzmax)
      character dummy2*32

c for flightleg  information

      parameter(maxdat=100000)
      parameter(maxdates=10000)
      real seconds(maxdat)
      real xlon(maxdat),ylat(maxdat),alt(maxdat),press(maxdat)
      real topo(maxdat)
      real dummy(50)
      double precision juldate,julstart,jul,juldata(maxdat)
      double precision julin(maxdates)
      integer inday(maxdates), intime(maxdates)
      integer jjjjmmdd(maxdat),ihmmss(maxdat)

C plotting grid
      real x(ntmax), y(0:nzmax), z(ntmax,0:nzmax)  ! irregularly spaced gridded data
      real y_meter(0:nzmax)
      real zreg(lzrg)                          ! calculated regular grid

      integer   i, hlab ,kkk
      integer   numzgrid, maxzgrid, minzgrid
      character maxvalue*10,minvalue*10,meanvalue*10

      real yra(5),xra(5)

      external cpdrpl, fill


C Read CONTROL file

      open(10,file='CONTROL')
      read(10,'(a120)') dirname
      read(10,'(a120)') datesname
      read(10,'(a120)') datafilename
      read(10,*) iswitch
      read(10,*) lconcormixratio
      read(10,*) decayconstant
      read(10,'(a10)') alspec
      read(10,*) do_cont
      read(10,*) do_fill
      read(10,*) do_color
      read(10,*) do_linlab
      read(10,*) do_gridpts
      read(10,*) do_lblbar
      read(10,*) scalefact
      read(10,*) zoverrulemax
      read(10,'(a10)') unitlabel
      read(10,*) minheight
      read(10,*) maxheight
      read(10,*) do_surface
      read(10,*) numcon
      read(10,*) ifact
      read(10,*) fixcont
      if (fixcont.eq.1) then
         do 9 i=1,numcon-1
            read(10,*) zlev(i)
9        continue
      endif
      if (maxheight.le.0.) maxheight=30000.
      maxheight=maxheight/1000.0
      if (minheight.le.0.) minheight=0.
      minheight=minheight/1000.0
      if (do_fill.eq.0) do_lblbar=0
      if (startlev.eq.1) do_surface=0
      if (numcon.gt.maxncon)  numcon=maxncon
      ifact=numcon/(ifact-1)
      close(10)

      len=index(dirname,' ')-1
      if (iswitch.eq.0) then
        filename=dirname(1:len)//'header_nest'
      else
        filename=dirname(1:len)//'header'
      endif



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

      call gopks(ierrf,0)
      call gopwk(idmeta,lunitmeta,iworkmeta)
      call gacwk(idmeta)
      call gstxfp(-13,2)


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

c read flightleg information
****************************

      eps=0.0001

c     write(*,*) datafilename, datesname

      open(20,file=datafilename,status='old')
      i=1
100     read(20,*,end=99) jjjjmmdd(i),ihmmss(i),ylat(i),xlon(i),
     +  alt(i)
        if (i.gt.1.and.jjjjmmdd(i).eq.jjjjmmdd(i-1)
     +            .and.ihmmss(i).eq.ihmmss(i-1)
     +            .and.abs(ylat(i)-ylat(i-1)).le.eps
     +            .and.abs(xlon(i)-xlon(i-1)).le.eps
     +            .and.abs(alt(i)-alt(i-1)).le.eps) goto 100
        juldata(i)=juldate(jjjjmmdd(i),ihmmss(i))
        if (i.eq.1) julstart=juldate(jjjjmmdd(1),0)
        seconds(i)=(juldata(i)-julstart)*dble(86400.)
        i=i+1
        goto 100

99    numb=i-1

       close(20)


        tstart=sngl(juldate(jjjjmmdd(1),ihmmss(1))-
     +    juldate(jjjjmmdd(1),0))*86400.
        tend=sngl(juldate(jjjjmmdd(numb),ihmmss(numb))-
     +    juldate(jjjjmmdd(1),0))*86400.

C Read available dates from model simulation
********************************************

      open(20,file=datesname,form='formatted',status='old')
      i=1
55    read(20,'(i8,i6)',end=49) inday(i),intime(i)
        julin(i)=juldate(inday(i),intime(i))
        i=i+1
        goto  55
49    numdates=i-1
      close(20)


C Read header file of model results
***********************************

      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) jjjmmdd,jhmmss
c          write(*,*) jjjmmdd,jhmmss
      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)(1:8),dummy2
        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)

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

      if (loutstep.lt.0) nspec=numpoint

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

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

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

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


c set filenames
***************

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

c initialize
************

      zmax=0.
      zmin=1.e30
      zmean=0.
      izaehl=0

C Loop over all available measurement data
******************************************

      do num=1,numb

c find appropriate model output
*******************************

      do iout=2,numdates
       if ((julin(iout-1).lt.juldata(num)).and.
     +    (julin(iout).ge.juldata(num))) then
          if (iout.ne.ioutold) then
            ioutold=iout
            write(aday,'(i8.8)') inday(iout)
            write(atime,'(i6.6)') intime(iout)
            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

            read(10,end=99) itime
c           write(*,*) aday//atime


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

      do 38 k=1,nspec

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
            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
c30            read(10) (wetgrid(ix,jy,k,nage),yuncert,jy=1,numygrid)
 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
c31            read(10) (drygrid(ix,jy,k,nage),yuncert,jy=1,numygrid)
 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
c39              read(10) (grid(ix,jy,kz,k,nage),yuncert,jy=1,numygrid)
 39              read(10) (grid(ix,jy,kz,k,nage),jy=1,numygrid)
          endif

78        continue


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

          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)*scalefact
                gridtotal(ix,jy,kz,k)=gridtotal(ix,jy,kz,k)+
     +          grid(ix,jy,kz,k,nage)
41              continue


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

38      continue

        close(10)

        endif ! (iout.ne.ioutold)
        endif ! ((julin(iout-1).lt.juldata(num)).and.(julin(iout).ge.juldata(num)))

        end do  ! (iout)

C Determine which point to be plotted
*************************************

      pktx = int((xlon(num)-outlon0) / dxout)+1
      pkty = int((ylat(num)-outlat0) / dyout)+1
      if (numxgrid.lt.pktx.or.pktx.lt.0)  then
       write(*,*) 'ERROR : aircraft positions not in output grid ! '
       write(*,*) 'time-height section not plotted ! '
       stop
      endif
      if (numygrid.lt.pkty.or.pkty.lt.0) then
       write(*,*) 'ERROR : aircraft positions not in output grid ! '
       write(*,*) 'time-height section not plotted ! '
       stop
      endif


C Determine altitude levels used to plot
****************************************

      y(0) = 0.
      y_meter(0) = 0.
      dist = maxheight/float(nzmax)
      do 135 kz=1,nzmax
        y_meter(kz)= (y(kz-1)+dist)*1000.
        y(kz)= y(kz-1)+dist
135    continue


      maxzgrid=nzmax
      minzgrid=0

      x(num)=seconds(num)/3600.
      topo(num)=heightnn(pktx,pkty,0)

      do lev = 1,maxzgrid
        do kz=1,numzgrid
         if (heightnn(pktx,pkty,kz).gt.y_meter(lev)) goto 223
        end do
 223    continue
        kz=min(kz,numzgrid)
        z(num,lev)= gridtotal(pktx,pkty,kz,lspec)
        if (heightnn(pktx,pkty,0).ge.y_meter(lev)) then    !below topo
         z(num,lev)=-1.0
        endif
        if (z(num,lev).gt.0.) then
         izaehl=izaehl+1
         zmean=zmean+z(num,lev)
c    +          *(heightnn(pktx,pkty,kz)-heightnn(pktx,pkty,kz-1))
        endif
        if (z(num,lev).ge.0..and.z(num,lev).lt.zmin) zmin=z(num,lev)
        if (z(num,lev).gt.zmax) zmax=z(num,lev)
      end do  !(lev loop)
      z(num,0) = -1.0

c end loop over measurement data
********************************

      end do   ! num loop

c determine zmean and contour levels
************************************
      
c     zmean=zmean/float(numb)
c    +(outheight(maxzgrid)-outheight(minzgrid))
      zmean=zmean/izaehl
c     write(*,*) 'zmean=',zmean

      zlev(0) = 0.
      zmax2= zmax   !save zmax field value
      if (zoverrulemax.gt.0.) zmax=zoverrulemax
      zlev(numcon) = zmax
      do 50 i=1,numcon
         lfin(i)=i+2
         if (fixcont.eq.0) then
            zlev(i)=zlev(0)+float(i)/float(numcon)*zmax
         endif
50    continue
      zlev(numcon) = max(zmax,zlev(numcon-1))


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

      call arinam(map,lmap)
      plim1=tstart/3600.0
      plim2=tend/3600.0
      plim3=minheight
      plim4 = maxheight
      call set(0.11,0.95,0.22,0.8,plim1,plim2,plim3,plim4,1)
      call cpseti('SET',0)
      call cpsetr('XC1',plim1)
      call cpsetr('XCM',plim2)
      call cpsetr('YC1',y(0))
      call cpsetr('YCN',y(maxzgrid))
c     call cpsetr('YC1',plim3)
c     call cpsetr('YCN',plim4)
      call cpseti('CLS - Contour Level Selection',0)
      call cpseti('NCL - Number of Contour Levels',numcon)

      call cpseti('RWC - Real Workspace for Contours',500)

      call cpsetr('SPV',-1.0)

      call cpsetc('ILT - Information Label Text String', ' ')
      call cpsetc('HLT - High/Low Label Text Strings', ' '' ')
      call cpseti('LLP - Line Labels', 2)
      if (do_linlab.eq.1.and.do_cont.eq.1)  hlab=3
      if (do_linlab.eq.1.and.do_cont.ne.1)  hlab=2
      if (do_linlab.ne.1.and.do_cont.eq.1)  hlab=1
      if (do_linlab.ne.1.and.do_cont.ne.1)  hlab=0
      do 11 i=1,numcon
        call cpseti('PAI - Parameter Array Index',i)
        call cpsetr('CLV - CONTOUR LEVEL VALUES',zlev(i))
        call cpseti('CLU - Contour Level Use Flag',hlab)
11    continue

c     write(*,*) 'vor cpsps2'
c     call cpseti('ZDS - ZDAT Dimension Selector',1)
c     CALL CPSPS2(X,Y,Z,ntmax,numb,maxzgrid+1,RWRK,LRWK,IWRK,LIWK,
c    +ZREG,LZRG)
c     write(*,*) 'nach cpsps2'
c     write(*,*) lzrg

c     write(*,*) 'vor cprect'
      call cprect(z,ntmax,numb,maxzgrid+1,rwrk,lrwk,iwrk,liwk)
      call cppkcl(z,rwrk,iwrk)
      call cpclam(z, rwrk, iwrk, map)
c     write(*,*) 'nach cprect'


C
C Draw Contours
C

c     call cpclam (zreg, rwrk, iwrk, map)

C Fillstyle solid
      call gsfais(1)
C Draw areas
      if (do_fill.eq.1) then
        call arscam(map, xwrk, ywrk, nwrk, iarea, igrp, nogrps, fill)
      endif
 
C Add Contour-lines, line-labels, surface-line
c     call cpcldm(zreg,rwrk,iwrk,map,cpdrpl)
c     call cpcldr(zreg,rwrk,iwrk)
      call cpcldm(z,rwrk,iwrk,map,cpdrpl)
      call cpcldr(z,rwrk,iwrk)
 
C Labeling: Contour Line Labels
c     call cplbdr(zreg,rwrk,iwrk)
      call cplbdr(z,rwrk,iwrk)

C mark gridpoints
      if (do_gridpts.eq.1) then
         do 45 I=1,it
            do 45 j = minzgrid, maxzgrid
         call points (x(i), y(j), 1, -4, 0)
45       continue
      endif

C surface-line and aircraft track
*********************************

        call setusv('LW',3000)
        call gsplci(0)             ! Plot bold white line
        do 53 i=2,numb
          xra(1)=seconds(i-1)/3600.
          xra(2)=seconds(i)/3600.
          yra(1)=alt(i-1)/1000.
          yra(2)=alt(i)/1000.
53        call curved(xra,yra,2)

        call setusv('LW',2000)
        call gsplci(15)             ! Plot coloured line
        do 57 i=2,numb
          xra(1)=seconds(i-1)/3600.
          xra(2)=seconds(i)/3600.
          yra(1)=alt(i-1)/1000.
          yra(2)=alt(i)/1000.
57        call curved(xra,yra,2)

        call gsplci(1)             ! Plot black line
        call setusv('LW',1000)

      if (do_surface.eq.1) then
        call setusv('LW',5000)
        do 58 i=2,numb
          xra(1)=seconds(i-1)/3600.
          xra(2)=seconds(i)/3600.
          yra(1)=topo(i-1)/1000.
          yra(2)=topo(i)/1000.
58        call curved(xra,yra,2)
        call setusv('LW',1000)

c        call cpseti('PAI - Parameter Array Index',-2)
c        call cpseti('CLU - Contour Level Use Flag',1)
      endif

C Axes
C*****
      call labmod('(f5.2)','(f4.1)',0,0,20,20,20,20,0)
      call gridal(6,6,6,6,1,1,5,0.0,0.0)


C Draw some lines marking the lidar data
****************************************

c     call gslwsc(3.)                    ! set line width
c     call line(199.00,1.1,229.00,1.1)
c     call line(199.00,1.1,199.00,plim4)
c     call line(229.00,1.1,229.00,plim4)
c     call line(199.00,plim4,229.00,plim4)

c     call gslwsc(4.)                    ! set line width
c     call line(130.18,y(0),156.18,y(0))

C Label-bar
C**********
      if (do_lblbar.eq.1) then
      if (((zlev(numcon)-zlev(0)).lt.10.).or.(zlev(numcon).lt.10.).or.
     +(zlev(numcon).gt.99999.)) then
         do 51 i=0,numcon
c1          write(lbls(i),'(pe9.3)') zlev(i)
51          write(lbls(i),'(f6.1)') zlev(i)
      else
         do 52 i=0,numcon
c52         write(lbls(i),*) nint(zlev(i))
 52         write(lbls(i),'(f6.1)') zlev(i)
      endif

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

      xlic1=0.11
      xlic2=0.95
      ylic1=0.07
      ylic2=0.12
      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.or.mod(i,ifact).eq.0.or.i.eq.numcon) then
        call wtstr(xlp,0.05,lbls(i),20,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 setusv('LW',1000)

C Unit-Label for Labelbar
      i1=index(unitlabel(1:10),' ')
      if (i1.eq.0) then
         i1=10
      else
         i1=i1-1
      endif
      call wtstr((xlic2-xlic1)/2+xlic1,0.025,unitlabel(1:i1),20,0,0)
      call wtstr(0.53,0.15,'Time (h)',20,0,0)
      call wtstr(0.01,0.5,'Altitude (km asl)',20,90,0)

      endif


C Draw the release date and time and current time
*************************************************
      iday=(tstart/3600)/24
      iiday=nint(iday)
      ihour=(tstart/3600)-(iday*24) 
      iihour=nint(ihour)
      i2day=(tend/3600)/24
      ii2day=nint(i2day)
      i2hour=(tend/3600)-(i2day*24) 
      ii2hour=nint(i2hour)
      write(label,900) ibdate,ibtime
c     write(label,900) 20050208,ibtime
900   format('Date: ',i8,'  ',i6,'##')
      i1=index(label,'##')
      if (i1.eq.0) then
         i1=21
      else
         i1=i1-1
      endif
      call wtstr(.5,.93,label(1:i1),14,0,0)


      write(label,901) iiday,iihour,ii2day,ii2hour
901   format('Start: +',i3,' d ',i3,' h   End: +',i3,' d ',i3,' h')
      call wtstr(.5,.90,label(1:39),14,0,0)


      write(latstr,'(SP,f5.1)') pktfi
      i1=index(latstr,'+')
      if (i1.eq.0) then
         i1=index(latstr,'-')
         if (i1.eq.0) then
            i1=1
         endif
      else
         i1=i1+1
      endif

      write(lonstr,'(SP,f6.1)') pktlambda
      i2=index(lonstr,'+')
      if (i2.eq.0) then
         i2=index(lonstr,'-')
         if (i2.eq.0) then
            i2=1
         endif
      else
         i2=i2+1
      endif
c     label='Position (lat/lon): '//latstr(i1: )//' / '//lonstr(i2: )//
      label='       Tracer profile '
c     i1=index(label, '##')
c     if (i1.eq.0) then
         i1=38
c     else
c        i1=i1-1
c     endif
cc    call wtstr(.5,.97,label(1:i1),14,0,0)
      call wtstr(.5,.97,label(1:i1),14,0,0)

      if(((zmax2-zmin).lt.10.).or.(zmax2.lt.10.).
     +or.(zmax2.gt.99999.)) then
	    write(maxvalue,'(1pe10.3)') zmax2
	    write(minvalue,'(1pe10.3)') zmin
	    write(meanvalue,'(1pe10.3)') zmean
      else
            write(maxvalue,'(f10.3)') zmax2
            write(minvalue,'(f10.3)') zmin
            write(meanvalue,'(f10.3)') zmean
      endif

      call set(0.,1.,0.,1.,0.,1.,0.,1.,1)
      call wtstr(.5,.87,'   Mean value: '//meanvalue,13,0,0)
      call wtstr(.5,.85,'Maximum value: '//maxvalue,13,0,0)
      call wtstr(.5,.83,'Minimum value: '//minvalue,13,0,0)



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

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


      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+1)
         CALL GFA(NWRK,XWRK,YWRK)
      ENDIF
C
C Otherwise, do nothing
C
      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
         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

      FUNCTION juldate(YYYYMMDD,HHMISS)
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
*                                                                             *
*     Calculates the Julian date                                              *
*                                                                             *
*     AUTHOR: Andreas Stohl (15 October 1993)                                 *
*                                                                             *
*     Variables:                                                              *
*     DD             Day                                                      *
*     HH             Hour                                                     *
*     HHMISS         Hour, minute + second                                    *
*     JA,JM,JY       help variables                                           *
*     JULDATE        Julian Date                                              *
*     JULDAY         help variable                                            *
*     MI             Minute                                                   *
*     MM             Month                                                    *
*     SS             Second                                                   *
*     YYYY           Year                                                     *
*     YYYYMMDDHH     Date and Time                                            *
*                                                                             *
*     Constants:                                                              *
*     IGREG          help constant                                            *
*                                                                             *
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 

      IMPLICIT NONE

      INTEGER YYYYMMDD,YYYY,MM,DD,HH,MI,SS,HHMISS
      INTEGER JULDAY,JY,JM,JA,IGREG
      DOUBLE PRECISION juldate
      PARAMETER (IGREG=15+31*(10+12*1582))

      YYYY=YYYYMMDD/10000
      MM=(YYYYMMDD-10000*YYYY)/100
      DD=YYYYMMDD-10000*YYYY-100*MM
      HH=HHMISS/10000
      MI=(HHMISS-10000*HH)/100
      SS=HHMISS-10000*HH-100*MI

      IF (YYYY.EQ.0) PAUSE 'There is no Year Zero.'
      IF (YYYY.LT.0) YYYY=YYYY+1
      IF (MM.GT.2) THEN
        JY=YYYY
        JM=MM+1
      ELSE
        JY=YYYY-1
        JM=MM+13
      ENDIF
      JULDAY=INT(365.25*JY)+INT(30.6001*JM)+DD+1720995
      IF (DD+31*(MM+12*YYYY).GE.IGREG) THEN
        JA=INT(0.01*JY)
        JULDAY=JULDAY+2-JA+INT(0.25*JA)
      ENDIF

      juldate=DBLE(FLOAT(JULDAY))+DBLE(FLOAT(HH)/24.)+
     +DBLE(FLOAT(MI)/1440.)+DBLE(FLOAT(SS)/86400.)

      RETURN
      END



      SUBROUTINE caldate(JULDATE,YYYYMMDD,HHMISS)
c                           i       o       o
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
*                                                                             *
*     Calculates the Gregorian date from the Julian date                      *
*                                                                             *
*     AUTHOR: Andreas Stohl (21 January 1994), adapted from Numerical Recipes *
*                                                                             *
*     Variables:                                                              *
*     DD             Day                                                      *
*     HH             Hour                                                     *
*     HHMISS         Hour, Minute, Second                                     *
*     JA,JB,JC,JD,JE help variables                                           *
*     JALPHA         help variable                                            *
*     JULDATE        Julian Date                                              *
*     JULDAY         help variable                                            *
*     MI             Minute                                                   *
*     MM             Month                                                    *
*     SS             Seconds                                                  *
*     YYYY           Year                                                     *
*     YYYYMMDD       Year, Month, Day                                         *
*                                                                             *
*     Constants:                                                              *
*     IGREG          help constant                                            *
*                                                                             *
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 

      IMPLICIT NONE

      INTEGER YYYYMMDD,YYYY,MM,DD,HHMISS,HH,MI,SS
      INTEGER JULDAY,JA,JB,JC,JD,JE,IGREG,JALPHA
      DOUBLE PRECISION JULDATE
      PARAMETER (IGREG=2299161)

      JULDAY=INT(JULDATE)
      IF(JULDAY.GE.IGREG)THEN
        JALPHA=INT(((JULDAY-1867216)-0.25)/36524.25)
        JA=JULDAY+1+JALPHA-INT(0.25*JALPHA)
      ELSE
        JA=JULDAY
      ENDIF
      JB=JA+1524
      JC=INT(6680.+((JB-2439870)-122.1)/365.25)
      JD=365*JC+INT(0.25*JC)
      JE=INT((JB-JD)/30.6001)
      DD=JB-JD-INT(30.6001*JE)
      MM=JE-1
      IF (MM.GT.12) MM=MM-12
      YYYY=JC-4715
      IF (MM.GT.2) YYYY=YYYY-1
      IF (YYYY.LE.0) YYYY=YYYY-1

      YYYYMMDD=10000*YYYY+100*MM+DD
      HH=INT(24.*(JULDATE-FLOAT(JULDAY)))
      MI=INT(1440.*(JULDATE-FLOAT(JULDAY))-60.*FLOAT(HH))
      SS=NINT(86400.*(JULDATE-FLOAT(JULDAY))-3600.*FLOAT(HH))
     +-60.*FLOAT(MI)
      IF (SS.EQ.60) THEN  ! 60 seconds = 1 minute
        SS=0
        MI=MI+1
      ENDIF
      IF (MI.EQ.60) THEN
        MI=0
        HH=HH+1
      ENDIF
      HHMISS=10000*HH+100*MI+SS

      RETURN
      END
