!*******************************************************************************
!$RCSfile: Class1answers.f90,v $
!$Revision: 1.0 $
!$Date: 2013/11/15 19:03:53 $
!*******************************************************************************
MODULE class1Tools

  USE tooluser_fd
  USE Extract_fi
  USE SCIPtool
  USE GetTimes_fi
  USE scipuff_fi, ONLY: create
  USE MPI_fi, ONLY: useMPI,myid
  IMPLICIT NONE

  TYPE an_output
     character (len=PATH_MAXLENGTH) :: outFile    ! output filename
     character (len=3)              :: ConDepVis  ! output type: con, dep, vis
     character (len=128)            :: chemSpecies! chemical species name
     integer                        :: avgPeriods ! number of hours to average
     integer                        :: rank       ! 1st or 4th or 8th high...
     integer                        :: lun        ! unit to write to
     real                           :: exceeds    ! print all vals over this
  end TYPE an_output

  REAL, DIMENSION(:,:,:),  ALLOCATABLE :: con       !index: iArea,iRec,ikind
  REAL, DIMENSION(:,:,:),  ALLOCATABLE :: dep       !index: iArea,iRec,ikind
  REAL, DIMENSION(:,:,:,:),ALLOCATABLE :: srfFldVal !index: iArea,iRec,ikind,2
  REAL, DIMENSION(:,:),    ALLOCATABLE :: lat,lon   !index: iArea,iRec
  REAL, DIMENSION(:,:),    ALLOCATABLE :: xRec,yRec !index: iArea,iRec
  REAL, DIMENSION(:,:,:,:),ALLOCATABLE :: max_val   !index: iArea,iRec,iOut,iRank
  REAL, DIMENSION(:,:,:),  ALLOCATABLE :: sum_val   !index: iArea,iRec,iOut
  REAL, DIMENSION(:,:),    ALLOCATABLE :: Bext_nat  !index: iArea, Month
  REAL, DIMENSION(:,:),    ALLOCATABLE :: NatCond, LargeRH, SmallRH, SeaSalt

  INTEGER,DIMENSION(:,:,:,:),ALLOCATABLE :: max_time !index: same as max_val
  INTEGER,DIMENSION(:,:),  ALLOCATABLE :: iHours    ! index: iArea,iOut
  INTEGER,DIMENSION(:),    ALLOCATABLE :: nRec      ! index: iArea

! The iKind index of the species we need to calc max concentrations

  INTEGER, DIMENSION(:),   ALLOCATABLE :: iCon

! The iKind index of the species we need to calc total S and N deposition

  INTEGER :: iNO, iNO2, iNO3, iN2O5, iHNO3, iHONO, iPNA, iPAN, iPANX, iNTR
  INTEGER :: iPNO3, iANO3K, iPNH4, iANH4K, iNH3
  INTEGER :: iSO2,  iSULF,  iPSO4, iASO4K

! The iKind index of the species we need to calculate visibility

  ! Already declared iPSO4, iASO4K, iPNO3, iANO3K above
  INTEGER :: iPSOA, iPOC, iPEC, iASOIL, iACORS

  CHARACTER (len=PATH_MAXLENGTH), allocatable, dimension(:) :: class1Filename
  CHARACTER (len=20),             allocatable, dimension(:) :: class1Name

  CHARACTER (len=PATH_MAXLENGTH)    :: flagDir, class1AreaDir
  CHARACTER (len=PATH_MAXLENGTH)    :: projName, controlFile
  CHARACTER (len=PATH_MAXLENGTH)    :: NatCondFile, LargeRHFile
  CHARACTER (len=PATH_MAXLENGTH)    :: SmallRHFile, SeaSaltFile

  CHARACTER (len=1),   dimension(3) :: comment   = (/ "#", ";", "!" /)
  CHARACTER (len=1),   dimension(5) :: delim     = (/ " ", ",", "#", ";", "!" /)
  CHARACTER (len=1),   dimension(3) :: datedelim = (/ "-", "_", ":" /)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
CONTAINS
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  SUBROUTINE Class1con( iArea, iKind, iOut, output, iymdh)

    IMPLICIT NONE

    TYPE(an_output)         :: output ! what to calculate
    integer                 :: iArea  ! which Class I area to process
    integer                 :: iKind  ! which chemical species in con() to use
    integer                 :: iOut   ! which output we're working on
    integer                 :: ir     ! rank of this average (H1H, H8H, etc.)
    integer                 :: iRec   ! Class I area receptor index
    integer                 :: iymdh  ! current time in YYYYMMDDHH format

!--- Find the n-periond average value, then find the rank of this value,
!    and insert the value at the right rank

!    print*,"iArea,iOut,iHours,sum_val= ",iArea,iOut,iHours(iArea,iOut),sum_val(iArea,10,iOut)

    if (iHours(iArea,iOut) == output%avgPeriods) then ! end averaging period
!       print*
!       print*,"averging for output ",iArea,iOut
!       print*,"sum_val = ",sum_val(iArea,10,iOut)
!       print*,"avg_val = ",sum_val(iArea,10,iOut)/iHours(iArea,iOut)
!       print*

       do iRec = 1, nRec(iArea)

          sum_val(iArea,iRec,iOut) = sum_val(iArea,iRec,iOut)/iHours(iArea,iOut)

          ir = find_rank( sum_val(iArea,iRec,iOut), iArea, iRec, iOut, &
               output%rank, .false. )

          call insert( sum_val(iArea,iRec,iOut), iArea, iRec, iOut,   &
               output%rank, ir, iymdh )  ! inserts into max_val at right rank

          if ( sum_val(iArea,iRec,iOut) > output%exceeds )             &
               write(output%lun,11) trim(class1Name(iArea)),           &
               trim(output%chemSpecies),                               &
               lat(iArea,iRec),  lon(iArea,iRec),                      &
!               xRec(iArea,iRec), yRec(iArea,iRec),                     &
               iymdh, output%avgPeriods,">",output%exceeds,            &
               sum_val(iArea,iRec,iOut)
11        format(2(a5,","),2(f12.6,","),i11.10,",",i6,",",a,g16.6,",",g16.6)

       end do

       iHours(iArea,iOut) = 0
       sum_val(iArea,:,iOut) = 0.

    end if ! if (iHours(iArea,iOut) == output%avgPeriods) then

!--- add this hour to the running average

    do iRec = 1, nRec(iArea)
       sum_val(iArea,iRec,iOut) =                                       &
            sum_val(iArea,iRec,iOut) + con(iArea,iRec,iKind)
    end do

    iHours(iArea,iOut) = iHours(iArea,iOut) + 1

    RETURN

  END SUBROUTINE Class1con
!
!*************************************************************************
!
  SUBROUTINE Class1dep( iArea, iOut, output, iymdh )

    IMPLICIT NONE

    TYPE(an_output)         :: output ! what to calculate
    integer                 :: iArea  ! which Class I area to process
    integer                 :: iOut   ! which output we're working on
    integer                 :: ir     ! rank of this average (H1H, H8H, etc.)
    integer                 :: iRec   ! Class I area receptor index
    integer                 :: iymdh  ! current time in YYYYMMDDHH format

!--- Find the n-periond average value, then find the rank of this value,
!    and insert the value at the right rank

    if (iHours(iArea,iOut) == output%avgPeriods) then ! end averaging period

       do iRec = 1,nRec(iArea)

          sum_val(iArea,iRec,iOut) = sum_val(iArea,iRec,iOut)/iHours(iArea,iOut)

          ir = find_rank( sum_val(iArea,iRec,iOut), iArea, iRec, iOut, &
               output%rank, .false. )

          call insert( sum_val(iArea,iRec,iOut), iArea, iRec, iOut,   &
               output%rank, ir, iymdh )    ! inserts into max_val at right rank

          if ( sum_val(iArea,iRec,iOut) > output%exceeds )             &
               write(output%lun,11) trim(class1Name(iArea)),           &
               trim(output%chemSpecies),                               &
               lat(iArea,iRec),  lon(iArea,iRec),                      &
!               xRec(iArea,iRec), yRec(iArea,iRec),                     &
               iymdh, output%avgPeriods,">",output%exceeds,            &
               sum_val(iArea,iRec,iOut)
11        format(2(a5,","),2(f12.6,","),i11.10,",",i6,",",a,g16.6,",",g16.6)

       end do

       iHours(iArea,iOut) = 0
       sum_val(iArea,:,iOut) = 0.

    end if ! if (iHours(iArea,iOut) == output%avgPeriods) then

!--- add this hour to the running total deposition of N and S

    do iRec = 1, nRec(iArea)

       if (output%chemSpecies == "N") then ! Total Nitrogen Deposition

          if (iNO   > 0) sum_val(iArea,iRec,iOut) = sum_val(iArea,iRec,iOut) &
               + 14./30.  * dep(iArea,iRec,iNO   ) ! NO
          if (iNO2  > 0) sum_val(iArea,iRec,iOut) = sum_val(iArea,iRec,iOut) &
               + 14./46.  * dep(iArea,iRec,iNO2  ) ! NO2
          if (iNO3   > 0) sum_val(iArea,iRec,iOut) = sum_val(iArea,iRec,iOut) &
               + 14./62.  * dep(iArea,iRec,iNO3  ) ! NO3
          if (iN2O5  > 0) sum_val(iArea,iRec,iOut) = sum_val(iArea,iRec,iOut) &
               + 28./102. * dep(iArea,iRec,iN2O5 ) ! N2O5
          if (iHNO3  > 0) sum_val(iArea,iRec,iOut) = sum_val(iArea,iRec,iOut) &
               + 14./63.  * dep(iArea,iRec,iHNO3 ) ! HNO3
          if (iHONO  > 0) sum_val(iArea,iRec,iOut) = sum_val(iArea,iRec,iOut) &
               + 14./47.  * dep(iArea,iRec,iHONO ) ! HONO
          if (iPNA   > 0) sum_val(iArea,iRec,iOut) = sum_val(iArea,iRec,iOut) &
               + 14./76.  * dep(iArea,iRec,iPNA  ) ! HNO4
          if (iPAN   > 0) sum_val(iArea,iRec,iOut) = sum_val(iArea,iRec,iOut) &
               + 14./121. * dep(iArea,iRec,iPAN  ) ! CH3COOONO2
          if (iPANX  > 0) sum_val(iArea,iRec,iOut) = sum_val(iArea,iRec,iOut) &
               + 14./121. * dep(iArea,iRec,iPANX ) ! Higher PAN
          if (iNTR   > 0) sum_val(iArea,iRec,iOut) = sum_val(iArea,iRec,iOut) &
               + 14./130. * dep(iArea,iRec,iNTR  ) ! Org. Nitr
          if (iPNO3  > 0) sum_val(iArea,iRec,iOut) = sum_val(iArea,iRec,iOut) &
               + 28./96.  * dep(iArea,iRec,iPNO3 ) ! NH4NO3
          if (iANO3K > 0) sum_val(iArea,iRec,iOut) = sum_val(iArea,iRec,iOut) &
               + 28./96.  * dep(iArea,iRec,iANO3K) ! NH4NO3
          if (iNH3   > 0) sum_val(iArea,iRec,iOut) = sum_val(iArea,iRec,iOut) &
               + 14./17.  * dep(iArea,iRec,iNH3  ) ! ammonia
          if (iPNH4  > 0) sum_val(iArea,iRec,iOut) = sum_val(iArea,iRec,iOut) &
               + 14./18.  * dep(iArea,iRec,iPNH4 ) ! PM Ammonia
          if (iANH4K > 0) sum_val(iArea,iRec,iOut) = sum_val(iArea,iRec,iOut) &
               + 14./18.  * dep(iArea,iRec,iANH4K) ! PM Ammonia

       elseif (output%chemSpecies == "S") then ! Total Sulfur Dep

          if (iSO2   > 0) sum_val(iArea,iRec,iOut) = sum_val(iArea,iRec,iOut) &
               + 32./64.  * dep(iArea,iRec,iSO2  )   ! SO2
          if (iSULF  > 0) sum_val(iArea,iRec,iOut) = sum_val(iArea,iRec,iOut) &
               + 32./98.  * dep(iArea,iRec,iSULF )   ! H2SO4
          if (iPSO4  > 0) sum_val(iArea,iRec,iOut) = sum_val(iArea,iRec,iOut) &
               + 32./96.  * dep(iArea,iRec,iPSO4 )   ! (NH4)2SO4
          if (iASO4K > 0) sum_val(iArea,iRec,iOut) = sum_val(iArea,iRec,iOut) &
               + 32./96.  * dep(iArea,iRec,iASO4K)   ! (NH4)2SO4

       end if

    end do  ! do iRec = 1, nRec(iArea)

    iHours(iArea,iOut) = iHours(iArea,iOut) + 1

    RETURN

  END SUBROUTINE Class1dep
