PROGRAM C2DATEM

!-------------------------------------------------------------------------------
! sample program to read hysplit4 concentration file and merge results
! with measurements in the DATEM file format
!-------------------------------------------------------------------------------
! Last Revised: 15 Jun 2001 - initial version from con2stn
!               31 Jul 2001 - optional diagnostics
!               15 Dec 2003 - more flexible input format
!               04 Mar 2005 - add data rotation option
!               30 Jun 2008 - header record optional
!               01 Aug 2011 - format change to -c option
!               06 Sep 2011 - documentation enhancement
!               04 Nov 2011 - pass supplemental information to output
!               17 Nov 2011 - syntax correction for PGF compiler
!               11 Dec 2011 - further mdification to supplemental header
!               18 Oct 2012 - expanded format for output concentration field
!-------------------------------------------------------------------------------

! output to diagnostic file
  LOGICAL        :: DIAG
  INTEGER        :: HEAD 

! output of supplemental lat/lon information
  LOGICAL        :: INFO

! identification strings
  CHARACTER(500) :: LINE1,LINE2
  CHARACTER(80)  :: FNAME,LABEL         
  CHARACTER(4)   :: PTYPE,MODEL
  INTEGER(4)     :: K,KLEN

! master measured concentration array
  INTEGER,      ALLOCATABLE :: SITE(:)
  INTEGER,      ALLOCATABLE :: MTIME(:)
  INTEGER,      ALLOCATABLE :: IYR(:),IMO(:),IDA(:),IHR(:),IDUR(:)
  REAL,         ALLOCATABLE :: CNTR(:),CONC(:),SLAT(:),SLON(:)
  REAL,         ALLOCATABLE :: PLAT(:),PLON(:)

! model results concentration arrays
  REAL,         ALLOCATABLE :: CSUM  (:,:)
  CHARACTER(4), ALLOCATABLE :: IDENT (:) 
  INTEGER,      ALLOCATABLE :: HEIGHT(:)

! special concentration packing variables
  INTEGER(2)    :: ip,jp
  INTEGER       :: cpack, nsta 

! optional data rotation for plume matching statistics
  REAL          :: ANGLE, DELTA  ! rotation angle 
  REAL          :: RLAT, RLON    ! rotation origin point
  REAL          :: XS,YS         ! distance components
  REAL          :: ANG1,ANG2     ! angle from abcissa
  REAL          :: DIST          ! distance from rotation point

!--------------------------------------------------------------------
! configure arrays and extraction method
!--------------------------------------------------------------------

! read any command line arguments and set default values
  CALL SETARG(CFACT,LINEAR,JLEV,JPOL,DIAG,INFO,ANGLE,RLAT,RLON,HEAD)

! diagnostic messages output file
  IF(DIAG) OPEN(40,FILE='c2datem.txt')

! open concentration file to determine array allocations
  CALL CONSET(diag,nloc,nlon,nlat,nlvl,ntyp)
  ALLOCATE (csum(nlon,nlat), STAT=kret)
  ALLOCATE (height(nlvl),    STAT=kret)
  ALLOCATE (ident(ntyp),     STAT=kret)
  IF(DIAG) WRITE(40,*)'Concentration file allocations' 

! open measurement data file to determine array allocations
  CALL SETOBS(diag,numobs)
  ALLOCATE (iyr  (numobs), STAT=kret)
  ALLOCATE (imo  (numobs), STAT=kret)
  ALLOCATE (ida  (numobs), STAT=kret)
  ALLOCATE (ihr  (numobs), STAT=kret)
  ALLOCATE (idur (numobs), STAT=kret)
  ALLOCATE (slat (numobs), STAT=kret)
  ALLOCATE (slon (numobs), STAT=kret)
  ALLOCATE (plat (numobs), STAT=kret)
  ALLOCATE (plon (numobs), STAT=kret)
  ALLOCATE (conc (numobs), STAT=kret)
  ALLOCATE (site (numobs), STAT=kret)
  ALLOCATE (cntr (numobs), STAT=kret)
  ALLOCATE (mtime(numobs), STAT=kret)
  IF(DIAG) WRITE(40,*)'Measurement data allocations'

! initialize calendar functions
  CALL tminit

!--------------------------------------------------------------------
! load measurement data
!--------------------------------------------------------------------

! identification records
  READ(20,'(A80)')LINE1(1:80)
  READ(20,'(A80)')LINE2(1:80)

  DO N=1,NUMOBS
