!--------------------------------------------------------------------
! Program to read merged DATEM format measured and calculated air 
! concentration files and convert to grads format for plotting 
!--------------------------------------------------------------------
! Unit 10 - measurement and model results merged data file
! Unit 20 - grads formatted data file
! Unit 30 - grads control file
! Unit 40 - sample grads display script
! Unit 50 - dummy grid definition file
!--------------------------------------------------------------------
! Last Revised: 10 Aug 2001
!--------------------------------------------------------------------

PROGRAM C2GRADS

  IMPLICIT NONE


  CHARACTER(8), ALLOCATABLE :: SITE(:)               ! station label
  INTEGER,      ALLOCATABLE :: IYR(:),IMO(:),IDA(:)  ! sample date
  INTEGER,      ALLOCATABLE :: IHR(:),IDUR(:)        ! time & duration 
  REAL,         ALLOCATABLE :: SLAT(:),SLON(:)       ! location
  REAL,         ALLOCATABLE :: CALC(:),MEAS(:)       ! concentrations

  INTEGER,      ALLOCATABLE :: INDX(:),TVAL(:)

  REAL          :: RTIME,DUMY,LOGM,LOGC,RCOM,BIAS
  CHARACTER(80) :: LABEL
  INTEGER       :: IBYR,IBMO,IBDA,IBHR
  INTEGER       :: K,M,N,KRET,NLEV,NREC,NUMOBS,NFLAG,DELTA

  character months(12)*3
  data months/'jan','feb','mar','apr','may','jun',      &
              'jul','aug','sep','oct','nov','dec'/

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

  INTERFACE

    SUBROUTINE SETARG (LABEL)
    IMPLICIT NONE
    CHARACTER(80), INTENT(OUT) :: LABEL
    END SUBROUTINE SETARG

    SUBROUTINE SORT(NUMPTS,VALUE,INDX)
    IMPLICIT NONE
    INTEGER,  INTENT(IN)  :: NUMPTS
    INTEGER,  INTENT(IN)  :: VALUE(:)
    INTEGER,  INTENT(OUT) :: INDX (:)
    END SUBROUTINE SORT

  END INTERFACE

!--------------------------------------------------------------------
! read any command line arguments and set default values
  CALL SETARG (LABEL)

! open data file to determine array allocations
  READ(10,'(a)')
  DO WHILE (kret.EQ.0)
     READ(10,'(a)',IOSTAT=kret)
     numobs=numobs+1
  END DO
  numobs=numobs-1
  REWIND(10)

  ALLOCATE (site (numobs), STAT=kret)
  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 (calc (numobs), STAT=kret)
  ALLOCATE (meas (numobs), STAT=kret)

  ALLOCATE (indx (numobs), STAT=kret)
  ALLOCATE (tval (numobs), STAT=kret)

!--------------------------------------------------------------------
! load measurement and model data and sort by time

  READ(10,'(a)')
  DO N=1,NUMOBS
     READ(10,'(I5,A,F7.2,F8.2,4I5,2F10.2)')         &
          K,SITE(N),SLAT(N),SLON(N),                &
	      IYR(N),IMO(N),IDA(N),IHR(N),MEAS(N),CALC(N)
!    compute psuedo time value for temporal sorting
     TVAL(N)=(IYR(N)*366+IMO(N)*31+IDA(N))*24+IHR(N)
  END DO
  CLOSE(10)

  CALL SORT(NUMOBS,TVAL,INDX)

!--------------------------------------------------------------------
! output grads data file 

! OPEN(20,FILE='merged.dat',ACCESS='direct',FORM='unformatted',RECL=28)
  OPEN(20,FILE='merged.dat',FORM='unformatted')

  M=0
  NREC=0
  NFLAG=1
  DELTA=0
  RTIME=0.0

  DO N=1,NUMOBS
     K=INDX(N)

!    start of new time group
     IF(N.GT.1)THEN
  	    IF(TVAL(K).GT.TVAL(INDX(N-1)))THEN
           M=M+1
		   NLEV=0
!          WRITE(20,REC=M)SITE(K),SLAT(K),SLON(K),RTIME,NLEV,NFLAG
           WRITE(20      )SITE(K),SLAT(K),SLON(K),RTIME,NLEV,NFLAG

           NREC=NREC+1
!          set time interval
           IF(DELTA.EQ.0)DELTA=TVAL(K)-TVAL(INDX(N-1))

        END IF
     ELSE
        IBYR=IYR(K)
        IBMO=IMO(K)
        IBDA=IDA(K)
        IBHR=IHR(K)/100
     END IF

!    header of length 28 bytes
     NLEV=1
     M=M+1
