!$$$  PROGRAM DOCUMENTATION BLOCK
!
! PROGRAM:     CONREAD          READS HYSPLIT GRIDDED CONCENTRATION FILE
!   PRGMMR:    ROLAND DRAXLER   ORG: R/ARL       DATE:96-06-01
!
! ABSTRACT:  THIS CODE WRITTEN AT THE AIR RESOURCES LABORATORY ...
!   SAMPLE PROGRAM TO READ THE GRIDDED HYSPLIT CONCENTRATION FILE AND
!   DUMP SELECTED STATISTICS EACH TIME PERIOD TO STANDARD OUTPUT
!
! PROGRAM HISTORY LOG:
!   LAST REVISED: 18 Feb 1997 (RRD)
!                 02 Oct 2000 (RRD) - fortran90 upgrade
!                 30 Nov 2000 (RRD) - concentration packing
!                 01 Mar 2002 (RRD) - minutes to output
!                 10 Apr 2003 (RRD) - extended diagnostics
!                 21 Nov 2003 (RRD) - added minutes to emission time
!                 22 Dec 2003 (RRD) - support for cpack=2
!                 27 Mar 2007 (BS)  - format change 
!                 28 Jul 2008 (RRD) - format change 
!
! USAGE:  CONREAD
!
!   INPUT PARAMETERS:
!     PROMPTS FOR ENTRY ON STANDARD INPUT AS REQUIRED
!   INPUT FILES:
!     UNIT 10 - binary concentration file as named on input
!   OUTPUT FILES:
!     NONE
!
! ATTRIBUTES:
!   LANGUAGE: FORTRAN 90
!   MACHINE:  IBM RS6000
!
!$$$

PROGRAM CONREAD

  IMPLICIT NONE

  CHARACTER(256)            :: LABEL
  CHARACTER(4)              :: PTYPE, MODEL

  REAL,         ALLOCATABLE :: CSUM  (:,:) ! master concentration array
  CHARACTER(4), ALLOCATABLE :: IDENT (:) 
  INTEGER,      ALLOCATABLE :: HEIGHT(:)
  REAL,         ALLOCATABLE :: LATSTN(:),LONSTN(:)

  LOGICAL    :: ftest
  INTEGER(2) :: ip,jp
  INTEGER    :: ich,ii,jj,kk,kl,kp,n,nloc,nlat,nlon,cpack,iargc 
  INTEGER    :: np,nxyp,kret,krec,ktest,ntyp,nlvl,level,narg,karg 
  INTEGER    :: iyr,imo,ida,ihr,imn,ifh,ibyr,ibmo,ibda,ibhr,ibmn
  REAL       :: cmax,cmin,clat,clon,dlat,dlon,olon,olat,olvl 
  REAL       :: latmax,lonmax,latmin,lonmin 

!======================================================================

  KARG=IARGC()
  IF(KARG.GE.1)THEN
     DO NARG=1,KARG
        CALL GETARG(NARG,LABEL)

!       help option (-H)
        IF(LABEL(1:2).EQ.'-H'.OR.LABEL(1:2).EQ.'-h')THEN
           WRITE(*,*)'SAMPLE PROGRAM TO READ THE GRIDDED HYSPLIT CONCENTRATION FILE AND'
           WRITE(*,*)'DUMP SELECTED STATISTICS EACH TIME PERIOD TO STANDARD OUTPUT'
           STOP
        END IF
     END DO
  END IF

  KTEST=0
  WRITE(*,*)'Enter concentration file name:'
  READ(*,'(A)')LABEL
  WRITE(*,*)'Extended diagnostics (0/1):',ktest
  READ(*,*)KTEST

  INQUIRE (FILE=LABEL, EXIST=FTEST)
  IF(FTEST)THEN
     OPEN(10,FILE=LABEL,FORM='UNFORMATTED',ACCESS='SEQUENTIAL')
  ELSE
     WRITE(*,*)'File not found: ',LABEL
     STOP
  END IF

! 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
  write(*,*)'Meteorology: ',model,iyr,imo,ida,ihr,ich
  write(*,*)'Packing: ',cpack

  IF(CPACK.EQ.2) ALLOCATE (latstn(1),lonstn(1))

  DO N=1,NLOC
     READ(10,IOSTAT=kret)IBYR,IBMO,IBDA,IBHR,OLAT,OLON,OLVL,IBMN
     IF(kret.ne.0)IBMN=0
     write(*,*)'Simulation start: ',ibyr,ibmo,ibda,ibhr,ibmn
     write(*,*)'Origin point    : ',olat,olon,olvl
  END DO

! horizontal grid index record
  READ(10) NLAT, NLON, DLAT, DLON, CLAT, CLON
  write(*,'(A20,2I10)')  ' Numb points lat/lon:',nlat,nlon
  write(*,'(A20,2F10.4)')'   Increment lat/lon:',dlat,dlon
  write(*,'(A20,2F10.4)')'      Corner lat/lon:',clat,clon
  ALLOCATE (csum(nlon,nlat), STAT=kret)