!    observational data in datem format
     READ(20,*)IYR(N),IMO(N),IDA(N),IHR(N),IDUR(N),     &
	           SLAT(N),SLON(N),CONC(N),SITE(N)

!    adjust hour-min
     LHRS=IHR(N)/100
     LMIN=IHR(N)-LHRS*100

!    convert time to minutes for computational purposes
     CALL TM2MIN(IYR(N),IMO(N),IDA(N),LHRS,LMIN,MTIME(N))

!    set concentration to zero for possible averaging (discard meas)
     CONC(N)=0.0
     CNTR(N)=0.0
  END DO
  CLOSE(20)

! save orignal positions
  PLAT=SLAT
  PLON=SLON

! plume rotation (positive angle moves plume clockwise!)
! note that the plume data remain fixed but the samplers are moved
  IF(ANGLE.NE.0.0)THEN
     FACT=COS(RLAT/57.3)
     DELTA=ANGLE/57.3

     DO N=1,NUMOBS
!       convert to polar coordinates
        YS=(SLAT(N)-RLAT)
        XS=(SLON(N)-RLON)*FACT
        DIST=SQRT(XS*XS+YS*YS)

!       apply angular shift
        ANG1=ATAN2(YS,XS)
        ANG2=ANG1+DELTA
        XS=DIST*COS(ANG2)
        YS=DIST*SIN(ANG2)

!       compute new location
        SLAT(N)=YS+RLAT
        SLON(N)=XS/FACT+RLON
     END DO
  END IF

!--------------------------------------------------------------------
! setup concentration grid information from index records
!--------------------------------------------------------------------

! meteo file information and calculation start data
  READ(10)MODEL,KYR,KMO,KDA,KHR,KCH,NLOC,CPACK
  DO N=1,NLOC
     READ(10)IBYR,IBMO,IBDA,IBHR,OLAT,OLON,OLVL
  END DO

! horizontal grid index record
  READ(10) NLAT, NLON, DLAT, DLON, CLAT, CLON

! vertical grid index record
  READ(10) NLVL, (HEIGHT(KK),KK=1,NLVL)

! pollutant identification record
  READ(10) NTYP, (IDENT(KK),KK=1,NTYP)

! origin point and filename for statistics plot
  IF(HEAD.EQ.1)THEN
     IF(INFO)THEN
        INQUIRE(UNIT=10,NAME=LABEL)
        KLEN=INDEX(LABEL,' ')-1
        DO K=KLEN,1,-1
!          replace forward (47) and back (92) slash with the ICHAR equivalent
           IF(ICHAR(LABEL(K:K)).EQ.47.OR.ICHAR(LABEL(K:K)).EQ.92) EXIT
        END DO
        K=MAX(1,K+1)
        KLEN=MIN(KLEN,59)
        FNAME=LABEL(K:KLEN)
        WRITE(30,'(2F10.3,2A)')OLAT,OLON,' ',FNAME
     ELSE
        KLEN=INDEX(LINE1,' ')-1
        WRITE(30,'(2A,F5.1)')LINE1(1:KLEN),' : Rotation=',ANGLE
     END IF
     WRITE(30,'(A)')LINE2(1:80)

! origin point and initial label information for summary file
  ELSEIF(HEAD.EQ.2)THEN
     KLEN=INDEX(LINE1,' ')-1
     IF(INFO)THEN
        WRITE(30,'(2F10.3,2A)')OLAT,OLON,' ',LINE1(1:KLEN)
     ELSE
        WRITE(30,'(2A,F5.1)')LINE1(1:KLEN),' : Rotation=',ANGLE
     END IF
     WRITE(30,'(A)')LINE2(1:80)

! no header information output
  ELSE
     CONTINUE
  END IF

!--------------------------------------------------------------------
! process gridded concentrations by time period
!--------------------------------------------------------------------

  KREC=0
mloop : DO WHILE (KREC.EQ.0)

!    load model average start and stop times
     READ(10,IOSTAT=KREC)KYR,KMO,KDA,KHR,KMN,KFH
     IF(KREC.NE.0)EXIT mloop
     READ(10)            JYR,JMO,JDA,JHR,JMN,JFH

!    load model concentration data
     DO KP=1,NTYP
     DO KL=1,NLVL