!    WRITE(20,REC=M)SITE(K),SLAT(K),SLON(K),RTIME,NLEV,NFLAG
     WRITE(20      )SITE(K),SLAT(K),SLON(K),RTIME,NLEV,NFLAG

!    process data
     DUMY=-999.0
     LOGM=-999.0
     LOGC=-999.0
     RCOM=-999.0
  
     BIAS=CALC(K)-MEAS(K)
     IF(MEAS(K).GT.0.0)RCOM=CALC(K)/MEAS(K)
     IF(CALC(K).GT.0.0)LOGC=ALOG10(CALC(K))
     IF(MEAS(K).GT.0.0)LOGM=ALOG10(MEAS(K))

!    surface data 
     M=M+1
!    WRITE(20,REC=M)MEAS(K),CALC(K),LOGM,LOGC,RCOM,BIAS,DUMY
     WRITE(20      )MEAS(K),CALC(K),LOGM,LOGC,RCOM,BIAS,DUMY
   END DO

!  end of file
   NLEV=0
   M=M+1
   NREC=NREC+1
!  WRITE(20,REC=M)'        ',0.0,0.0,0.0,NLEV,NFLAG
   WRITE(20      )'        ',0.0,0.0,0.0,NLEV,NFLAG
   CLOSE(20)

!--------------------------------------------------------------------
! output grads control file
 
  open(30,file='merged.ctl')
  write(30,'(a)')'dset merged.dat'
  write(30,'(a)')'dtype station'
  write(30,'(a)')'stnmap merged.map'
  k=index(label,' ')-1
  write(30,'(a)')'title '//label(1:k)
  write(30,'(a)')'undef -999.0'
  write(30,'(a)')'options sequential'
  write(30,'(a,i4,a,i2.2,a1,i2.2,a3,2i4,a)') 'tdef ',nrec,   & 
           ' linear ',ibhr,'z',ibda,months(ibmo),ibyr,delta,'hr'
  write(30,'(a)')'vars 7'
  write(30,'(a)')'meas 0 99 Measured Concentration'
  write(30,'(a)')'calc 0 99 Calculated Concentration'
  write(30,'(a)')'logm 0 99 LOG10 Measured'
  write(30,'(a)')'logc 0 99 LOG10 Calculated'
  write(30,'(a)')'rcom 0 99 Ratio Calculated to Measured'
  write(30,'(a)')'bias 0 99 Bias Calculated-Measured'
  write(30,'(a)')'dumy 0 99 Dummy Field'
  write(30,'(a)')'endvars'
  close(30)

!--------------------------------------------------------------------
! sample grads display script

  open(40,file='merged.gs')
  write(40,'(a)' )'''set gxout shaded'''
  write(40,'(a)' )'''set clevs 10.0 20.0 50.0 100.0 200.0'''
  write(40,'(a)' )'''set ccols 5 10 8 7 2'''
  write(40,'(a)' )'''d oacres(dummy,meas.2)'''
  write(40,'(a)' )'''set gxout contour'''
  write(40,'(a)' )'''set cthick 8'''
  write(40,'(a)' )'''set clevs 10.0 20.0 50.0 100.0 200.0'''
  write(40,'(a)' )'''ccols 1 1 1 1 1'''
  write(40,'(a)' )'''d oacres(dummy,calc.2)'''
  write(40,'(3a)')'''draw title ',label(1:k),''''
  close(40)

!--------------------------------------------------------------------
! dummy grid definition control file
 
  open(50,file='dummy.ctl')
  write(50,'(a)')'dset dummy.dat'
  write(50,'(a)')'options template'
  write(50,'(a)')'undef -999.0'
  write(50,'(a)')'title dummy control file'
  write(50,'(a)')'xdef  75 linear  -120.0 1.0'
  write(50,'(a)')'ydef  50 linear    20.0 1.0'
  write(50,'(a)')'zdef   1 levels  0'
  write(50,'(a,i2.2,a1,i2.2,a3,2i4,a)') 'tdef 1 linear ',        & 
                        ibhr,'z',ibda,months(ibmo),ibyr,delta,'hr'
  write(50,'(a)')'vars   1'
  write(50,'(a)')'dummy  0 99 dummy variable'
  write(50,'(a)')'endvars'
  close(50)

END PROGRAM C2GRADS

!--------------------------------------------------------------------
! Process command line arguments and open all files for either input
! or output.
!--------------------------------------------------------------------
! Last Revised: 09 Aug 2001
!--------------------------------------------------------------------