! vertical grid index record
  READ(10) NLVL
  BACKSPACE(10)
  ALLOCATE (height(nlvl), STAT=kret)
  READ(10) NLVL, (HEIGHT(KK),KK=1,NLVL)
  write(*,'(A20,20I6)')' Levels and heights:',nlvl,(height(kk),kk=1,nlvl)

! pollutant identification record
  READ(10) NTYP
  BACKSPACE(10)
  ALLOCATE (IDENT(ntyp), STAT=kret) 
  READ(10) NTYP, (IDENT(KK),KK=1,NTYP)
  write(*,'(A20,I5,20A5)')' Number and ID:',ntyp,(ident(kk),kk=1,ntyp)

! pause to examine header information
  write(*,*)'Enter to continue ...'
  read(*,*)

  NXYP=0
  KREC=0
  DO WHILE (KREC.GE.0)

     READ(10,END=900)IYR,IMO,IDA,IHR,IMN,IFH
     write(*,'(A15,6I5)')' Sample start:',iyr,imo,ida,ihr,imn,ifh

     READ(10)IYR,IMO,IDA,IHR,IMN,IFH
     write(*,'(A15,6I5)')' Sample stop: ',iyr,imo,ida,ihr,imn,ifh

     DO KP=1,NTYP
     DO KL=1,NLVL

        IF(CPACK.EQ.1)THEN
           CSUM = 0.0
           READ(10)PTYPE,LEVEL,NXYP,(IP,JP,CSUM(IP,JP),NP=1,NXYP)

        ELSEIF(CPACK.EQ.2)THEN
           IF(NXYP.EQ.0)THEN
!             first time through
              READ(10)PTYPE,LEVEL,NXYP,LONSTN,LATSTN,CSUM(1,1)
              IF(NXYP.GT.1)THEN
                 DEALLOCATE (LONSTN,LATSTN,CSUM)
                 ALLOCATE (LONSTN(NXYP),LATSTN(NXYP),CSUM(NXYP,1))
                 NLON=NXYP
              END IF
              BACKSPACE(10)
           END IF
           READ(10)PTYPE,LEVEL,NXYP,                         &
                  (LONSTN(II),LATSTN(II),CSUM(II,1),II=1,NXYP)
             
        ELSEIF(CPACK.EQ.0)THEN
           READ(10)PTYPE,LEVEL,((CSUM(II,JJ),II=1,NLON),JJ=1,NLAT)

        ELSE
           WRITE(*,*)'Unsupported concentration packing, CPACK=',CPACK
           STOP
        END IF

!       find max/min value and location
        CMAX=0.0
        CMIN=1.0E+25
        IF(KTEST.EQ.1)WRITE(*,*)
        DO JJ=1,NLAT
        DO II=1,NLON
           IF(CSUM(II,JJ).GT.CMAX)THEN
              CMAX=CSUM(II,JJ)
              IF(CPACK.NE.2)THEN
                 lonmax=(II-1)*DLON+CLON
                 latmax=(JJ-1)*DLAT+CLAT
              ELSE
                 lonmax=lonstn(ii)
                 latmax=latstn(ii)
              END IF
           END IF
           IF(CSUM(II,JJ).GT.0.0.AND.CSUM(II,JJ).LT.CMIN)THEN
              CMIN=CSUM(II,JJ)
              IF(CPACK.NE.2)THEN
                 lonmin=(II-1)*DLON+CLON
                 latmin=(JJ-1)*DLAT+CLAT
              ELSE
                 lonmin=lonstn(ii)
                 latmin=latstn(ii)
              END IF
           END IF
           IF(KTEST.EQ.1.AND.CSUM(II,JJ).GT.0.0)WRITE(*,*)II,JJ,CSUM(II,JJ)
        END DO
        END DO
        IF(KTEST.EQ.1)WRITE(*,*)

        if(cmax.eq.0.0)then
           write(*,'(12X,3A,I5,5X,A)')                     &
          'POLID=',ptype,' LEVEL=',level,'  CON= 0.0'  
        else
           write(*,'(12X,3A,I5,5X,A,E10.3,A,2F9.3)')      &
          'POLID=',ptype,' LEVEL=',level,'  MAX=',cmax,'   AT:',latmax,lonmax
           write(*,'(12X,3A,I5,5X,A,E10.3,A,2F9.3)')      &
          'POLID=',ptype,' LEVEL=',level,'  MIN=',cmin,'   AT:',latmin,lonmin
        end if

     END DO
     END DO
     WRITE(*,'(A)')' '
  END DO

900 DEALLOCATE (csum,ident,height,STAT=kret)

END PROGRAM conread