!       test if pollutant and level of record matches request
        IF(KP.EQ.JPOL.AND.KL.EQ.JLEV)THEN
           IF(CPACK.EQ.1)THEN
              CSUM = 0.0
              READ(10)PTYPE,LEVEL,NXYP,(IP,JP,CSUM(IP,JP),NP=1,NXYP)
           ELSE
              READ(10)PTYPE,LEVEL,((CSUM(II,JJ),II=1,NLON),JJ=1,NLAT)
           END IF
        ELSE
           READ(10)PTYPE,LEVEL
        END IF
     END DO
     END DO

!    sample duration in hours
     CALL TM2MIN(KYR,KMO,KDA,KHR,KMN,JET1)
     CALL TM2MIN(JYR,JMO,JDA,JHR,JMN,JET2)
     IF(DIAG) WRITE(40,*)'Processing model results: ',KMO,KDA,KHR,KMN

     NSTA=0
!    match observation to model result times
     DO N=1,NUMOBS

!       adjust sample duration times (permits duration = DDHHMM) for
!       samples that are longer than 99 hours
        LDAY=IDUR(N)/10000
        LDHR=(IDUR(N)-LDAY*10000)/100
        LDMN=IDUR(N)-(LDAY*10000+LDHR*100)

!       compute end of sample time
        NTIME=MTIME(N)+(LDAY*24+LDHR)*60+LDMN

!       assume that model durations always <= measured durations           
        IF(NTIME.GT.JET1.AND.MTIME(N).LT.JET2)THEN

!          model result falls within sample start and ending times
!          summation of modeling results into sampling window

!          convert location to grid point
           XI=(SLON(N)-CLON)/DLON+1.0
           YJ=(SLAT(N)-CLAT)/DLAT+1.0
           II=INT(XI)
           JJ=INT(YJ)

!          estimate concentration at measurement location
           IF(II.GE.1.AND.II.LT.NLON.AND.JJ.GE.1.AND.JJ.LE.NLAT)THEN
              NSTA=NSTA+1
              IF(LINEAR.EQ.0)THEN
!                nearest neighbor concentration extraction
                 II=NINT(XI)
                 JJ=NINT(YJ)
                 XCON=CSUM(II,JJ)*CFACT
              ELSE
!                linear interpolation extraction
                 XF=XI-II
                 YF=YJ-JJ
                 CTOP=(CSUM(II+1,JJ+1)-CSUM(II,JJ+1))*XF+CSUM(II,JJ+1)
                 CBOT=(CSUM(II+1,JJ  )-CSUM(II,JJ  ))*XF+CSUM(II,JJ)
                 XCON=((CTOP-CBOT)*YF+CBOT)*CFACT
              END IF
           ELSE
!             at edge or outside of grid
              XCON=-1.0
              IF(DIAG) WRITE(40,*)'Station outside of grid: ',SITE(N)
           END IF

!          sum into the measurement array element
           IF(XCON.GE.0.0)THEN
              CONC(N)=CONC(N)+XCON
              CNTR(N)=CNTR(N)+1.0
           END IF

        END IF

        IF(NTIME.EQ.JET2)THEN
!          when this sample is also the end of the summation
           IF(CNTR(N).GT.0.0)THEN
              CONC(N)=CONC(N)/CNTR(N)
           ELSE
              CONC(N)=-1.0
           END IF                     

!          results in standard datem format
           WRITE(30,'(I4,2I3,I5.4,I7.4,2F10.4,F12.1,I8)')     &
                 iyr(n),imo(n),ida(n),ihr(n),idur(n),         &
                 plat(n),plon(n),conc(n),site(n)

           CONC(N)=0.0
           CNTR(N)=0.0
        END IF

     END DO
     IF(DIAG) WRITE(40,*)'   stations matched: ',NSTA

! model concentrations time loop
  END DO mloop

  CLOSE(30)

END PROGRAM c2datem

!---------------------------------------------------------------------

SUBROUTINE CONSET(diag,nloc,nlon,nlat,nlvl,ntyp)

  LOGICAL      :: diag
  INTEGER(2)   :: ip,jp
  INTEGER      :: cpack, height 
  CHARACTER(4) :: IDENT, PTYPE, MODEL

! meteo file information and calculation start data
  READ(10,IOSTAT=kret)MODEL,IYR,IMO,IDA,IHR,ICH, NLOC, CPACK
  IF(kret.ne.0)CPACK=0   ! old format files don't have cpack

! calculation start for each location
  DO N=1,NLOC
     READ(10)IBYR,IBMO,IBDA,IBHR,OLAT,OLON,OLVL
  END DO

