C Plot FLEXPART tracer age spectra along aircraft flight leg, or similar
************************************************************************

      program agespectra
      parameter(maxdata=60000,maxspec=3,maxageclass=32)
      parameter(nageclassmax=20)
      integer jjjjmmdd(maxdata),ihmmss(maxdata),lage(0:maxageclass)
      real tracer(maxdata,maxspec,maxageclass),tracegas(maxdata)
      real xlon(maxdata),ylat(maxdata),alt(maxdata)
      real hour(maxdata),yt(maxdata),hourtracegas(maxdata)
      real dummy(50)
      double precision juldate
      character filename*120,filename2*120,fmt*23
      

      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=49000,liwk=49000,lmap=900000)
      parameter(nwrk=49000,isiz=5)
      real xwrk(nwrk),ywrk(nwrk),rwrk(lrwk),iwrk(liwk)
      integer map(lmap),iarea(isiz),igrp(isiz)
      integer dst(8),ind(12)
      real yra(5),xra(5)
      character lbls(0:maxcon)*2,date*8,panel(5)*2
      character*10 label1,label2,unitlabel1,unitlabel2
      character title*80

      data panel/'a)','b)','c)','d)','e)'/




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

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

      call setusv('LW',1000)
      call gstxfp(-13,2)
      call pcseti('FN',21)



      plotheightmin=0.16
      plotheight=0.95-plotheightmin

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

      open(10,file='CONTROL_agespectra',status='old')
      read(10,'(a)') title
      read(10,*) overrulemax
      read(10,*) lcolor
      read(10,*) numplots
      read(10,*)
      panelheight=plotheight/float(numplots)
      do 400 nplot=1,numplots
        read(10,*)
        read(10,'(a120)') filename
        read(10,*) kspec
        read(10,*) scaleconc
        read(10,'(a)') label1
        read(10,'(a)') unitlabel1
        read(10,'(a)') label2
        read(10,'(a)') unitlabel2


C Read tracer data
******************

        rmaxalt=0.

c       write(*,'(a)') filename
        open(20,file=filename,status='old')
        read(20,*) nageclass,(lage(i),i=1,nageclass)
        read(20,*) nspec
        fmt='(i8,1x,i6,  (1x,e13.5))'
        write(fmt(11:12),'(i2.2)') 3+nspec*nageclass
        num=1
100       read(20,fmt,end=99) jjjjmmdd(num),ihmmss(num),ylat(num),
     +    xlon(num),alt(num),((tracer(num,k,j),j=1,nageclass),k=1,nspec)
          if (alt(num).gt.rmaxalt) rmaxalt=alt(num)

          num=num+1
        goto 100

99      num=num-1
        close(20)

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

        nageplot=min(nageclass,nageclassmax)
        xc=1.0
        xb=0.0
        xg=1.0
        xl=0.0
        xh=0.84
        jr=0
        jn=1
        do 10 i=1,nageplot+2
          call colspec(lcolor,i-1,nageplot+2,xl,xh,jr,jn,xc,xb,xg,0,i+1)
10      continue

C Determine maximum values
**************************

        ymax=0.
        ymin=0.
        timespan=juldate(jjjjmmdd(num),ihmmss(num))-
     +  juldate(jjjjmmdd(1),ihmmss(1))
        do 20 i=1,num 
          do 30 k=1,nspec
            do 30 j=1,nageclass
30            tracer(i,k,j)=tracer(i,k,j)*scaleconc    ! scale the tracer concentrations
          if (tracer(i,kspec,nageclass).gt.ymax)
     +    ymax=tracer(i,kspec,nageclass)
          hour(i)=sngl(juldate(jjjjmmdd(i),ihmmss(i))-
     +    juldate(jjjjmmdd(1),0))
          if (timespan.lt.1.) hour(i)=hour(i)*24.
20        continue
        ymax=1.01*ymax

        ymax=max(ymax,4.)

      if (overrulemax.gt.0.) then
        ymax=overrulemax
      endif

c             write(*,*) 'ymax ',ymax