!
!*************************************************************************
!
  SUBROUTINE Class1vis( iArea, iOut, mo, output, iymdh )

    IMPLICIT NONE

    TYPE(an_output)         :: output ! what to calculate
    integer                 :: iArea  ! which Class I area to process
    integer                 :: iOut   ! which output we're working on
    integer                 :: mo     ! month of current time-step
    integer                 :: ir     ! rank of this average (H1H, H8H, etc.)
    integer                 :: iRec   ! Class I area receptor index
    integer                 :: iymdh  ! current time in YYYYMMDDHH format

    real                    :: TotSulfate, LgSulfate, SmSulfate
    real                    :: TotNitrate, LgNitrate, SmNitrate
    real                    :: TotOrganic, LgOrganic, SmOrganic
    real                    :: Bext, dBext
    real                    :: con_PEC, con_ASOIL, con_ACORS, con_NO2

!--- Find the n-periond average value, then find the rank of this value,
!    and insert the value at the right rank

    if (iHours(iArea,iOut) == output%avgPeriods) then ! end averaging period

       do iRec = 1, nRec(iArea) ! sum_val is dBext here

          sum_val(iArea,iRec,iOut) = sum_val(iArea,iRec,iOut)/iHours(iArea,iOut)

          ir = find_rank( sum_val(iArea,iRec,iOut), iArea, iRec, iOut, &
               output%rank, .false. )

          call insert( sum_val(iArea,iRec,iOut), iArea, iRec, iOut,   &
               output%rank, ir, iymdh )  ! inserts into max_val at right rank

          if ( sum_val(iArea,iRec,iOut) > output%exceeds )             &
               write(output%lun,11) trim(class1Name(iArea)),           &
               trim(output%chemSpecies),                               &
               lat(iArea,iRec),  lon(iArea,iRec),                      &
!               xRec(iArea,iRec), yRec(iArea,iRec),                     &
               iymdh, output%avgPeriods,">",output%exceeds,            &
               sum_val(iArea,iRec,iOut)
11        format(2(a5,","),2(f12.6,","),i11.10,",",i6,",",a,g16.6,",",g16.6)

       end do

       iHours(iArea,iOut) = 0
       sum_val(iArea,:,iOut) = 0.

    end if ! if (iHours(iArea,iOut) == output%avgPeriods) then

!--- Calculate this hour's visibility obscuration

    do iRec = 1, nRec(iArea)

!--- Calculate Total Sulfate, and Large Sulfate fraction, according to
!    Figure 5 from FLAG (2010).  Ditto for Nitate and Organic Mass.
!    If some component of these are not found in the DOS file, then set
!    the concentrations to zero, but continue with the calculation.

       TotSulfate = NatCond(iArea,1)
       if (iPSO4  > 0) TotSulfate = TotSulfate + con(iArea,iRec,iPSO4)
       if (iASO4K > 0) TotSulfate = TotSulfate + con(iArea,iRec,iASO4K)
       if (TotSulfate < 20.) then              ! 20 ug/m3
          LgSulfate = TotSulfate / 20. * TotSulfate
       else
          LgSulfate = TotSulfate
       end if
       SmSulfate = TotSulfate - LgSulfate

       TotNitrate = NatCond(iArea,2)
       if (iPNO3  > 0) TotSulfate = TotSulfate + con(iArea,iRec,iPNO3)
       if (iANO3K > 0) TotSulfate = TotSulfate + con(iArea,iRec,iANO3K)
       if (TotNitrate < 20.) then              ! 20 ug/m3
          LgNitrate = TotNitrate / 20. * TotNitrate
       else
          LgNitrate = TotNitrate
       end if
       SmNitrate = TotNitrate - LgNitrate

       TotOrganic = NatCond(iArea,3)
       if (iPSOA  > 0) TotSulfate = TotSulfate + con(iArea,iRec,iPSOA)
       if (iPOC   > 0) TotSulfate = TotSulfate + con(iArea,iRec,iPOC)
       if (TotOrganic < 20.) then              ! 20 ug/m3
          LgOrganic = TotOrganic / 20. * TotOrganic
       else
          LgOrganic = TotOrganic
       end if
       SmOrganic = TotOrganic - LgOrganic

       con_PEC   = 0.  ! If these values not found in this run, then set
       con_ASOIL = 0.  ! them to zero, but continue with the calculation.
       con_ACORS = 0.
       con_NO2   = 0.
       if (iPEC > 0)   con_PEC   = con(iArea,iRec,iPEC)
       if (iASOIL > 0) con_ASOIL = con(iArea,iRec,iASOIL)
       if (iACORS > 0) con_ACORS = con(iArea,iRec,iACORS)
       if (iNO2   > 0) con_NO2   = con(iArea,iRec,iNO2)

!--- Now calculate the total (project + natural background) extinction, and
!    the "delta Bext" percent change to natural background.  Note that
!    Sm/Lg Sulfate/Nitrate/Organic already have the Project Component added
!    above.

       Bext =                                      &
            2.2 * SmallRH(iArea,mo) * SmSulfate  + &       ! small sulfate
            4.8 * LargeRH(iArea,mo) * LgSulfate  + &       ! large sulfate
            2.4 * SmallRH(iArea,mo) * SmNitrate  + &       ! small nitrate
            5.1 * LargeRH(iArea,mo) * LgNitrate  + &       ! large nitrate
            2.8 *                     SmOrganic  + &       ! small organic mass
            6.1 *                     LgOrganic  + &       ! large organic mass
            10. * (NatCond(iArea,4) + con_PEC)   + &       ! elemental carbon
            1.0 * (NatCond(iArea,5) + con_ASOIL) + &       ! fine soil
            0.6 * (NatCond(iArea,6) + con_ACORS) + &       ! coarse mass
            1.7 * SeaSalt(iArea,mo) * NatCond(iArea,7) + & ! sea salt
            NatCond(iArea,8)                     + &       ! Rayleigh
            0.1755 * con_NO2                               ! NO2

       dBext = (Bext - Bext_nat(iArea,mo)) / Bext_nat(iArea,mo) * 100. ! in %

!--- Finally, save the change (delta) in extinction for this hour, as a
!    running sum.

       sum_val(iArea,iRec,iOut) = sum_val(iArea,iRec,iOut) + dBext

    end do   ! do iRec = 1,nRec

    iHours(iArea,iOut) = iHours(iArea,iOut) + 1

    RETURN

  END SUBROUTINE Class1vis
!
!*************************************************************************
!
  SUBROUTINE read_control(nAreas, nOuts, output)

    implicit none

    integer                                 :: nAreas, nOuts
    integer                                 :: ios, i
    integer                                 :: num_words
    integer, parameter                      :: lun = 1
    TYPE(an_output), dimension(nOuts)       :: output
    character (len=99)                      :: line, word, word2

    nAreas  = 0 ! initialize
    nOuts   = 0

99  format(a99)

    write(*,*)
    write(*,'(3a,$)') " Parsing ",trim(controlfile),", "
    call flush(6)
    open(lun,file=controlFile,status='old',err=999)

    read(lun,99,iostat=ios) line
    call count_words(line,num_words)    ! read to first non-comment non-blank
    do while ((any(line(1:1) == comment) .or. num_words == 0) .and. ios == 0)
       read(lun,99,iostat=ios) line
       call count_words(line,num_words)
    end do

    do while (ios == 0)                 ! loop over lines in control file

       call get_word(line,1,word)
       call CUPPER( word )

       if ( word(1:7) == "PROJECT" ) call get_word(line,2,projName)

       if ( word(1:5) == "AREA" ) then
          nAreas = nAreas + 1
          call get_word( line,2,class1Filename(nAreas) )
          call ADDPATH( class1Filename(nAreas),class1AreaDir)
       end if

       if (word(1:3) == "CON" .or. word(1:3) == "DEP" .or.               &
           word(1:3) == "VIS") then
          nOuts = nOuts + 1

          output(nOuts)%ConDepVis = word(1:3)

          call get_word( line, 2, output(nOuts)%chemSpecies ) ! 2nd word = species
          call CUPPER(output(nOuts)%chemSpecies)         ! must be uppercase

          call get_word( line, 3, word2 )                ! 3rd word is period
          read(word2,*) output(nOuts)%avgPeriods
          call get_word( line, 4, word2 )                ! 4th word is hr,mo,yr
          if (word2 == 'yr') output(nOuts)%avgPeriods = &
               output(nOuts)%avgPeriods * 8760            ! number of hours per yr

          call get_word( line, 5, word2 )                ! 5th word is H1H or >5%
          call CUPPER( word2 )
          if (word2(1:1) == "H") then             ! a word like H1H or H8H or H2H
             i = 2
             do while (word2(i:i) /= "H" .and. i < len_trim(word2))
                i = i + 1
             end do
             read(word2(2:i-1),*) output(nOuts)%rank
          else if (word2(1:1) == ">") then
             i = 2
             do while (word2(i:i) /= " " .and. i < len_trim(word2))
                i = i + 1
             end do
             read(word2(2:i),*) output(nOuts)%exceeds
          else
             write(*,*) "Error parsing 5th word of the following line:"
             write(*,*) trim(line)
             write(*,*) "5th word must start with 'H' or '>'."
             stop
          end if

          call get_word( line, 6, output(nOuts)%outFile )      ! 6th word is filename
       end if

       read(lun,99,iostat=ios) line
       call count_words(line,num_words)     ! non-blank line
       do while (num_words == 0 .and. ios == 0)
          read(lun,99,iostat=ios) line
          call count_words(line,num_words)
       end do

    end do ! do while (ios == 0) ! loop over lines in control file

    write(*,'(i3,a)') nOuts," outputs to be written."
    write(*,*) ! blank line to make output look nice

    close(lun)
    return

999 nError = UK_ERROR
    write(eMessage,*) "Error opening control file "//trim(controlFile)

  END SUBROUTINE read_control
!
!*************************************************************************
!
  SUBROUTINE control_counts(nAreas,nOuts)

    implicit none
    character (len=99) :: line, word
    integer            :: nAreas, nOuts, ios, num_words
    integer, parameter :: lun = 1

    open(lun,file=controlFile,status='old',err=999)

    nAreas  = 0 ! initialize
    nOuts   = 0

99  format(a99)

    read(lun,99,iostat=ios) line
    call count_words(line,num_words)    ! read to first non-comment non-blank
    do while ((any(line(1:1) == comment) .or. num_words == 0) .and. ios == 0)
       read(lun,99,iostat=ios) line
       call count_words(line,num_words)
    end do

    do while (ios == 0)

       call get_word(line,1,word)
       call CUPPER( word )

       if (word == "AREA") nAreas = nAreas + 1
       if (word == "CON" .or. word == "DEP" .or. word == "VIS") &
            nOuts = nOuts + 1

       read(lun,99,iostat=ios) line ! read to the next non-comment,
       call count_words(line,num_words)     ! non-blank line
       do while (num_words == 0 .and. ios == 0)
          read(lun,99,iostat=ios) line
          call count_words(line,num_words)
       end do

    end do

    close(lun)
    return

999 nError = UK_ERROR
    write(eMessage,*) "Error opening control file "//trim(controlFile)

  END SUBROUTINE control_counts
!
!*************************************************************************
!
  SUBROUTINE read_class1_rec(nAreas, MaxRec)

    IMPLICIT NONE

    integer                 :: nAreas ! the number of Class I areas
    integer                 :: iArea  ! which Class I area to process
    integer                 :: MaxRec ! max( nRec(:) )
    integer                 :: ios, i, astat
    integer, parameter      :: lun = 10
    real                    :: xlon, xlat
    character (len=1)       :: slash, backslash
    logical                 :: lexist

    slash     = char(47) ! AKA "forward slash"
    backslash = char(92)