! horizontal grid index record
  READ(10) NLAT, NLON, DLAT, DLON, CLAT, CLON

! vertical grid index record
  READ(10) NLVL, (HEIGHT,KK=1,NLVL)

! pollutant identification record
  READ(10) NTYP, (IDENT,KK=1,NTYP)

  REWIND(10)
  IF(DIAG) WRITE(40,*)'Model grid size: ',NLAT,NLON

END SUBROUTINE conset  

!---------------------------------------------------------------------

SUBROUTINE SETARG(CFACT,LINEAR,JLEV,JPOL,DIAG,INFO,ANGLE,RLAT,RLON,HEAD)

  CHARACTER(80) :: FNAME,LABEL   ! identification strings
  LOGICAL       :: INFO
  LOGICAL       :: FTEST, DIAG
  INTEGER       :: HEAD
  REAL          :: ANGLE         ! rotation angle degrees
  REAL          :: RLAT, RLON    ! rotation origin point

! command line argument list

  KARG=IARGC()

  IF(KARG.EQ.0)THEN
     WRITE(*,*)'Program to read hysplit4 concentration file and merge'
     WRITE(*,*)'results with measurements in the DATEM file format'
     WRITE(*,*)
     WRITE(*,*)'USAGE: c2datem -[arguments]'
     WRITE(*,*)'-c[input to output concentration multiplier]'
     WRITE(*,*)'-d[write to diagnostic file]'
     WRITE(*,*)'-h[header info: 2=input text, (1)=file name, 0=skip]'
     WRITE(*,*)'-i[input concentration file name]'
     WRITE(*,*)'-m[measurement data file name]'
     WRITE(*,*)'-o[output concentration file name]'
     WRITE(*,*)'-p[pollutant index select for multiple species]'
     WRITE(*,*)'-r[rotation angle:latitude:longitude]'
     WRITE(*,*)'-s[supplemental lat-lon information to output]'
     WRITE(*,*)'-x[n(neighbor) or i(interpolation)]'
     WRITE(*,*)'-z[level index select for multiple levels]'
     WRITE(*,*)'ENTER to continue ...'
     READ(*,*)
     STOP
  END IF

! set default values

  DIAG=.FALSE.
  INFO=.FALSE.
  CFACT=1.0
  LINEAR=0
  JLEV=1
  JPOL=1
  ANGLE=0.0
  HEAD=1

! go through each argument
  DO WHILE (KARG.GT.0)

     CALL GETARG(KARG,LABEL)
     SELECT CASE (LABEL(1:2))

     CASE ('-c','-C')
!       concentration units multiplier
        KLEN=INDEX(LABEL,' ')-1
        READ(LABEL(3:),'(F10.0)')CFACT

     CASE ('-d','-D')
!       write to diagnostic file
        DIAG=.TRUE.

     CASE ('-h','-H')
!       write header record
        READ(LABEL(3:),'(I10)')HEAD
        HEAD=MAX(0, MIN(HEAD,2))

     CASE ('-i','-I')
!       UNIT-10 input concentration file
        KLEN=INDEX(LABEL,' ')-1
        FNAME=LABEL(3:KLEN)
        INQUIRE (FILE=FNAME,EXIST=FTEST)
        IF(FTEST)THEN
           OPEN(10,FILE=FNAME,FORM='UNFORMATTED',ACCESS='SEQUENTIAL')
        ELSE
           WRITE(*,*)'File not found:',FNAME
           STOP
        END IF

     CASE ('-m','-M')
!       UNIT-20 measurement concentration file
        KLEN=INDEX(LABEL,' ')-1
        FNAME=LABEL(3:KLEN)
        OPEN(20,FILE=FNAME,ACTION='READ')

     CASE ('-o','-O')
!       UNIT-30 output concentration file
        KLEN=INDEX(LABEL,' ')-1
        FNAME=LABEL(3:KLEN)
        OPEN(30,FILE=FNAME)

     CASE ('-p','-P')
!       pollutant number
        KLEN=INDEX(LABEL,' ')-1
        READ(LABEL(3:),'(I10)')JPOL

     CASE ('-r','-R')