C Determine optimal labels
**************************

        if (hour(num)-hour(1).gt.8) then   ! for long flights: tick marks every hour
          xmin=float(int(hour(1)))
          xmax=float(int(hour(num)+1))
          jxen=nint(xmax-xmin)
          indlab=1
        else                               ! for short flights: tick marks every 1/2 hour
          xmin=float(int(hour(1)*2))/2.
          xmax=float(int(hour(num)*2+1))/2.
          jxen=nint(2*(xmax-xmin))
          indlab=2
        endif

        range=ymax-ymin
        aux=float(int(log10(range)))
        if (range/10**aux.le.3) aux=aux-1
        if (range.le.6.*10**aux) then
          jyen=nint(range)/10**aux+1
          range=10.**aux*float(jyen)
        else if (range.le.12.*10**aux) then
          jyen=nint(range)/2/10**aux+1
          range=2.*10.**aux*float(jyen)
        else
          jyen=nint(range)/4/10**aux+1
          range=4.*10.**aux*float(jyen)
        endif
        tickint=range/float(jyen)
        naux=int(ymin/tickint)
        ymin=naux*tickint
        prelimmax=ymin+jyen*tickint
        if (prelimmax.lt.ymax) then
          ymax=prelimmax+tickint
          jyen=jyen+1
        else
          ymax=prelimmax
        endif

        deltat=hour(2)-hour(1)


C Plot age spectrum
*******************

        call set(0.12,0.91,plotheightmin+float(nplot-1)*panelheight,
     +  plotheightmin+float(nplot)*panelheight,xmin,xmax,ymin,ymax,1)

        do 341 n=nageplot,1,-1
          icolor=int(float(n)/float(nageplot)*float(nageplot))+3
          call gsplci(icolor)
          do 342 i=1,num
            yt(i)=tracer(i,kspec,n)
            yra(1)=0.
            yra(2)=tracer(i,kspec,n)
            yra(3)=tracer(i,kspec,n)
            yra(4)=0.
            yra(5)=0.
            xra(3)=hour(i)+deltat/2.
            xra(4)=hour(i)+deltat/2.
            if (i.eq.1) then
              xra(1)=hour(i)-deltat
              xra(2)=hour(i)-deltat
              xra(5)=hour(i)-deltat
            else
              xra(1)=hour(i-1)
              xra(2)=hour(i-1)
              xra(5)=hour(i-1)
            endif
            call sfsgfa(xra,yra,4,dst,nst,ind,nnd,icolor)
342         continue
            call gsplci(1)
341       continue

        call gsplci(1)
        if (indlab.eq.1) then
          call labmod('(I3)','(f6.1)',0,0,14,14,14,14,0)
        else
          call labmod('(f4.1)','(f6.1)',0,0,14,14,14,14,0)
        endif
        call setusv('LW',1000)
        CALL TICK4 (10,6,10,6)
        if (nplot.ne.1) then
          CALL GRIDAL (jxen,1,jyen,1,0,1,10,xmin,ymin)
          CALL GRIDAL (jxen,1,jyen,1,0,1,6,xmin,ymax)
        else
          CALL GRIDAL (jxen,1,jyen,1,1,1,10,xmin,ymin)
          CALL GRIDAL (jxen,1,jyen,1,1,1,6,xmin,ymax)
        endif


C Read from CONTROL file which trace gas species shall be plotted
*****************************************************************

        read(10,*) ntotnumbplot
        do 400 nn=1,ntotnumbplot
          read(10,'(a120)') filename2
          read(10,*) ntotfile
          read(10,*) nnumbplot
          read(10,*) scaletracegas


**************************************************************************************
C Read trace gas measurements; if altitude is selected, take altitude from tracer file
**************************************************************************************

        if (filename2(1:8).eq.'Pressure') then
          numtracegas=num
          do 333 i=1,num
            hourtracegas(i)=hour(i)