!
!--- Allocate the arrays to nrec
!
    allocate( nRec(nAreas), STAT=astat )
    if ( astat > 0 ) then
       nError = UK_ERROR
       WRITE(eMessage,*) 'Error allocating nRec with size ',nAreas*4
       return
    end if

    write(*,*)
    do iArea = 1,nAreas
!
!--- First extract the Class I area code from the filename (first 4 characters)
!
       i = index(class1Filename(iArea), slash, .true.)                 ! Linux
       if (i == 0) i = index(class1Filename(iArea), backslash, .true.) ! DOS
       i = i + 1                  ! works if it's a filename with no path too
       class1Name(iArea) = class1Filename(iArea)(i:i+3)
!
!--- Count the number of receptors, and allocate the global vars lat,lon
!
       inquire(file=TRIM(class1Filename(iArea)),exist=lexist)
       if (.not. lexist) goto 99
       write(*,*)'Reading receptor locations from ',TRIM(class1Filename(iArea))
       open(lun+iArea,file=class1Filename(iArea), status='old',err=99,iostat=ios)

       nRec(iArea) = 0
       read(lun+iArea,*,iostat=ios) xlon, xlat ! just make sure they're readable
       do while (ios == 0)
          nRec(iArea) = nRec(iArea) + 1
          read(lun+iArea,*,iostat=ios) xlon, xlat
       end do

    end do ! do iArea = 1,nAreas
    write(*,*)

    MaxRec = maxval(nRec(:))
!
!--- Allocate the arrays to hold the lat-lon for ALL Class I areas
!
    allocate( lat((iArea),MaxRec), lon((iArea),MaxRec), STAT=astat )
    if ( astat > 0 ) then
       nError = UK_ERROR
       WRITE(eMessage,*) 'Error allocating coordinates (lat-lon)'
       return
    end if
    lat = -999. ! fill with a "missing value" flag, since these will be sparse
    lon = -999.
!
!--- Now read the Class I receptors, but don't read the elev (not used)
!
    do iArea = 1,nAreas

       rewind(lun+iArea)
       do i = 1,nRec(iArea)
          read(lun+iArea,*,iostat=ios) lon(iArea,i), lat(iArea,i)
          if (ios /= 0) goto 999
       end do
       close(lun+iArea)

    end do ! do iArea = 1,nAreas

    RETURN ! no errors

99  CONTINUE
    nError = UK_ERROR
    WRITE(eMessage,*) 'Error opening file ',trim(class1Filename(iArea))
    RETURN

999 CONTINUE
    nError = UK_ERROR
    WRITE(eMessage,*) 'Error ',ios,'reading file ',trim(class1Filename(iArea))
    RETURN

  END SUBROUTINE read_class1_rec
!
!*************************************************************************
!
  SUBROUTINE init_visibility(nAreas)

    IMPLICIT NONE

    integer                 :: nAreas ! number of Class I areas to consider
    integer                 :: iArea  ! loop over Class I Area
    integer                 :: mo     ! month
    integer                 :: ios, astat
    integer, parameter      :: lun = 9

    real                    :: TotSulfate, LgSulfate, SmSulfate
    real                    :: TotNitrate, LgNitrate, SmNitrate
    real                    :: TotOrganic, LgOrganic, SmOrganic

    character (len=256)     :: line   ! one line from a FLAG(2010) CSV table
    character (len=99)      :: word   ! one word from a FLAG(2010) CSV table

    logical                 :: lexist

!--- Allocate the FLAG(2010) table containers

    allocate( NatCond(nAreas,8), STAT=astat )
    IF( astat > 0 )THEN
       nError = UK_ERROR
       WRITE(eMessage,*) 'Error allocating NatCond'
       GO TO 9999
    END IF
    NatCond = 0.
    allocate( LargeRH(nAreas,12), STAT=astat )
    IF( astat > 0 )THEN
       nError = UK_ERROR
       WRITE(eMessage,*) 'Error allocating LargeRH'
       GO TO 9999
    END IF
    LargeRH = 0.
    allocate( SmallRH(nAreas,12), STAT=astat )
    IF( astat > 0 )THEN
       nError = UK_ERROR
       WRITE(eMessage,*) 'Error allocating SmallRH'
       GO TO 9999
    END IF
    SmallRH = 0.
    allocate( SeaSalt(nAreas,12), STAT=astat )
    IF( astat > 0 )THEN
       nError = UK_ERROR
       WRITE(eMessage,*) 'Error allocating SeaSalt'
       GO TO 9999
    END IF
    SeaSalt = 0.

    allocate( Bext_nat(nAreas,12), STAT=astat )
    IF( astat > 0 )THEN
       nError = UK_ERROR
       WRITE(eMessage,*) 'Error allocating Bext_nat'
       GO TO 9999
    END IF
    Bext_nat = 0.

!--- Read the FLAG look-up tables and extract the values for each Class I area

    NatCondFile = "flag.2010.table.6.annual.csv"
    CALL ADDPATH(NatCondFile,trim(flagDir))
    INQUIRE( file=NatCondFile,EXIST=lexist )
    if (.not. lexist) then
       nError = UK_ERROR
       WRITE(eMessage,*) 'Error opening ',trim(NatCondFile)
       GO TO 9999
    end if
    write(*,*) "Reading ",trim(NatCondFile)          ! table 5 (best 20%)
    open(lun, file=NatCondFile, status='old',err=91) ! or table 6 (annual)
    read(lun,*,iostat=ios) ! header line
    do while (ios == 0)
       read(lun,'(A256)',iostat=ios) line
       call get_csv(line,1,word)
       do iArea = 1, nAreas
          if (trim(word) == trim(class1Name(iArea))) then
             call get_csv(line, 3,word)
             read(word,*) NatCond(iArea,1) ! (NH4)2SO4
             call get_csv(line, 4,word)
             read(word,*) NatCond(iArea,2) ! NH4NO3
             call get_csv(line, 5,word)
             read(word,*) NatCond(iArea,3) ! Organic matter
             call get_csv(line, 6,word)
             read(word,*) NatCond(iArea,4) ! Elem carbon
             call get_csv(line, 7,word)
             read(word,*) NatCond(iArea,5) ! Soil
             call get_csv(line, 8,word)
             read(word,*) NatCond(iArea,6) ! Coarse Mass
             call get_csv(line, 9,word)
             read(word,*) NatCond(iArea,7) ! Sea Salt
             call get_csv(line,10,word)
             read(word,*) NatCond(iArea,8) ! Rayleigh
          end if
       end do
    end do
    close(lun)

    do iArea = 1, nAreas
       if ( NatCond(iArea,1) == 0.) then
         nError = UK_ERROR
         WRITE(eMessage,*)"Error in init_visibility(): can't find Class I area '",  &
               trim(class1Name(iArea)),"' in NatCond file ",trim(NatCondFile)
         GO TO 9999
       endif
    end do

    LargeRHFile = "flag.2010.table.7.FlargeRH.csv"
    CALL ADDPATH(LargeRHFile,trim(flagDir))
    INQUIRE( file=LargeRHFile,EXIST=lexist )
    if (.not. lexist) then
       nError = UK_ERROR
       WRITE(eMessage,*) 'Error opening ',trim(LargeRHFile)
       GO TO 9999
    end if
    write(*,*) "Reading ",trim(LargeRHFile)
    open(lun, file=LargeRHFile, status='old',err=92) ! table 7 (large RH factors)
    read(lun,*,iostat=ios) ! header line
    do while (ios == 0)
       read(lun,'(A256)',iostat=ios) line
       call get_csv(line,1,word)
       do iArea = 1, nAreas
          if (word == class1Name(iArea)) then
             do mo = 1, 12 ! read 12 months of LargeRH values
                call get_csv(line, 2+mo,word)
                read(word,*) LargeRH(iArea,mo)
             end do
          end if
       end do
    end do
    close(lun)

    SmallRHFIle = "flag.2010.table.8.FsmallRH.csv"
    CALL ADDPATH(SmallRHFIle,trim(flagDir))
    INQUIRE( file=SmallRHFile,EXIST=lexist )
    if (.not. lexist) then
       nError = UK_ERROR
       WRITE(eMessage,*) 'Error opening ',trim(SmallRHFile)
       GO TO 9999
    end if
    write(*,*) "Reading ",trim(SmallRHFile)
    open(lun, file=SmallRHFile, status='old',err=93) ! table 8 (small RH factors)
    read(lun,*,iostat=ios) ! header line
    do while (ios == 0)
       read(lun,'(A256)',iostat=ios) line
       call get_csv(line,1,word)
       do iArea = 1, nAreas
          if (word == class1Name(iArea)) then
             do mo = 1, 12 ! read 12 months of SmallRH values
                call get_csv(line, 2+mo,word)
                read(word,*) SmallRH(iArea,mo)
             end do
          end if
       end do
    end do
    close(lun)

    SeaSaltFile = "flag.2010.table.9.FseasaltRH.csv"
    CALL ADDPATH(SeaSaltFile,trim(flagDir))
    INQUIRE( file=SeaSaltFile,EXIST=lexist )
    if (.not. lexist) then
       nError = UK_ERROR
       WRITE(eMessage,*) 'Error opening ',trim(SeaSaltFile)
       GO TO 9999
    end if
    write(*,*) "Reading ",trim(SeaSaltFile)
    open(lun, file=SeaSaltFile, status='old',err=94) ! table 9 (sea salt factors)
    read(lun,*,iostat=ios) ! header line
    do while (ios == 0)
       read(lun,'(A256)',iostat=ios) line
       call get_csv(line,1,word)
       do iArea = 1, nAreas
          if (word == class1Name(iArea)) then
             do mo = 1, 12 ! read 12 months of SeaSalt values
                call get_csv(line, 2+mo,word)
                read(word,*) SeaSalt(iArea,mo)
             end do
          end if
       end do
    end do
    close(lun)
!
!--- Calculate the monthly natural background extinction, bext_natcond
!
    do iArea = 1, nAreas  ! loop over Class I Areas
       do mo = 1, 12      ! 12 months
!
!--- Calculate Total Sulfate, and large sulfate fraction, according to
!    Figure 5 from FLAG (2010).  Ditto for Nitate and Organic Mass
!
          TotSulfate = NatCond(iArea,1)
          if (TotSulfate < 20.) then ! 20 ug/m3
             LgSulfate = TotSulfate/20.*TotSulfate
          else
             LgSulfate = TotSulfate
          end if
          SmSulfate = TotSulfate - LgSulfate

          TotNitrate = NatCond(iArea,2)
          if (TotNitrate < 20.) then ! 20 ug/m3
             LgNitrate = TotNitrate/20.*TotNitrate
          else
             LgNitrate = TotNitrate
          end if
          SmNitrate = TotNitrate - LgNitrate

          TotOrganic = NatCond(iArea,3)
          if (TotOrganic < 20.) then ! 20 ug/m3
             LgOrganic = TotOrganic/20.*TotOrganic
          else
             LgOrganic = TotOrganic
          end if
          SmOrganic = TotOrganic - LgOrganic

!--- Calculate natural background visibility obscuration, by Area and Month

          Bext_nat(iArea,mo) =                              &
               2.2 * SmallRH(iArea,mo) * SmSulfate        + & ! small sulfate
               4.8 * LargeRH(iArea,mo) * LgSulfate        + & ! large sulfate
               2.4 * SmallRH(iArea,mo) * SmNitrate        + & ! small nitrate
               5.1 * LargeRH(iArea,mo) * LgNitrate        + & ! large nitrate
               2.8 *                     SmOrganic        + & ! sm org mass
               6.1 *                     LgOrganic        + & ! lg org mass
               10. *                     NatCond(iArea,4) + & ! elem. carbon
               1.0 *                     NatCond(iArea,5) + & ! fine soil
               0.6 *                     NatCond(iArea,6) + & ! coarse mass
               1.7 * SeaSalt(iArea,mo) * NatCond(iArea,7) + & ! sea salt
               NatCond(iArea,8)                           + & ! Rayleigh scat
               0.1755 * 0.                                    ! NO2

       end do  ! do mo = 1, 12 ! 12 months
    end do     ! do iArea = 1, nAreas  ! loop over Class I Areas

    write(*,*) ! blank line to make it look nice

9999 RETURN

91 WRITE(eMessage,*)"SciClass1Flag: error opening file ",trim(NatCondFile)
   GO TO 9998

92 WRITE(eMessage,*)"SciClass1Flag: error opening file ",trim(LargeRHFile)
   GO TO 9998