!       plume rotation from a point
        K1=INDEX(LABEL(1:),':')
        K2=INDEX(LABEL(K1+1:),':')+K1
        IF(K1.EQ.0.OR.K2.EQ.0)THEN
           ANGLE=0.0
        ELSE
           READ(LABEL(3:K1-1),'(F10.0)')ANGLE
           READ(LABEL(K1+1:K2-1),'(F10.0)')RLAT
           READ(LABEL(K2+1:),'(F10.0)')RLON
           IF(ANGLE.GT.45.0)THEN
              WRITE(*,*)'Plume rotation limited to 45 deg'
              STOP
           END IF
        END IF

     CASE ('-s','-S')
!       write supplemental info to output 
        INFO=.TRUE.

     CASE ('-x','-X')
!       concentration extraction method
        KLEN=INDEX(LABEL,' ')-1
        IF(LABEL(3:3).EQ.'i'.OR.LABEL(3:3).EQ.'I')LINEAR=1

     CASE ('-z','-Z')
!       level number
        KLEN=INDEX(LABEL,' ')-1
        READ(LABEL(3:),'(I10)')JLEV

     END SELECT
     KARG=KARG-1
  END DO

END SUBROUTINE setarg

!---------------------------------------------------------------------

SUBROUTINE SETOBS(diag,numobs)

  LOGICAL :: diag

  numobs=0
  kret=0

  READ(20,*)
  READ(20,*)

  DO WHILE (kret.EQ.0)
     READ(20,*,IOSTAT=kret)
     numobs=numobs+1
  END DO

  numobs=numobs-1
  REWIND(20)
  IF(DIAG) WRITE(40,*)'Number of observation records: ',numobs

END SUBROUTINE setobs

!-------------------------------------------------------------------------------

SUBROUTINE TMINIT

  IMPLICIT NONE

  INTEGER  :: K, IYEAR, NDPY, NDTOT
  INTEGER  :: NADPY (200)     ! accumulated days table from Jan 1, 1900        
                              ! where the first value is for Jan 1, 1901
  COMMON /TMVALS/ NADPY

  NDTOT=0
  DO IYEAR=1900,2099

     K=IYEAR-1900+1
     NDPY=365

     IF(MOD(IYEAR,100).NE.0)THEN
!       regular year
        IF(MOD(IYEAR,4).EQ.0)   NDPY=366 ! leap year
     ELSE
!       century
        IF(MOD(IYEAR,400).EQ.0) NDPY=366 ! leap century
     END IF

     NDTOT=NDTOT+NDPY
     NADPY(K)=NDTOT 

  END DO

END SUBROUTINE tminit

!-------------------------------------------------------------------------------

SUBROUTINE TM2MIN(IY,IM,ID,IH,MN,MACC)

  IMPLICIT NONE

!-------------------------------------------------------------------------------
! argument list variables
!-------------------------------------------------------------------------------

  INTEGER,  INTENT(IN)   :: iy,im,id,ih,mn       ! date and time
  INTEGER,  INTENT(OUT)  :: macc                 ! minutes accumulated 

!-------------------------------------------------------------------------------
! internal variables
!-------------------------------------------------------------------------------

  INTEGER  :: iyear,idt

! default number of accumulated days in each month (non leap year)
  INTEGER  :: NADPM (12) = (/0,31,59,90,120,151,181,212,243,273,304,334/)

! accumulated days table from Jan 1, 1900 where the first value is Jan 1, 1901
  INTEGER  :: NADPY (200)
  COMMON /TMVALS/ NADPY

  SAVE NADPM

!-------------------------------------------------------------------------------

! check for initialization
  IF(NADPY(1).NE.365)THEN
     WRITE(*,*)'*ERROR* tm2min: calendar routines not initialized (tminit)'
     STOP 900
  END IF

! compute four digit year
  IF(IY.LE.100)THEN
     IYEAR=MOD(IY+99,100)+1901
  ELSEIF(IY.LT.1941)THEN
     WRITE(*,*)'*ERROR* tm2min: year out of valid range (1941-2040) - ',IY
     STOP 900
  ELSE
     IYEAR=IY
  END IF

! assume that no data prior to 1940 hence those years represents 2000-2040
  IF(IYEAR.LT.1940)IYEAR=IYEAR+100

! number of accumulated days until this year
  IDT=NADPY(IYEAR-1900)

! add accumulated days for this year
  IDT=IDT+NADPM(IM)+(ID-1)

! adjust accumulated days for leap year
! does not account for centuries where mod(iyear,400)=0
  IF(MOD(IYEAR,4).EQ.0.AND.IM.GT.2)IDT=IDT+1

! convert to minutes
  MACC=MN+(IH+IDT*24)*60

END SUBROUTINE tm2min