!!! Use standard atmosphere to convert to altitude
            if (alt(i).gt.899.) then
              alt(i)=max(1013.-alt(i),0.)/111.
            else if (alt(i).gt.795) then
              alt(i)=1.+(899.-alt(i))/101.
            else if (alt(i).gt.701) then
              alt(i)=2.+(795.-alt(i))/94.
            else if (alt(i).gt.617) then
              alt(i)=3.+(701.-alt(i))/84. 
            else if (alt(i).gt.540) then
              alt(i)=4.+(617.-alt(i))/77. 
            else if (alt(i).gt.472) then
              alt(i)=5.+(540.-alt(i))/68. 
            else if (alt(i).gt.356) then
              alt(i)=6.+(472.-alt(i))/58.
            else if (alt(i).gt.265) then
              alt(i)=8.+(356.-alt(i))/45.5
            else if (alt(i).gt.194) then
              alt(i)=10.+(265.-alt(i))/35.5
            else
              alt(i)=12.+(194.-alt(i))/24.
            endif
333         tracegas(i)=alt(i)
        else if (filename2(1:8).eq.'Altitude') then
          numtracegas=num
          do 33 i=1,num
            hourtracegas(i)=hour(i)
            if (rmaxalt.lt.1000.) then
             tracegas(i)=alt(i)
            else
             tracegas(i)=alt(i)/1000.
            endif
33        continue
        else
          open(20,file=filename2,status='old')
          read(20,*) n
          do 34 i=1,n-1
34          read(20,*)
          numtracegas=1
!!!! CHANGE seconds from integer to real format
35          read(20,*,end=36) isec,(dummy(l),l=1,ntotfile)
c5          read(20,*,end=36) sec,(dummy(l),l=1,ntotfile)
            tracegas(numtracegas)=dummy(nnumbplot)*scaletracegas
            hourtracegas(numtracegas)=float(isec)/3600.
c           hourtracegas(numtracegas)=sec/3600.
C Check for missing values
            if (tracegas(numtracegas).gt.999999.)
     +      tracegas(numtracegas)=-99999.
            numtracegas=numtracegas+1
            goto 35
36        numtracegas=numtracegas-1
          close(20)
        endif


C Determine maximum and minimum values (only for first species to be plotted)
*****************************************************************************

        if (nn.eq.1) then
          tracegasmax=0.
          tracegasmin=9.e29
          do 21 i=1,numtracegas 
            if (tracegas(i).gt.tracegasmax) tracegasmax=tracegas(i)
            if ((tracegas(i).lt.tracegasmin).and.(tracegas(i).gt.0.))
     +      tracegasmin=tracegas(i)
21          continue
          tracegasmax=tracegasmax+0.01*(tracegasmax-tracegasmin)
          tracegasmin=tracegasmin-0.01*(tracegasmax-tracegasmin)

C CHANGE: Use this for scaling
           if (numtracegas.ne.num) then
             tracegasmin=80.
             tracegasmax=150.
           endif


C Determine optimal labels
**************************

          goto 4567

          tracegasrange=tracegasmax-tracegasmin
          aux=float(int(log10(tracegasrange)))
          if (tracegasrange/10**aux.le.3) aux=aux-1
          if (tracegasrange.le.6.*10**aux) then
            jyen2=nint(tracegasrange)/10**aux+1
            tracegasrange=10.**aux*float(jyen2)
          else if (tracegasrange.le.12.*10**aux) then
            jyen2=nint(tracegasrange)/2/10**aux+1
            tracegasrange=2.*10.**aux*float(jyen2)
          else
            jyen2=nint(tracegasrange)/4/10**aux+1
            tracegasrange=4.*10.**aux*float(jyen2)
          endif
          tickint=tracegasrange/float(jyen2)
          naux=int(tracegasmin/tickint)
          tracegasmin=naux*tickint
          prelimmax=tracegasmin+jyen2*tickint
          if (prelimmax.lt.tracegasmax) then
            tracegasmax=prelimmax+tickint
            jyen2=jyen2+1
          else
            tracegasmax=prelimmax
          endif
        endif

4567    continue

            write(*,*) tracegasmin,tracegasmax