SUBROUTINE SETARG (LABEL)

  IMPLICIT NONE

  CHARACTER(80), INTENT(OUT) :: LABEL

  CHARACTER(80) :: DNAME,MODEL
  CHARACTER(1)  :: TRACER
  LOGICAL       :: FTEST
  INTEGER       :: KLEN,KLEN1,KLEN2,KARG,IARGC

  KARG=IARGC()
  IF(KARG.LT.1)THEN
     WRITE(*,*)'USAGE: c2grads -d -t -m'
     WRITE(*,*)'-d[directory name of experiment]'
     WRITE(*,*)'-t[tracer number: 1,2,3, etc]'
     WRITE(*,*)'-m[model variation string]'
     READ(*,*)
     STOP
  END IF

! default values if not set on command line
  TRACER='1'
  MODEL='_001'
  DNAME='experiment'
 
! go through each argument
  DO WHILE (KARG.GT.0)

     CALL GETARG(KARG,LABEL)

     SELECT CASE (LABEL(1:2))

!    experiment directory
     CASE ('-d','-D')
        KLEN1=INDEX(LABEL,' ')-1
        DNAME=LABEL(3:KLEN1)
        KLEN1=INDEX(DNAME,' ')-1

!    model variation
     CASE ('-m','-M')
        KLEN2=INDEX(LABEL,' ')-1
        MODEL='_'//LABEL(3:KLEN2)
        KLEN2=INDEX(MODEL,' ')-1

!    tracer number
     CASE ('-t','-T')
        TRACER=LABEL(3:3)

     END SELECT
     KARG=KARG-1
  END DO

! open the merged data file
  LABEL='../mdl_stat/'//DNAME(1:KLEN1)//'/dat'//TRACER//  &
        MODEL(1:KLEN2)//'.txt'
  INQUIRE (FILE=LABEL,EXIST=FTEST)
  IF(FTEST)THEN
     OPEN(10,FILE=LABEL,ACTION='READ')
  ELSE
     WRITE(*,*)'File not found:',LABEL
     READ(*,*)
     STOP
  END IF

! return label string for identification
  LABEL=DNAME(1:KLEN1)//'/dat'//TRACER//MODEL(1:KLEN2)//'.txt'

END SUBROUTINE setarg

!-------------------------------------------------------------------- 
! Ascending order binary sort which returns the sorted values in
! array DSORT by cumulative percent from 1 to 99.  For instance, 
! DSORT(50) returns the median value, while DSORT(25 & 75) return the
! quartile values. Routine requres 2**N passes through the data, 
! where N equals the number of data points.
!--------------------------------------------------------------------
! Last Revised: 18 Jul 2001
!--------------------------------------------------------------------

SUBROUTINE SORT(NUMPTS,VALUE,INDX)

  IMPLICIT NONE

  INTEGER,  INTENT(IN)  :: NUMPTS
  INTEGER,  INTENT(IN)  :: VALUE(:)
  INTEGER,  INTENT(OUT) :: INDX (:)

  INTEGER, ALLOCATABLE :: TEMP (:)

  INTEGER :: I,I1,I2,I3,L1,L2,KRET,NLEN,LEN1,LEN2,INDX1,INDX2

  IF(NUMPTS.LE.1)THEN
     WRITE(*,*)'Too few points to sort: ',NUMPTS
	 RETURN
  END IF

  ALLOCATE (TEMP(NUMPTS), STAT=KRET)

! initialize index array
  DO I=1,NUMPTS
     INDX(I)=I
  END DO

    NLEN=1
100 IF(NLEN.GE.NUMPTS)GO TO 800

    I1=1
    I2=I1+NLEN
    I3=1

200 L1=0
    L2=0
    LEN1=NUMPTS-I1+1
    IF(LEN1.GT.NLEN)LEN1=NLEN
    LEN2=NUMPTS-I2+1
    IF(LEN2.GT.NLEN)LEN2=NLEN
250 INDX1=INDX(I1+L1)
    INDX2=INDX(I2+L2)

    IF(VALUE(INDX1).GT.VALUE(INDX2))GO TO 400

300 TEMP(I3)=INDX(I1+L1)
    L1=L1+1
    I3=I3+1
    GO TO 500

400 CONTINUE
    TEMP(I3)=INDX(I2+L2)
    L2=L2+1
    I3=I3+1

500 IF(I3.GT.NUMPTS)GO TO 550
    IF(L1.LT.LEN1.AND.L2.LT.LEN2)GO TO 250
    IF(L1.EQ.LEN1.AND.L2.LT.LEN2)GO TO 400
    IF(L1.LT.LEN1.AND.L2.EQ.LEN2)GO TO 300
    I1=I3
    I2=I1+NLEN
    IF(I2.LE.NUMPTS)GO TO 200
    L1=0
    GO TO 300

550 DO I=1,NUMPTS
       INDX(I)=TEMP(I)
    END DO
    NLEN=NLEN*2
    GO TO 100

800 CONTINUE

DEALLOCATE (temp)

END SUBROUTINE SORT