93 WRITE(eMessage,*)"SciClass1Flag: error opening file ",trim(SmallRHFile)
   GO TO 9998

94 WRITE(eMessage,*)"SciClass1Flag: error opening file ",trim(SeaSaltFile)
   GO TO 9998

9998 nError = UK_ERROR

     RETURN

  END SUBROUTINE init_visibility
!
!*************************************************************************
!
  SUBROUTINE get_word(line,which_word,word)

    IMPLICIT NONE

    character (len=*) :: line       ! input string
    character (len=*) :: word       ! output string
    integer           :: which_word ! which word to return

    integer           :: i,j, this_word, this_len  ! local variables
!
!-----Entry point
!
    i = 1
    this_len = len_trim(line)
    do while (i < this_len)
       if (any(line(i:i) == comment)) this_len = i - 1 ! ignore commented parts
       i = i + 1
    end do

    this_word = 0
    i = 1

    do while (this_word < which_word)
       if (any(line(i:i) == delim)) then ! advance to next non-space
          do while (any(line(i:i) == delim) .and. i < this_len)
             i = i + 1
          end do
       end if                           ! line(i:i) now the beginnig of next word
       if (i < this_len) then
          j = i + 1                     ! next character
       else
          j = i                         ! this character
       end if

       if (line(i:i) == '"') then       ! advance to matching double quote
          do while (line(j:j) /= '"' .and. j < this_len)
             j = j + 1                  ! advance to end of this word (next ")
          end do
       else if (line(i:i) == "'") then  ! advance to matching single quote
          do while (line(j:j) /= "'" .and. j < this_len)
             j = j + 1                  ! advance to end of this word (next ')
          end do
       else                             ! word doesn't start with a quote
          do while ((.not. any(line(j:j) == delim)) .and. j < this_len)
             j = j + 1                  ! advance to end of this word
          end do
          if (j > i .and. any(line(j:j) == delim)) j = j - 1
       end if
       this_word = this_word + 1
       if (this_word == which_word) then
          word = trim(line(i:j))
          if (word(1:1) == '"' .or. word(1:1) == "'") then
             word = word(2:(len_trim(word)-1)) ! remove the quotes
          end if
          return                        ! done, got the right word, exit
       else if (j == this_len) then
          write(*,*) "*** Error: Asked for word ",which_word,", but only ", &
               this_word," words found."
          write(*,*) "           Problematic line is:"
          write(*,*) trim(line)
          nError = UK_ERROR
          WRITE(eMessage,*)"*** Error: Asked for word ",which_word,", but only ", &
               this_word," words found."
          WRITE(eInform,*) "           Problematic line is:"//TRIM(line)

          RETURN
       end if
       i = j + 1                        ! start at next character
    end do

  END SUBROUTINE get_word
!
!*************************************************************************
!
  SUBROUTINE get_csv( line, which_word, word )

    implicit none
    character (len=*) :: line        ! input string
    character (len=*) :: word        ! output string
    character (len=1) :: comma = "," ! delimiter between words
    integer           :: which_word  ! which word to return
    integer           :: i,j, this_word, this_len  ! local variables
!
!-----Entry point
!
    this_len = len_trim(line)
    this_word = 0
    i = 1

    do while (this_word < which_word)
       if (line(i:i) == comma) then     ! advance to next non-comma
          do while (line(i:i) == comma .and. i < this_len)
             i = i + 1
          end do
       end if                           ! line(i:i) now the beginnig of next word
       if (i < this_len) then
          j = i + 1                     ! next character
       else
          j = i                         ! this character
       end if

       if (line(i:i) == '"') then       ! advance to matching double quote
          do while (line(j:j) /= '"' .and. j < this_len)
             j = j + 1                  ! advance to end of this word (next ")
          end do
       else if (line(i:i) == "'") then  ! advance to matching single quote
          do while (line(j:j) /= "'" .and. j < this_len)
             j = j + 1                  ! advance to end of this word (next ')
          end do
       else                             ! word doesn't start with a quote
          do while ((.not. line(j:j) == comma) .and. j < this_len)
             j = j + 1                  ! advance to end of this word
          end do
          if (j > i .and. line(j:j) == comma) j = j - 1
       end if
       this_word = this_word + 1
       if (this_word == which_word) then
          word = trim(line(i:j))
          if (word(1:1) == '"' .or. word(1:1) == "'") then
             word = word(2:(len_trim(word)-1)) ! remove the quotes
          end if
          return                        ! done, got the right word, exit
       else if (j == this_len) then
          write(*,*) "*** Error: Asked for word ",which_word,", but only ", &
               this_word," words found."
          write(*,*) "           Problematic line is:"
          write(*,*) trim(line)
          stop
       end if
       i = j + 1                        ! start at next character
    end do

  END SUBROUTINE get_csv
!
!*************************************************************************
!
  SUBROUTINE count_words(line,num_words)
!
!------------------------------------------------------------------------------
!     MESOSCALE MODEL INTERFACE PROGRAM (MMIF)
!     VERSION 3.0 2013-09-30
!
!     This subroutine counts the number of words in a line.
!
!     Development History:
!     2013-09-05  Original Development (ENVIRON International Corp.)

    implicit none

    character (len=*) :: line       ! input string
    integer           :: num_words  ! the number of words found

    integer           :: i,j,this_len
!
!-----Entry point
!
    i = 1
    this_len = len_trim(line)
    do while (i < this_len)
       if (any(line(i:i) == comment)) this_len = i - 1 ! ignore commented parts
       i = i + 1
    end do

    num_words = 0
    i = 1

    do while (i < this_len)
       if (any(line(i:i) == delim)) then ! advance to next non-space
          do while (any(line(i:i) == delim) .and. i < this_len)
             i = i + 1
          end do
          ! what if there's a delimiter at the end of the line?
          if (i == this_len .and. any(line(i:i) == delim) ) return
       end if                           ! line(i:i) now the beginnig of nextword

       j = i + 1                        ! next character
       if (line(i:i) == '"') then       ! advance to matching double quote
          do while (line(j:j) /= '"' .and. j < this_len)
             j = j + 1                  ! advance to end of this word (next ")
          end do
       else if (line(i:i) == "'") then  ! advance to matching single quote
          do while (line(j:j) /= '"' .and. j < this_len)
             j = j + 1                  ! advance to end of this word (next ')
          end do
       else                             ! word doesn't start with a quote
          do while ((.not. any(line(j:j) == delim)) .and. j < this_len)
             j = j + 1                  ! advance to end of this word
          end do
       end if
       num_words = num_words + 1
       i = j                            ! start at next character
    end do

  END SUBROUTINE count_words
!
!*************************************************************************
!
  SUBROUTINE usage

    write(*,*) "Usage: class1answer [-h | --help] [--sample] [-I:scpuff.ini] -i ctrl.inp"
    write(*,*) "Options:"
    write(*,*) "   -i file    Specify control file to use"
    write(*,*) "   -I:file    Specify the location of scipuff.ini to use"
    write(*,*) "   --sample   Print a sample control file to the screen"
    write(*,*) "   --help     Print this help message"
    write(*,*) "   -h         Print this help message"

    stop

  END SUBROUTINE usage
!
!*************************************************************************
!
  SUBROUTINE sample_control_file

1   format(a)

    write(*,1) "; Comment characters are #, ;, and !.  Blank lines are ignored."
    write(*,1) "; Keywords are case-insenstive, filenames are not."
    write(*,1)
    write(*,1) "# Download the Class I area receptors from the following URL:"
    write(*,1) "# http://www.nature.nps.gov/air/maps/receptors/index.cfm"
    write(*,1)
    write(*,1) "# Class I areas to include in this run (repeat AREA keyword)"
    write(*,1) "Area    /home/data/ClassI/meve-recep.dat   ! Mesa Verde National Park"
    write(*,1) "AREA    /home/data/ClassI/wemi-recep.dat   ! Weminuche Wilderness"
    write(*,1)
    write(*,1) "# Outputs:"
    write(*,1) "con  no2    1 yr  H1H   max_concentrations.csv"
    write(*,1) "con  no2    1 hr  >7.5  1hr_no2_over_sils.csv"
    write(*,1) "con  pm10  24 hr  H1H   max_concentrations.csv"
    write(*,1) "con  pm10   1 yr  H1H   max_concentrations.csv"
    write(*,1) "con  pm25  24 hr  H1H   max_concentrations.csv"
    write(*,1) "con  pm25   1 yr  H1H   max_concentrations.csv"
    write(*,1) "con  so2    3 hr  H1H   max_concentrations.csv"
    write(*,1) "con  so2   24 hr  H1H   max_concentrations.csv"
    write(*,1) "con  so2    1 yr  H1H   max_concentrations.csv"
    write(*,1) "dep  N      1 yr  H1H   max_deposition.csv"
    write(*,1) "dep  S      1 yr  H1H   max_deposition.csv"
    write(*,1) "vis  dBext 24 hr  H8H   visibility_results.csv"
    write(*,1) "vis  dBext 24 hr  >5    visibility_exceedences.csv"
    write(*,1) "vis  dBext 24 hr  >10   visibility_exceedences.csv"

    stop

  END SUBROUTINE sample_control_file
!
!*************************************************************************
!
  integer function find_rank(val,iArea,iRec,iOut,rank,debug)

    integer iArea, iRec, iOut, rank, r
    real    val
    logical debug

    r = 1
    if (debug) then
       print 10,"iRec,r,rank,val,max_val:",iRec,r,rank,val,  &
            max_val(iArea,iRec,iOut,r)
10     format(a,3i4,2g18.8)
       call flush(6)
    endif
    do while (val < max_val(iArea,iRec,iOut,r) .and. r < rank)
       if (debug) then
          print*,"r,val,max_val:",r,val,max_val(iArea,iRec,iOut,r)
          call flush(6)
       endif
       r = r + 1
    end do
    if (val < max_val(iArea,iRec,iOut,r)) r = r + 1 ! correct for r == rank
    if (debug) print*,"found rank: ",r
    find_rank = r

    return
  end function find_rank
!
!*************************************************************************
!
  SUBROUTINE insert(val, iArea, iRec, iOut, rank, ind, iymdh) ! real

! insert the value val in max_val(,) as the ind-th value, for a given iRec

    integer iArea, iRec, iOut, rank, r, ind
    integer iymdh  ! current time in YYYYMMDDHH format
    real    val

    if (ind > rank) return    ! nothing to be done

    if (ind < rank) then
       do r = rank, ind+1, -1 ! move those lower than ind down one notch
          max_val(iArea,iRec,iOut,r)  = max_val(iArea,iRec,iOut,r-1)
          max_time(iArea,iRec,iOut,r) = max_time(iArea,iRec,iOut,r-1)
       end do
    end if

    max_val(iArea,iRec,iOut,ind)  = val
    max_time(iArea,iRec,iOut,ind) = iymdh

    return
  END SUBROUTINE insert
!
!*************************************************************************
!
  REAL FUNCTION ppm2ugm3(Species)
!
!---- Convert from PPM (gasses) to ug/m3 (micro-grams per cubic meter).
!     Here are some verification values, taken from EPA definitions, which
!     are not necessarily the most accurate scientifically...
!          NO :   0.1 PPM = 123 ug/m3
!          NO2: 0.053 PPM = 100 ug/m3
!          NO2:   0.1 PPM = 188 ug/m3
!          SO2: 0.075 PPM = 196 ug/m3
!          CO :     9 PPM =  40 mg/m3 (millegrams, not micrograms)
!          O3 : 0.075 PPM = 147 ug/m3
!          Pb :     1 PPM = 8.5 mg/m3 (millegrams, not micrograms)
!     See http://www.skcinc.com/converter/converter.asp, for example.
!
    character (len=*) :: Species ! input:  "NO" or "NO2" or "O3", etc.
    real              :: factor  ! output: conversion factor
    real              :: MW      ! molecular weight of species

    if (Species == "") then
       MW = -999. ! flag for not found
    elseif (Species == "N")   then
       MW = 14.01
    elseif (Species == "N2")   then
       MW = 28.01
    elseif (Species == "NO") then
       MW = 30.01
    elseif (Species == "NO2") then
       MW = 46.01
    elseif (Species == "N2O") then
       MW = 44.01
    elseif (Species == "N2O5") then
       MW = 95.01
    elseif (Species == "HNO3") then
       MW = 63.01
    elseif (Species == "HONO") then
       MW = 47.01
    elseif (Species == "PNA") then ! HNO4
       MW = 79.01
    elseif (Species == "PAN") then ! CH3COOONO2
       MW = 121.05
    elseif (Species == "NH3")  then
       MW = 17.03
    elseif (Species == "S") then
       MW = 32.06
    elseif (Species == "SO2") then
       MW = 64.06
    elseif (Species == "CO")  then
       MW = 28.01
    elseif (Species == "CO2")  then
       MW = 44.01
    elseif (Species == "O3")  then
       MW = 48.00
    elseif (Species == "O2")  then
       MW = 32.00
    elseif (Species == "O")   then
       MW = 16.00
    elseif (Species == "H2O")  then
       MW = 18.02
    elseif (Species == "Pb") then
       MW = 207.20
    else
       MW = -999. ! flag for not found
    end if

    if (MW == -999.) then
       PPM2ugm3 = 1.                 ! default to no conversion
    else
       PPM2ugm3 = 1000. * MW / 24.45 ! 1000 is millegrams to micrograms
    endif

    RETURN

  END FUNCTION ppm2ugm3
!
!*************************************************************************
!
  SUBROUTINE find_units_conv(units_in, units_out, avg_time, units_conv)

    implicit none
! inputs:
    character (len=*) :: units_in, units_out
    integer           :: avg_time   ! in seconds

! output:
    real              :: units_conv ! output conversion factor

! local
    integer           :: i,j,k,l,m
    character (len=9) :: numer_in, numer_out, denom_in, denom_out, time_in

    i = index(units_in, "-")
    j = index(units_in, "/")
    k = len_trim(units_in)

    l = index(units_out,"/")
    m = len_trim(units_out)

    units_conv = 1.   ! default

    numer_out = units_out(1:l-1)
    denom_out = units_out(l+1:m)

    if ( i > 0 ) then ! either kg-s/m3 or g-s/m3, Tracer (not mult-component)
       time_in  = units_in(i+1:j-1)
       if ( TRIM(time_in) == "s" ) then
          units_conv = units_conv / avg_time ! s/s
       else if ( TRIM(time_in) == "hr" ) then
          units_conv = units_conv / avg_time * 3600 ! hr/s * s/hr
       else
          write(*,*) "Error in find_units_conv(): Unknown time units:", &
               trim(time_in)
          write(*,*) "  units_in  = ",trim(units_in)
          write(*,*) "  units_out = ",trim(units_out)
          stop
       endif
       numer_in = units_in(1:i-1)
       denom_in = units_in(j+1:k)
    else              ! any of kg/m2, g/m2, ug/m3 or ug/m2
       numer_in = units_in(1:j-1)
       denom_in = units_in(j+1:k)
    end if

    if ( numer_in == numer_out) then
       ! nothing to do here
    else if ( numer_in == "kg" .and. numer_out == "ug") then
       units_conv = units_conv * 1.e9 ! ug/kg
    else if ( numer_in == "g"  .and. numer_out == "ug") then
       units_conv = units_conv * 1.e6 ! ug/g
    end if

    if ( denom_in == denom_out ) then
       ! nothing to do here
    else if ( denom_in == "m2" .and. denom_out == "ha") then
       units_conv = units_conv * 1.e4 ! m^2/hectare
    else if ( denom_in == "ha" .and. denom_out == "m2") then
       units_conv = units_conv / 1.e4 ! hectare/m^2
    else if ( denom_in == "m2" .and. denom_out == "m3") then
       write(*,*) "Error in find_units_conv(): tried to conver /m2 to /m3"
       write(*,*) "  units_in  = ",trim(units_in)
       write(*,*) "  units_out = ",trim(units_out)
       stop
    end if

    RETURN

  END SUBROUTINE find_units_conv
!
!*************************************************************************
!
  LOGICAL FUNCTION leap_yr(Yr) ! Y2K correct.

    IMPLICIT NONE

    integer yr

! If integer divide and float divide give the same result, then it's
! evenly divisible.  Could have used mod() here, but I think this is
! more portable.  Might give problems on some old Pentium systems
! with that pesky 4./2. = 1.9999998 error.

! The year is 365.2422 days long.
! Rules for leap years: if year is evenly divisible by 4, then it's
! a leap year, except those evenly divisible by 100, but there is a
! leap year in those evenly divisible by 400.  This will give a mean
! year 365.2425 days long.  Error is .0003 days (25.92 seconds) so
! it will take 3333 years for the calendar to get off by one day.

    if ((float(Yr/4) .eq. Yr/4.) .and. ((float(Yr/100) .ne. Yr/100.) &
         .or. (float(Yr/400) .eq. Yr/400.))) then
       leap_yr = .true.
    else
       leap_yr = .false.
    end if

    RETURN

  END FUNCTION leap_yr

!======================================================================
! DataInitialization
!======================================================================

  SUBROUTINE InitSciClass1Flag()

! Intializes all common variables

    USE tooluser_fd
    USE Extract_fi
    USE GetTimes_fi
    USE contour_fd

    IMPLICIT NONE

!--- Initialize

    callerID = USER_ID
    maxTry   = 5
    nFields  = 0

!Error handling

    nError      = NO_ERROR
    eMessage    = ' '
    eInform     = ' '
    eAction     = ' '
    eRoutine    = ' '

!Flags

    extractMethod = GRID

    Project%name = ' '
    Project%path = ' '

!Plot type

    plotType%type     = HP_MEAN
    plotType%data     = 0.0

!Contour

    contourHead%number    = 1
    contourHead%DrawMode  = PLOT_OFF
    contourHead%LabelMode = PLOT_OFF
    contourHead%Scale     = 1.0
    contourHead%Unit      = 'default'
    contourMode           = CLOSE_CONTOUR

!Plot category data

    sliceHeight    = 0.0             !Horizontal slice height
    sliceXmin      = DEF_VAL_R       !Minimum X slice value - vertical slice
    sliceXmax      = DEF_VAL_R       !Maximum X slice value - vertical slice
    sliceYmin      = DEF_VAL_R       !Minimum Y slice value - vertical slice
    sliceYmax      = DEF_VAL_R       !Maximum Y slice value - vertical slice
    sliceZmin      = 0.0             !Minimum Z slice value - vertical slice
    sliceZmax      = 2500.           !Maximum Z slice value - vertical slice
    sliceNz        = 21              !Slice Z points - Vertical slice

    iround         = -999
    riskLevel      = 10.0

!Plot data sizes
    nClass    = 0
    nChoice   = 0
    nKind     = 0
    nTimePuff = 0
    nTimeSrf  = 0
    nTimeMet  = 0
    nTimeOut  = 0
    nTable    = 0
    nCol      = 0
    nRow      = 0
    nTriangle = 0
    nNode     = 0
    nPoint    = 0
    nLine     = 0

!Field max values

    FldMax = -HUGE(1.0)
    FldMin = HUGE(1.0)

!Command line strings

    ini_file    = ''
    PrjName     = ''

!Grid

    maxGrd         = 0                      !Maximum grid cells
    GridType       = NATIVE
    ContourType    = AUTO
    extractMethod  = GRID

    xMin =  HUGE(1.0)
    xMax = -HUGE(1.0)
    yMin = xMin
    yMax = xMax
    zMin = sliceZmin
    zMax = SliceZmax
    dx   = 0.0
    dy   = 0.0
    dz   = 0.0

    nx = 0
    ny = 0
    nz = 0

!Unit numbers
    lun_in  = 5
    lun_out = 6
    lun_grd = 0
    lun_ter = 0

    Field%category   = HP_SSLICE              !Field category index
    Field%class      = -1
    Field%choice     = 0
    Field%kind       = 0
    Field%timeID     = 0
    Field%userTime   = 0.0
    Field%iNotUsed   = SCIPoff
    Field%maxCells   = 25000
    Field%maxLev     = 99
    Field%fldLev     = 99
    Field%resolution = 1.0
    Field%interpType = 0.0
    Field%units      = 'XX'
    Field%project    = ' '
    Field%path       = ' '

! Fill up SCIPFieldCoordinateT structure
    Field%coordinate%mode                 = HD_LATLON
    Field%coordinate%UTMZone              = NOT_SET_I
    Field%coordinate%reference%x          = NOT_SET_R
    Field%coordinate%reference%y          = NOT_SET_R
    Field%coordinate%reference%lat        = NOT_SET_R
    Field%coordinate%reference%lon        = NOT_SET_R
    Field%coordinate%vertSlice%resolution = NOT_SET_I
    Field%coordinate%vertSlice%startPt%x  = NOT_SET_R
    Field%coordinate%vertSlice%startPt%y  = NOT_SET_R
    Field%coordinate%vertSlice%endPt      = Field%coordinate%vertSlice%startPt
    Field%coordinate%horzSlice%height     = NOT_SET_R
    Field%coordinate%horzSlice%mode       = NOT_SET_I

    NULLIFY( TimePuff )
    NULLIFY( TimeSrf )
    NULLIFY( TimeMet )
    NULLIFY( TimeOut )

    contourExt = '.txt'

    RETURN
  END SUBROUTINE InitSciClass1Flag

!======================================================================
! ToolInitialization
!======================================================================

  SUBROUTINE InitToolSciClass1Flag()

    USE basic_fd
    USE tooluser_fd
    USE SCIPtool
    USE Extract_fi
    USE GetTimes_fi

    IMPLICIT NONE

    INTEGER irv
    INTEGER request
    INTEGER(LEN_ADDRESS) CallBackAddress
    INTEGER idefault
    LOGICAL lExist

    CHARACTER(32)  mode
    CHARACTER(128) string1,string2,string3,string4

    TYPE( fileNameT ) :: file

    CHARACTER(PATH_MAXLENGTH), EXTERNAL :: StripNull, AddNull

    INTEGER,              EXTERNAL :: sysGetProfileString
    INTEGER,              EXTERNAL :: sysGetProfileInt
    INTEGER(LEN_ADDRESS), EXTERNAL :: ADDRESSOF
    INTEGER,              EXTERNAL :: sppCallback

    CallBackAddress = ADDRESSOF( sppCallback )

    request =  0

!--- Start up : Initialize

    IF( LEN_TRIM(ini_file) <= 0 )THEN
       WRITE(6,*)
       WRITE(6,'(''PATH for scipuff.ini file(press Enter for default): '',$)')
       READ(lun_in,'(A)')string2
       IF ( LEN_TRIM(string2) == 0 ) THEN
          ini_file = 'scipuff.ini'
          INQUIRE( file=ini_file,EXIST=lexist )
          IF( .NOT.lexist )THEN
             ini_file = 'DEFAULT'
          END IF
       ELSE
          ini_file = TRIM(ADJUSTL(string2))//'/scipuff.ini'
          INQUIRE( file=ini_file,EXIST=lexist )
          IF( .NOT.lexist )THEN
             WRITE(6,*)'Cannot find ',TRIM(ini_file)
             ini_file = 'DEFAULT'
          END IF
       END IF
    ELSE
       ini_file = ADJUSTL(ini_file)
    END IF

    IF( TRIM(ini_file) == 'DEFAULT' )THEN
       WRITE(6,*)'Using default scipuff.ini'
       ini_file = 'scipuff.ini'
    END IF

    INQUIRE(file=ini_file,EXIST=lexist)
    IF( .NOT.lexist )THEN
       nError      = UK_ERROR
       eMessage    = 'Unable to initialize the SCIPtool library'
       eInform     = 'Specified initialization file not found'
       eAction     = 'INI file='//TRIM(ini_file)
       eRoutine    = 'ToolInitialzaion'
       GOTO 9999
    END IF
    WRITE(6,*)
    WRITE(6,*)'Initializing SCIPtool from '//TRIM(ini_file)
    WRITE(6,*)

    string4 = AddNull( TRIM(ini_file) )
    string1 = AddNull( 'Paths' )
    string2 = AddNull( 'SciDataDir' )
    string3 = '.'
    irv = sysGetProfileString( string1,string2,string3,flagDir,string4 )
    IF( TRIM(flagDir) /= '.' )THEN
      string3 = flagDir
      flagDir = 'Flag2010'
      CALL ADDPATH(flagDir,string3)
      class1AreaDir = 'Class1Area'
      CALL ADDPATH(class1AreaDir,string3)
    ELSE
      class1AreaDir = '.'
    END IF
    WRITE(6,*)
    WRITE(6,*)'flagDir = ',TRIM(flagDir)
    WRITE(6,*)'class1AreaDir = ',TRIM(class1AreaDir)
    WRITE(6,*)

    string4 = AddNull( TRIM(ini_file) )
    string1 = AddNull( 'SCIPMode' )
    string2 = AddNull( 'GUIMode' )
    string3 = 'Standard'
    irv = sysGetProfileString( string1,string2,string3,mode,string4 )

    CALL cupper( mode )

    SELECT CASE( TRIM(mode) )
    CASE( 'STANDARD' )
       limit%met1D       = HUGE(1)
       limit%puffs       = 20000
       limit%surfaceGrid = 25000

    CASE( 'OPERATIONAL' )
       limit%met1D       = HUGE(1)
       limit%puffs       = 20000
       limit%surfaceGrid = 25000

    CASE( 'EXTENDED' )
       limit%met1D       = HUGE(1)
       limit%puffs       = 40000
       limit%surfaceGrid = 85000

    CASE( 'ULTIMATE' )
       limit%met1D       = HUGE(1)
       limit%puffs       = 60000
       limit%surfaceGrid = 100000

    CASE DEFAULT
       limit%met1D       = HUGE(1)
       limit%puffs       = 20000
       limit%surfaceGrid = 25000
       string2 = AddNull( 'MaxMet1D' )
       idefault = limit%met1D
       irv = sysGetProfileInt( string1,string2,idefault,limit%met1D,string4 )
       string2 = AddNull( 'MaxPuff' )
       idefault = limit%puffs
       irv = sysGetProfileInt( string1,string2,idefault,limit%puffs,string4 )
       string2 = AddNull( 'MaxGrid' )
       idefault = limit%surfaceGrid
       irv = sysGetProfileInt( string1,string2,idefault,limit%surfaceGrid,string4 )
    END SELECT

!==== Initialize PlotStub

    WRITE(6,'("Running SciClass1Flag in ",A," mode")')TRIM(mode)
    WRITE(6,'("  MaxPuff = ",I15)')limit%puffs
    WRITE(6,'("  MaxGrid = ",I15)')limit%surfaceGrid

    file%string = TRIM(ini_file)

    irv = SCIPInitTool( callerID,CallBackAddress,request,limit,file )
    IF( irv == SCIPfailure )THEN
       CALL toolError ('Failed to initialize the SCIPtool library' )
       GOTO 9999
    END IF

    lInit = .TRUE.

9999 CONTINUE

    RETURN
  END SUBROUTINE InitToolSciClass1Flag
!
!*************************************************************************
!
END MODULE class1Tools
!
!*************************************************************************
!
PROGRAM runSciClass1Flag

USE DefSize_fd
USE Extract_fi
USE tooluser_fd
USE SCIPtool
USE GetTimes_fi
USE sagdef_fd
USE sagstr_fd
USE PtrGrdStrItf
USE scipuff_fi, ONLY: create,ntypm,material,mat_mc
USE MPI_fi, ONLY: useMPI,myid
USE class1Tools
USE met_fi, ONLY:PrjCoord

IMPLICIT NONE

TYPE(an_output), allocatable, dimension(:) :: output
TYPE( SAGgrid_str ), POINTER               :: grdC, grdD

CHARACTER (len=8)                :: maxPath
CHARACTER (len=20)               :: units
CHARACTER (len=128)              :: string1
CHARACTER (len=PATH_MAXLENGTH)   :: dosFile, depFile
CHARACTER (len=256)              :: Command

INTEGER             :: nOuts            ! number of Output values to find
INTEGER             :: iOut             ! index for Output
INTEGER             :: nAreas           ! number of Class I Areas to process
INTEGER             :: iArea            ! index for Class I Areas
INTEGER             :: iRec             ! index for receptor within Class I areas
INTEGER             :: MaxRec           ! max( nRec(:) )
INTEGER             :: MaxRank          ! max( output(:)%rank )
INTEGER             :: avg_time         ! output frequency = integration time (s)
INTEGER             :: astat            ! short for alloc_stat
INTEGER             :: irv              ! return value for SCIP calls
INTEGER             :: ij               ! index of lat,lon where answer is found
INTEGER             :: iy,im,id,ih,iymdh! year, month, day, hour, YYYYMMDDHH
INTEGER             :: iClass, iChoice, iKind, iDOS, iDEP
INTEGER             :: i, j, iTime, ios, iArg
INTEGER             :: cGrdI, dGrdI
INTEGER             :: lun_dos,lun_dep,lun
INTEGER             :: numCommand
INTEGER(2)          :: ios2

LOGICAL             :: lexist, lopen    ! does a file exist, open ?
LOGICAL             :: isMULT = .false. ! multi-component or tracer run?

REAL                :: con_units_conv, dep_units_conv ! units conversion factors
REAL                :: answer           ! what we seek

INTEGER, EXTERNAL   :: SCIPGetPrjCoord
INTEGER, EXTERNAL   :: SAG_InitList, SAG_ClearList, SAG_SetSpecialValue
INTEGER, EXTERNAL   :: SAG_OpenID, SAG_CloseID, SAG_RmvGrdStr
INTEGER, EXTERNAL   :: SAG_NewGrdStr, SAG_InitGridID, SAG_ReadHeaderID
INTEGER, EXTERNAL   :: sysNumArgs
INTEGER, EXTERNAL   :: sysGetArg

CHARACTER(128), EXTERNAL :: AddExtension

!*************************************************************************
!--- BEGIN MAIN PROGRAM
!*************************************************************************

useMPI = .FALSE.
myid   = 0
create = .FALSE.

!--- Start up : Initialize common

CALL InitSciClass1Flag()

!--- Get the input file from the command line

controlFile = 'sciClass1Flag.inp' ! default
numCommand  = sysNumArgs()
iArg = 0
do while (iArg < numCommand)
   iArg = iArg + 1
   ios2 = sysGetArg(iArg,Command)
   if (Command(1:3) == '-I:') then
      read(Command(4:len_trim(Command)),'(a256)') ini_file
   elseif (trim(Command) == '-i') then
      iArg = iArg + 1
      ios2 = sysGetArg(iArg,controlFile)
   else if (trim(Command) == '-h' .or. trim(Command) == '--help') then
      call usage
   else if (trim(Command) == '--sample') then
      call sample_control_file
   else
      write(*,*) "Error: unknown switch ",trim(Command)
      call usage
   end if
end do

!--- Start up : Write Screen Header

WRITE(maxPath,'(I0)') PATH_MAXLENGTH

WRITE(6,'(/,A)')'Program : SciClass1Flag v'//TRIM(CODE_VERSION)//' &
     &(maxPath='//ADJUSTL(TRIM(maxPath))//')'

irv = SCIPGetVersionString(USER_ID,toolString)
IF( irv == SCIPsuccess )THEN
  irv = SCIPGetPathMaxLength()
  IF( irv == PATH_MAXLENGTH )THEN
    WRITE(6,'(A)')'   uses : SCIPtool v'//TRIM(toolString%string)
  ELSE
    nError = UK_ERROR
    eMessage = 'Error:SCIPtool PATH_MAXLENGTH mismatch'
    eInform ='SCIPUFFPostProcess='//ADJUSTL(TRIM(maxPath))//' : SCIPtool='
    WRITE(maxPath,'(I0)')irv
    eInform =TRIM(eInform)//ADJUSTL(TRIM(maxPath))
    GO TO 9999
  END IF
ELSE
  nError = UK_ERROR
  eMessage = 'Error checking SCIPtool version'
  GO TO 9999
END IF

!--- Start up : Initialize SCIPtool library

CALL InitToolSciClass1Flag()
IF ( nError /= NO_ERROR ) GOTO 9999

!--- Count the number of occurences of a few keywords

call control_counts(nAreas,nOuts)

!--- Allocate the arrays, based on the counts

allocate( output(nOuts),STAT=astat )
IF( astat > 0 )THEN
  nError = UK_ERROR
  WRITE(eMessage,*) 'Error allocating outputs(nOuts)'
  GO TO 9999
END IF
output(1:nOuts)%exceeds = HUGE(0.)  ! largest number possible
output(1:nOuts)%rank    = 1         ! default to output rank 1 (H1H)
output(1:nOUts)%lun     = 0         ! default logical unit number

allocate( class1Name(nAreas), class1Filename(nAreas), STAT=astat )
IF( astat > 0 )THEN
  nError = UK_ERROR
  WRITE(eMessage,*) 'Error allocating class1Name/class1Filename'
  GO TO 9999
END IF

!--- Parse the control file

call read_control( nAreas, nOuts, output )
IF( nError /= NO_ERROR ) goto 9999

!--- Read in the lat-lon receptors for the Class I areas, allocates lat and lon

call read_class1_rec( nAreas, MaxRec)
IF( nError /= NO_ERROR ) goto 9999

!--- Initialize the visibility calcuations, if we'll need them

if ( any(output(:)%ConDepVis == "VIS") ) then
   call init_visibility(nAreas)
   IF( nError /= NO_ERROR ) goto 9999
end if

!--- First delete all the output files, so we start fresh

do iOut = 1, nOuts
   open (1, file=output(iOut)%outFile)
   close(1, status='delete')
end do

!--- Open the output files, and write headers

lun = 10
do iOut = 1, nOuts

   inquire(file=output(iOut)%outFile, exist=lexist)
   if (.not. lexist) then ! not found, must be new, write a header
      lun = lun + 1
      output(iOut)%lun = lun
      open(lun, file=output(iOut)%outFile, status="new",iostat=ios)
      if (ios /= 0) then
         nError   = UK_ERROR
         eRoutine = 'SciClass1Flag'
         write(eMessage,*) 'Error ',ios,'opening output file.'
         WRITE(eInform,*)  'File = ',TRIM(output(iOut)%outFile)
         goto 9999
      end if

      write(lun,*) &
           "Area, Species, lat, lon, Time, Period, Rank, Value"
   else
      output(iOut)%lun = lun
   end if

end do

! -- Get Plot Times

WRITE(*,*)'Reading project files for ',TRIM(ProjName)
CALL SplitName( ProjName,Project%name,Project%path ) !project name & path
Project%ID      = 0
Project%version = 0

!==== Action : Retrieve output times

irv = SCIPNumPlotTimes( CallerID,Project,nTimePuff,nTimeSrf,nTimeMet,nNotUsed )
IF( irv /= SCIPsuccess )THEN
   CALL toolError( 'Failed retrieving number of project plot times' )
   GOTO 9999
END IF

CALL allocatePlotTimes( )

irv = SCIPGetPlotTimes( callerID,Project,TimePuff,TimeSrf,TimeMet )
IF( irv /= SCIPsuccess )THEN
   CALL toolError( 'Failed retrieving project plot times' )
   GOTO 9999
END IF

TimeOut  => TimeSrf
nTimeOut = nTimeSrf

irv = SCIPNumPlotClasses( callerID,Project,nClass,nChoice,nKind )
IF( irv /= SCIPsuccess )THEN
   CALL toolError( 'Failed retrieving number of project plot classes' )
   GOTO 9999
END IF

CALL allocatePlotLists( )

irv = SCIPGetPlotClasses( callerID,Project,ClassStr,ChoiceStr,KindStr, &
     CatClassArray,ClassChoiceArray,ProjectCoordinate )
IF( irv /= SCIPsuccess )THEN
   CALL toolError( 'Failed retrieving project plot classes' )
   GOTO 9999
END IF

!--- Create receptor grid in project coordinates

ALLOCATE( xRec(nAreas,maxRec),yRec(nAreas,maxRec),STAT=irv )
xRec = 0.0
yRec = 0.0

DO iArea = 1,nAreas
  irv = SCIPGetPrjCoord( callerID,nRec(iArea),               &
       lon(iArea,1:nRec(iArea)),lat(iArea,1:nRec(iArea)),    &
       xRec(iArea,1:nRec(iArea)),yRec(iArea,1:nRec(iArea)) )
  IF( irv /= SCIPsuccess )THEN
    CALL toolError( 'Failed converting to project coordinates' )
    GOTO 9999
  END IF
END DO

!==== Action : Set Plot Field

Field%maxCells = limit%surfaceGrid
Field%project  = TRIM(Project%name)
Field%path     = TRIM(Project%path)

Field%class    = -1
Field%choice   = -1
Field%category = HP_SURF
PlotType%type  = HP_MEAN

!--- Find the iClass for the DOSAGE (con) file, call it iDOS

iDOS = -1
DO i = 1,nClass
   string1 = ClassStr(i)%string
   CALL CUPPER( string1 )
   IF( INDEX(string1,'SURFACE DOSAGE') > 0 )THEN
      iDOS = i
      EXIT
   END IF
END DO
IF( iDOS == -1 )THEN
  WRITE(*,*)'SURFACE DOSAGE not available'
  GOTO 9999
END IF
Field%class = iDOS

!--- Find the Choice, which determines multi-component vs. tracer runs

iChoice = -1
DO i = 1,nChoice
   string1 = ChoiceStr(i)%string
   CALL CUPPER( string1 )
   IF( INDEX(string1,':COMPONENTS') > 0 )THEN
      iChoice = i
      EXIT
   END IF
END DO

isMult = .FALSE.
IF( iChoice == -1 )THEN
   WRITE(*,*)'This is a Single-Component run'
   iChoice = 1
   nKind   = 1
ELSE
   WRITE(*,*)'This is a Multi-Component run'
   isMult = .TRUE.
   nKind  = ClassChoiceArray(iDOS,iChoice)%nkind
END IF
write(*,*) ! blank line after multi-single component message

Field%choice = iChoice
CALL setPlotClassData(.FALSE.,.FALSE.)

!--- Find the index of each chemical species we need to use later

allocate(iCon(nKind)) ! This triggers whether or not to process in class1con
iCon = 0              ! Fill with zeroes, non-zero means 'process this one'

!--- Initialize to -1 so we can detect if they are not found

iPSO4  = -1 ; iASO4K = -1 ; iPNO3  = -1 ; iANO3K = -1 ; iPSOA  = -1
iPOC   = -1 ; iPEC   = -1 ; iASOIL = -1 ; iACORS = -1 ; iNO2   = -1
iNO    = -1 ; iNO3   = -1 ; iN2O5  = -1 ; iHNO3  = -1 ; iHONO  = -1
iPNA   = -1 ; iPAN   = -1 ; iPANX  = -1 ; iNTR   = -1 ; iNH3   = -1
iPNH4  = -1 ; iANH4K = -1 ; iSO2   = -1 ; iSULF  = -1 ;

do iKind = 1, nKind

   if (isMULT) then                        ! multi-species output
      Field%kind = iKind + ClassChoiceArray(iDos,iChoice)%ikind - 1
   else                                    ! Only Tracer output
      Field%kind = 0
   endif

   if (isMult) then                        ! multi-species output
      string1 = KindStr(Field%kind)%string
   else                                    ! Only Tracer output
      string1 = ChoiceStr(iChoice)%string
   endif
   CALL CUPPER( string1 )

! These concentrations are to be averaged/ranked and output

   do iOut = 1, nOuts
      if (TRIM(string1) == output(iOut)%chemSpecies) iCon(iOut) = iKind
   end do

! These species are needed for visilbity

   if (TRIM(string1) == 'PSO4' ) iPSO4  = iKind ! small sulfate
   if (TRIM(string1) == 'ASO4K') iASO4K = iKind ! large sulfate
   if (TRIM(string1) == 'PNO3' ) iPNO3  = iKind ! small nitrate
   if (TRIM(string1) == 'ANO3K') iANO3K = iKind ! large nitrate
   if (TRIM(string1) == 'PSOA' ) iPSOA  = iKind ! small organic mass component
   if (TRIM(string1) == 'POC'  ) iPOC   = iKind ! small organic mass component
   if (TRIM(string1) == 'PEC'  ) iPEC   = iKind ! elemental carbon
   if (TRIM(string1) == 'ASOIL') iASOIL = iKind ! soil
   if (TRIM(string1) == 'ACORS') iACORS = iKind ! coarse mass
   if (TRIM(string1) == 'NO2'  ) iNO2   = iKind ! nitrogen dioxide

! These additional species are needed for deposition
! Note: NO2, PNO3, ANO3K, PSO4, ASO4K, already found above

   if (TRIM(string1) == 'NO'   ) iNO    = iKind
   if (TRIM(string1) == 'NO3'  ) iNO3   = iKind
   if (TRIM(string1) == 'N2O5' ) iN2O5  = iKind
   if (TRIM(string1) == 'HNO3' ) iHNO3  = iKind
   if (TRIM(string1) == 'HONO' ) iHONO  = iKind
   if (TRIM(string1) == 'PNA'  ) iPNA   = iKind
   if (TRIM(string1) == 'PAN'  ) iPAN   = iKind
   if (TRIM(string1) == 'PANX' ) iPANX  = iKind
   if (TRIM(string1) == 'NTR'  ) iNTR   = iKind
   if (TRIM(string1) == 'NH3'  ) iNH3   = iKind
   if (TRIM(string1) == 'PNH4' ) iPNH4  = iKind
   if (TRIM(string1) == 'ANH4K') iANH4K = iKind
   if (TRIM(string1) == 'SO2'  ) iSO2   = iKind
   if (TRIM(string1) == 'SULF' ) iSULF  = iKind

end do  ! do i = 1, nKind

!--- Warn the user about missing chemical species.  We'll just skip their
!    contributions

if ( iPSO4   == -1) write(*,*) "Warning: PSO4  not found in DOS/DEP files, ", &
	"answers may be wrong!!!"
if ( iASO4K  == -1) write(*,*) "Warning: ASO4K not found in DOS/DEP files, ", &
	"answers may be wrong!!!"
if ( iPNO3   == -1) write(*,*) "Warning: PNO3  not found in DOS/DEP files, ", &
	"answers may be wrong!!!"
if ( iANO3K  == -1) write(*,*) "Warning: ANO3K not found in DOS/DEP files, ", &
	"answers may be wrong!!!"
if ( iPSOA   == -1) write(*,*) "Warning: PSOA  not found in DOS/DEP files, ", &
	"answers may be wrong!!!"
if ( iPOC    == -1) write(*,*) "Warning: POC   not found in DOS/DEP files, ", &
	"answers may be wrong!!!"
if ( iPEC    == -1) write(*,*) "Warning: PEC   not found in DOS/DEP files, ", &
	"answers may be wrong!!!"
if ( iASOIL  == -1) write(*,*) "Warning: ASOIL not found in DOS/DEP files, ", &
	"answers may be wrong!!!"
if ( iACORS  == -1) write(*,*) "Warning: ACORS not found in DOS/DEP files, ", &
	"answers may be wrong!!!"
if ( iNO2    == -1) write(*,*) "Warning: NO2   not found in DOS/DEP files, ", &
	"answers may be wrong!!!"
if ( iNO     == -1) write(*,*) "Warning: NO    not found in DOS/DEP files, ", &
	"answers may be wrong!!!"
if ( iNO3    == -1) write(*,*) "Warning: NO3   not found in DOS/DEP files, ", &
	"answers may be wrong!!!"
if ( iN2O5   == -1) write(*,*) "Warning: N2O5  not found in DOS/DEP files, ", &
	"answers may be wrong!!!"
if ( iHNO3   == -1) write(*,*) "Warning: HNO3  not found in DOS/DEP files, ", &
	"answers may be wrong!!!"
if ( iHONO   == -1) write(*,*) "Warning: HONO  not found in DOS/DEP files, ", &
	"answers may be wrong!!!"
if ( iPNA    == -1) write(*,*) "Warning: PNA   not found in DOS/DEP files, ", &
	"answers may be wrong!!!"
if ( iPAN    == -1) write(*,*) "Warning: PAN   not found in DOS/DEP files, ", &
	"answers may be wrong!!!"
if ( iPANX   == -1) write(*,*) "Warning: PANX  not found in DOS/DEP files, ", &
	"answers may be wrong!!!"
if ( iNTR    == -1) write(*,*) "Warning: NTR   not found in DOS/DEP files, ", &
	"answers may be wrong!!!"
if ( iNH3    == -1) write(*,*) "Warning: NH3   not found in DOS/DEP files, ", &
	"answers may be wrong!!!"
if ( iPNH4   == -1) write(*,*) "Warning: PNH4  not found in DOS/DEP files, ", &
	"answers may be wrong!!!"
if ( iANH4K  == -1) write(*,*) "Warning: ANH4K not found in DOS/DEP files, ", &
	"answers may be wrong!!!"
if ( iSO2    == -1) write(*,*) "Warning: SO2   not found in DOS/DEP files, ", &
	"answers may be wrong!!!"
if ( iSULF   == -1) write(*,*) "Warning: SULF  not found in DOS/DEP files, ", &
	"answers may be wrong!!!"

!--- Find the units of this DOSAGE file, and get the conversion factor

Field%timeID   = 1
Field%userTime = TimeOut(1)%time%runTime        ! integration time in hours
avg_time       = TimeOut(1)%time%runTime * 3600 ! integration time in seconds

units = TRIM(ClassChoiceArray(iDOS,iChoice)%units)
call CLOWER( units )
call find_units_conv(units,"ug/m3", avg_time, con_units_conv)
!print*,"units, con_units_conv = ",trim(units),con_units_conv

!--- Find the iClass for the DEPOSITION (drydep+wetdep) file, call it iDEP

iDEP = -1
DO i = 1,nClass
   string1 = ClassStr(i)%string
   CALL CUPPER( string1 )
   IF( INDEX(string1,'SURFACE DEPOSITION') > 0 )THEN
      iDEP = i
      EXIT
   END IF
END DO
IF( iDEP == -1 )THEN
  WRITE(*,*)'SURFACE DOSAGE not available'
  GOTO 9999
END IF

!--- Find the units of this DEPOSITION file, and get the conversion factor

Field%class = iDEP
units = TRIM(ClassChoiceArray(iDep,iChoice)%units)

call CLOWER( units )
call find_units_conv( units, "kg/ha", avg_time, dep_units_conv )
!print*,"units, dep_units_conv = ",trim(units),dep_units_conv

!--- Allocate vars to hold the concentration and deposition, at each Class I
!    Area, for nKind chemical species.  Use MaxRec = max(nRec(:))

ALLOCATE( con(nAreas,MaxRec,nKind), STAT=irv )
IF( irv /= 0 )THEN
   nError = UK_ERROR
   eRoutine = 'sciClass1Flag'
   eMessage = 'Error allocating con'
   WRITE(eInform,'(A,I0,A,I0)')'Request =',MaxRec,' : Error=',irv
   GO TO 9999
END IF
con = 0.

ALLOCATE( dep(nAreas,MaxRec,nKind), STAT=irv )
IF( irv /= 0 )THEN
   nError = UK_ERROR
   eRoutine = 'sciClass1Flag'
   eMessage = 'Error allocating dep'
   WRITE(eInform,'(A,I0,A,I0)')'Request =',MaxRec,' : Error=',irv
   GO TO 9999
END IF
dep = 0.

ALLOCATE( srfFldVal(nAreas,MaxRec,nKind,2), STAT=irv ) ! used temporarily
IF( irv /= 0 )THEN
   nError = UK_ERROR
   eRoutine = 'sciClass1Flag'
   eMessage = 'Error allocating srfFldVal'
   WRITE(eInform,'(A,I0,A,I0)')'Request =',MaxRec,' : Error=',irv
   GO TO 9999
END IF

MaxRank = maxval( output(:)%rank )

ALLOCATE( max_val(nAreas,MaxRec,nOuts,MaxRank), STAT=irv )
IF( irv /= 0 )THEN
   nError = UK_ERROR
   eRoutine = 'sciClass1Flag'
   eMessage = 'Error allocating max_val'
   WRITE(eInform,'(A,I0,A,I0)')'Request =',nAreas*MaxRec*nOuts*MaxRank, &
        ' : Error=',irv
   GO TO 9999
END IF
max_val = -HUGE(0.) ! most negative number

ALLOCATE( max_time(nAreas,MaxRec,nOuts,MaxRank), STAT=irv )
IF( irv /= 0 )THEN
   nError = UK_ERROR
   eRoutine = 'sciClass1Flag'
   eMessage = 'Error allocating max_time'
   WRITE(eInform,'(A,I0,A,I0)')'Request =',nAreas*MaxRec*nOuts*MaxRank, &
        ' : Error=',irv
   GO TO 9999
END IF
max_time = 0

ALLOCATE( iHours(nAreas,nOuts), STAT=irv )
IF( irv /= 0 )THEN
   nError = UK_ERROR
   eRoutine = 'sciClass1Flag'
   eMessage = 'Error allocating iHours'
   WRITE(eInform,'(A,I0,A,I0)')'Request =',nAreas*nOuts, &
        ' : Error=',irv
   GO TO 9999
END IF
iHours = 0

ALLOCATE( sum_val(nAreas,MaxRec,nOuts), STAT=irv )
IF( irv /= 0 )THEN
   nError = UK_ERROR
   eRoutine = 'sciClass1Flag'
   eMessage = 'Error allocating sum_val'
   WRITE(eInform,'(A,I0,A,I0)')'Request =',nAreas*MaxRec*nOuts,' : Error=',irv
   GO TO 9999
END IF
sum_val = 0.

CALL init_error()

irv = SAG_SetSpecialValue( .TRUE.,HP_SPV )

irv = SAG_InitList()

Field%project  = TRIM(Project%name)
Field%path     = TRIM(Project%path)

OPEN(24,FILE='sciClass1Flag.log')

CALL ReadProject( Project )

CALL init_srf_blocks( ntypm )

!--- Dos File

dosFile = TRIM( AddExtension( Project%name,'dos' ) )

!------ Get a new SAG structure

irv = SAG_NewGrdStr( cgrdI )
IF( irv /= SAG_OK )THEN
  nError   = UK_ERROR
  eRoutine = 'InitGrid'
  eMessage = 'Error creating SAG surface grid'
  CALL ReportFileName( eInform,'File=',dosFile )
  GOTO 9999
END IF

!------ Initialize SAG structure

lun_dos = 100 + cgrdI
irv = SAG_InitGridID( dosFile,lun_dos,SAG_GRID_BOTH,0,0,0,cgrdI )
IF( irv /= SAG_OK )THEN
  nError   = UK_ERROR
  eRoutine = 'InitGrid'
  eMessage = 'Error initializing SAG surface grid'
  CALL ReportFileName( eInform,'File=',dosFile )
  GOTO 9999
END IF

!------ Open surface file; read nvart for block version

irv = SAG_OpenID( cgrdI )
IF( irv /= SAG_OK )THEN
  nError   = UK_ERROR
  eRoutine = 'InitGrid'
  eMessage = 'Error opening SAG file'
  CALL ReportFileName( eInform,'File=',dosFile )
  GOTO 9999
END IF

grdC => SAG_PtrGrdStr( cgrdI ) ! Associate "local" grid structure pointer

grdC%mxnam = grdC%nvart

! Deposition file

depFile = TRIM( AddExtension( Project%name,'dep' ) )

!------ Get a new SAG structure

irv = SAG_NewGrdStr( dgrdI )
IF( irv /= SAG_OK )THEN
  nError   = UK_ERROR
  eRoutine = 'InitGrid'
  eMessage = 'Error creating SAG surface grid'
  CALL ReportFileName( eInform,'File=',depFile )
  GOTO 9999
END IF

!------ Initialize SAG structure

lun_dep = 100 + dgrdI
irv = SAG_InitGridID( depFile,lun_dep,SAG_GRID_BOTH,0,0,0,dgrdI )
IF( irv /= SAG_OK )THEN
  nError   = UK_ERROR
  eRoutine = 'InitGrid'
  eMessage = 'Error initializing SAG surface grid'
  CALL ReportFileName( eInform,'File=',depFile )
  GOTO 9999
END IF

!------ Open surface file; read nvart for block version

irv = SAG_OpenID( dgrdI )
IF( irv /= SAG_OK )THEN
  nError   = UK_ERROR
  eRoutine = 'InitGrid'
  eMessage = 'Error opening SAG file'
  CALL ReportFileName( eInform,'File=',depFile )
  GOTO 9999
END IF

grdD => SAG_PtrGrdStr( dgrdI ) ! Associate "local" grid structure pointer

grdD%mxnam = grdD%nvart

!*****************************************************************
!   LOOP OVER TIME STAMPS
!*****************************************************************

do iTime = 1, nTimeOut

! ---Get the concentrations and deposition (wet & dry) at all rececptors,
!    for this time-stamp

   Field%timeID   = iTime
   Field%userTime = TimeOut(iTime)%time%runTime  ! integration time in hours

   iy = TimeOut(iTime)%time%year  ! this year
   im = TimeOut(iTime)%time%month ! this month
   id = TimeOut(iTime)%time%day   ! this month
   ih = TimeOut(iTime)%time%hour  ! this month
   iymdh = iy*1000000 + im*10000 + id*100 + ih

   if ((id == 1 .and. ih == 1) .or. iTime == 1) then
      write(*,*) ! CR
      write(*,'(a,i4,"-",i2.2,a)') "Processing ",iy,im," Day:"
   endif
   if (ih == 1) then
      write(string1,*) id
      write(*,'(a," ",$)') trim(adjustl(string1))
   end if

!--- Retrieve the concentration for this time-stamp

   srfFldVal = 0.
   CALL GetAllFieldValues( cgrdI,nAreas,maxRec,nRec,xRec,yRec, srfFldVal,  &
        .FALSE.,isMult )
   IF( nError /= NO_ERROR )GOTO 9999

   where (srfFldVal == 1.e-30) srfFldVal = 0.
   con = srfFldVal(:,:,:,1)* con_units_conv
   !if( maxval(srfFldVal(:,:,iNO2,1)) > 0. )&
   !write(*,*)ih,maxval(srfFldVal(:,:,iNO2,1)),'(',maxloc(srfFldVal(:,:,iNO2,1)),')'
   !print*,"min con = ",minval(srfFldVal(:,:,:,1))

!--- Retrieve the deposition for this time-stamp

   srfFldVal = 0.
   CALL GetAllFieldValues( dgrdI,nAreas,maxRec,nRec,xRec,yRec, srfFldVal,   &
        .TRUE.,isMult )
   IF( nError /= NO_ERROR )GOTO 9999

   where (srfFldVal == 1.e-30) srfFldVal = 0.
   dep(:,:,:) = (srfFldVal(:,:,:,1) + srfFldVal(:,:,:,2)) * dep_units_conv
!   print*,"max dep = ",maxval(srfFldVal(:,:,:,1))
!   print*,"min dep = ",minval(srfFldVal(:,:,:,1))

!*****************************************************************
!   LOOP OVER CLASS I AREAS
!*****************************************************************

!   print*,"TIME = ",iymdh

   do iOut = 1, nOuts
      do iArea = 1, nAreas

         if (output(iOut)%avgPeriods == 8760 .and. iHours(iArea,iOut) == 0 &
              .and. leap_yr(iy)) output(iOut)%avgPeriods = 8784

         if (output(iOut)%ConDepVis == "CON") then ! iCon is like iKind here

            CALL Class1con( iArea, iCon(iOut), iOut, output(iOut), iymdh)
            IF( nError /= NO_ERROR ) goto 9999

         else if (output(iOut)%ConDepVis == "DEP") then

            CALL Class1dep( iArea, iOut, output(iOut), iymdh)
            IF( nError /= NO_ERROR ) goto 9999

         else if (output(iOut)%ConDepVis == "VIS") then

            CALL Class1vis( iArea, iOut, im, output(iOut), iymdh)
            IF( nError /= NO_ERROR ) goto 9999

         end if

      end do ! do iArea = 1, nAreas
   end do    ! do iOut  = 1, nOuts
end do       ! do iTime = 1, nTimeOut

write(*,*)   ! CR on last write(*,'(1x,i2.2,$)') TimeOut(iTime)%time%day
write(*,*)   ! an extra, just for looks

! Done with calculations, now write out the final answers

do iArea = 1, nAreas
   do iOut = 1, nOuts

!--- Skip the exceedence outputs: we wrote those in class1con() etc.

      if ( output(iOut)%exceeds == HUGE(0.) ) then

!--- find the maximum (over receptors) of the right rank

         ij = -1            ! detect not found
         answer = -HUGE(0.) ! most negative number
         do iRec = 1, nRec(iArea)
            if ( max_val(iArea,iRec,iOut,output(iOut)%rank) > answer) then
               answer = max_val(iArea,iRec,iOut,output(iOut)%rank)
               iymdh  = max_time(iArea,iRec,iOut,output(iOut)%rank)
               ij = iRec
            end if
         end do

!--- Write out the "answer"

         write(string1,*) output(iOut)%rank ! so H1H and H12H both work
         string1 = "H" // trim(adjustl(string1)) // "H"

         if (ij > 0) then
            write(output(iOut)%lun,11) trim(class1Name(iArea)),   &
                 trim(output(iOut)%chemSpecies),                  &
                 lat(iArea,ij),  lon(iArea,ij),                   &
!                 xRec(iArea,ij), yRec(iArea,ij),                  &
                 iymdh,output(iOut)%avgPeriods,trim(string1),     &
                 answer
         else
            write(output(iOut)%lun,11) trim(class1Name(iArea)),   &
                 trim(output(iOut)%chemSpecies),                  &
                 -999., -999.,                                    &
!                 -999., -999.,                                    &
                 -999,output(iOut)%avgPeriods,trim(string1),      &
                 -999.
         end if
11       format(2(a5,","),2(f12.6,","),i11.10,",",i6,",",a,",",g16.6)

         write(*,12) "Wrote ",output(iOut)%ConDepVis," for ",     &
              trim(Class1Name(iArea))," to output",iOut," ",      &
              trim(output(iOut)%outFile)
12       format(5a,i3,2a)

      end if ! if ( output(iOut)%exceeds < HUGE(0.) ) then
   end do    ! do iOut  = 1, nOuts
end do       ! do iArea = 1, nAreas

!--- Close the output files

do iOut = 1, nOuts
   inquire( output(iOut)%lun, opened=lexist)
   if ( lexist ) close( output(iOut)%lun )
end do ! do iOut  = 1, nOuts

!--- WE'RE DONE!!!

9999 CONTINUE

IF( nError /= NO_ERROR )THEN
  WRITE(*,*)TRIM(eMessage)
  WRITE(*,*)TRIM(eInform)
  WRITE(*,*)TRIM(eAction)
END IF

IF( ALLOCATED(output)         ) DEALLOCATE(output)
IF( ALLOCATED(con)            ) DEALLOCATE(con)
IF( ALLOCATED(dep)            ) DEALLOCATE(dep)
IF( ALLOCATED(srfFldVal)      ) DEALLOCATE(srfFldVal)
IF( ALLOCATED(lat)            ) DEALLOCATE(lat,lon)
IF( ALLOCATED(xRec)           ) DEALLOCATE(xRec,yRec)
IF( ALLOCATED(nRec)           ) DEALLOCATE(nRec)
IF( ALLOCATED(max_val)        ) DEALLOCATE(max_val)
IF( ALLOCATED(max_time)       ) DEALLOCATE(max_time)
IF( ALLOCATED(iHours)         ) DEALLOCATE(iHours)
IF( ALLOCATED(sum_val)        ) DEALLOCATE(sum_val)
IF( ALLOCATED(NatCond)        ) DEALLOCATE(NatCond)
IF( ALLOCATED(LargeRH)        ) DEALLOCATE(LargeRH)
IF( ALLOCATED(SmallRH)        ) DEALLOCATE(SmallRH)
IF( ALLOCATED(SeaSalt)        ) DEALLOCATE(SeaSalt)
IF( ALLOCATED(Bext_nat)       ) DEALLOCATE(Bext_nat)
IF( ALLOCATED(iCon)           ) DEALLOCATE(iCon)
IF( ALLOCATED(class1Name)     ) DEALLOCATE(class1Name)
IF( ALLOCATED(class1Filename) ) DEALLOCATE(class1Filename)

!------ Close file and deallocate grid structure

INQUIRE(UNIT=24,OPENED=lopen,IOSTAT=ios)
IF( lopen )CLOSE(24)

irv = SAG_CloseID( cgrdI )
irv = SAG_RmvGrdStr( cgrdI )
irv = SAG_CloseID( dgrdI )
irv = SAG_RmvGrdStr( dgrdI )

!-- Release memory
CALL ClearMemory()

!-- Clean up : Exit SCIPtool
IF ( lInit ) irv = SCIPExitTool()

END PROGRAM runSciClass1Flag