C Plot additional trace gas data as lines
*****************************************

        call set(0.12,0.91,plotheightmin+float(nplot-1)*panelheight,
     +  plotheightmin+float(nplot)*panelheight,xmin,xmax,tracegasmin,
     +  tracegasmax,1)

        call setusv('LW',4000)
        call gsplci(0)             ! Plot bold white line
        do 77 i=2,numtracegas
          if ((tracegas(i-1).gt.0.).and.(tracegas(i).gt.0.)) then  ! do not plot missing data
            xra(1)=hourtracegas(i-1)
            xra(2)=hourtracegas(i)
            yra(1)=tracegas(i-1)
            yra(2)=tracegas(i)
            call curved(xra,yra,2)
          endif
77        continue
        call setusv('LW',2000)
        call gsplci(2*nn-1)          ! Plot color line
        do 78 i=2,numtracegas
          if ((tracegas(i-1).gt.0.).and.(tracegas(i).gt.0.)) then  ! do not plot missing data
            xra(1)=hourtracegas(i-1)
            xra(2)=hourtracegas(i)
            yra(1)=tracegas(i-1)
            yra(2)=tracegas(i)
            call curved(xra,yra,2)
          endif
78        continue



        if (indlab.eq.1) then
          call labmod('(I3)','(I4)',0,0,14,14,-6,14,0)
        else
          call labmod('(f4.1)','(I4)',0,0,14,14,-6,14,0)
        endif
        call gsplci(1)
        call setusv('LW',1000)
        CALL TICK4 (10,6,-10,6)
        if (nplot.ne.1) then
          CALL GRIDAL (jxen,1,jyen2,1,0,1,10,xmax,tracegasmin)
        else
          CALL GRIDAL (jxen,1,jyen2,1,1,1,10,xmax,tracegasmin)
        endif


C Write axis labels
*******************

        call set(0.,1.,0.,1.,0.,1.,0.,1.,1)
        ypos=plotheightmin+(float(nplot-1)+.55)*panelheight

        call plchhq(0.04,ypos+0.35*panelheight,
     +  panel(numplots-nplot+1),.014,0.,0)
    
        call plchhq(0.04,ypos-0.12,label1,.013,90.,0)
        call plchhq(0.04,ypos-0.03,unitlabel1,.013,90.,0)
        call plchhq(0.97,ypos-0.08,label2,.013,90.,0)
        call plchhq(0.97,ypos-0.03,unitlabel2,.013,90.,0)

C Hide uppermost label for all but the top plot
        if (nplot.ne.numplots)
     +  call ngdots(0.935,ypos+0.45*panelheight,1,0.04,0)
        if (nplot.ne.numplots)
     +  call ngdots(0.095,ypos+0.45*panelheight,1,0.041,0)

400   continue




***************************************************
C Plot everything not belonging to a specific panel
***************************************************

      write(date,'(i8)') jjjjmmdd(1)
      if (timespan.lt.1.) then
        call plchhq(0.5,0.11,'Hours past 0 UTC on '//date,0.014,0.,0)
      else
        call plchhq(0.5,0.11,'Days past 0 UTC on '//date,0.014,0.,0)
      endif
      call plchhq(0.6,0.97,title,0.016,0.,0)


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

      lage(0)=0
      do 24 i=0,nageplot
24      write(lbls(i),'(i2)') lage(i)/86400


C Number of text labels


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


      xlic1=0.02
      xlic2=0.98
      ylic1=0.060
      ylic2=0.085
      do i=1,nageplot
        xxo=xlic1+(xlic2-xlic1)*float(i-1)/float(nageplot)
        xxpp=xlic1+(xlic2-xlic1)*float(i)/float(nageplot)
        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+3)
        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,nageplot
      xlp=xlic1+(xlic2-xlic1)*float(i)/float(nageplot)
      call plchhq(xlp-0.003,ylic1-.015,lbls(i),0.014,0.,0)
      end do
      call setusv('LW',2000)
      do i=1,nageplot-1
      if (mod(i,1).eq.0) then
        xlp=xlic1+(xlic2-xlic1)*float(i)/float(nageplot)
        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 plchhq(0.5,0.02,'Age (days)',.014,0.,0)

      call frame


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

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


      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

      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

