      program flexpart44_dataextract

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

      integer inday(maxdates),intime(maxdates)
      double precision julin(maxdates)


      real outheight(nzmax)
      integer ireleasestart(maxpoint),ireleaseend(maxpoint)
      integer itime
      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 oro(nxmax,nymax)
      character compoint(maxpoint)*10,species(maxspec)*7
      character*120 dirname,filename,datesname,datafilename,outfilename
      character fmt*23,aday*8,atime*6
      character alspec*10


      parameter(maxdata=50000,maxcolumns=10)
      integer jjjjmmdd(maxdata),ihmmss(maxdata)
      real xlon(maxdata),ylat(maxdata),alt(maxdata)
      real data(maxdata,maxcolumns),tracer(maxdata,maxspec,maxageclass)
      double precision juldate,juldata(maxdata)





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

      open(10,file='CONTROL')
      read(10,'(a120)') dirname
      read(10,'(a120)') datesname
      read(10,*) lconcormixratio
      read(10,*) decayconstant
      read(10,'(a10)') alspec               ! kind of species to be plotted
      read(10,'(a120)') datafilename
      read(10,'(a120)') outfilename
      read(10,*) ncolumns
      read(10,*) iswitch
      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     write(*,*) dirname 
c     write(*,*) datesname
c     write(*,*) datafilename
c     write(*,*) outfilename


C Read datafile
***************

      open(20,file=datafilename)
      n=1
200   read(20,*,end=199) jjjjmmdd(n),ihmmss(n),ylat(n),xlon(n),alt(n),
     +(data(n,j),j=1,ncolumns)
ccc  +(data(n,j),j=1,ncolumns-1)
      juldata(n)=juldate(jjjjmmdd(n),ihmmss(n))
      do 11 nage=1,maxageclass
        do 11 k=1,maxspec
11        tracer(n,k,nage)=-99.99
      n=n+1
      goto 200
199   continue
      close(20)

      numdata=n-1
      nummin=1

c      write(*,*) numdata, nummin


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

c      write(*,*) datesname

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


C Read header for 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) iiiimmdd,iimmss
      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)

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

      len=index(dirname,' ')-1
      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 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 Open output file
******************

      open(30,file=outfilename)
      write(30,*) nageclass,(lage(i),i=1,nageclass)
      write(30,*) lspec



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

      do 150 num=nummin,numdata


C Find appropriate model output file
************************************

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

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

            do 38 k=1,nspec

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

32              continue

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

                read(10) lval
                if (lval.eq.1) then           !sparse matrix
36                read(10) 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)
c31           read(10) (drygrid(ix,jy,k,nage),yuncert,jy=1,numygrid)
                endif

35              continue

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

                read(10) lval
                if (lval.eq.1) then           !sparse matrix
40                read(10) n,xx
                  if (n.lt.0) goto 78
                  kz=n/(numxgrid*numygrid)
                  jy=(n-kz*numxgrid*numygrid)/numxgrid
                  ix=n-numxgrid*numygrid*kz-numxgrid*jy
                  ix=ix+1
                  jy=jy+1
                  grid(ix,jy,kz,k,nage)=xx
                  goto 40
                else
                  do 39 kz=1,numzgrid
                    do 39 ix=1,numxgrid
 39                    read(10) (grid(ix,jy,kz,k,nage),
     +                jy=1,numygrid)
c39                    read(10) (grid(ix,jy,kz,k,nage),yuncert,
c    +                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)
41                    continue


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

38            continue
          endif

C Assign tracer data to measurement locations
*********************************************

          indx=int((xlon(num)-outlon0)/dxout)+1
          indy=int((ylat(num)-outlat0)/dyout)+1
          if ((indx.ge.1).and.(indy.ge.1).and.(indx.le.numxgrid).and.
     +    (indy.le.numygrid)) then
            do 82 l=1,numzgrid
              if ((outheight(l)+oro(indx,indy)).gt.alt(num)) then
                do 182 k=1,nspec
182               tracer(num,k,1)=grid(indx,indy,l,k,1)
                do 79 nage=2,nageclass
                  do 79 k=1,nspec
79                  tracer(num,k,nage)=tracer(num,k,nage-1)+
     +              grid(indx,indy,l,k,nage)
                goto 83
              endif
82            continue

83          continue
          endif

          endif

160       continue


C Dump tracer data to file
**************************

        fmt='(i8,1x,i6,  (1x,e13.5))'
        write(fmt(11:12),'(i2.2)') 3+ncolumns+nspec*nageclass

        write(30,fmt) jjjjmmdd(num),ihmmss(num),ylat(num),xlon(num),
     +  alt(num),(data(num,j),j=1,ncolumns),
ccc  +  alt(num),(data(num,j),j=1,ncolumns-1),
     +  (tracer(num,lspec,j),j=1,nageclass)

         
150     continue

      close(30)

      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
      REAL*8 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
