!*******************************************************************************
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                        :: iHours     ! running total
     integer                        :: rank       ! 1st or 4th or 8th high...
     real                           :: exceeds    ! print all vals over this
     real                           :: answer     ! highest value over avg per
  end TYPE an_output

  REAL, DIMENSION(:,:,:),  ALLOCATABLE :: con       ! index = iArea,iRec,ikind
  REAL, DIMENSION(:,:,:),  ALLOCATABLE :: ddep,wdep ! 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 :: NatCond, LargeRH, SmallRH, SeaSalt
  REAL, DIMENSION(:,:),    ALLOCATABLE :: Bext_nat  ! index = iArea, Month

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

  INTEGER, DIMENSION(:),   ALLOCATABLE :: iCon
  INTEGER, DIMENSION(:),   ALLOCATABLE :: nRec      ! index = iArea

! 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 :: iANO3I, iANO3J, iANO3K, iNH3, iANH4I, iANH4J, iANH4K, iSO2
  INTEGER :: iSULF, iASO4I, iASO4J, iASO4K

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

  ! Already declared iASO4I, iASO4J, iASO4K iANO3I, iANO3J, iANO3K above
  INTEGER :: iPSOA, iPOC, iAECI, iAECJ, iASOIL, iACORS, iANAK, iACLK

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

  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 = (/ "-", "_", ":" /)

END MODULE class1Tools

!
!*************************************************************************
!
  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 == "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 InitClass1answers()

! 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 InitClass1answers

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

  SUBROUTINE InitToolClass1answers()

    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 = '/home/user/sid/scipuff/SCIP_Working/bin/linux/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

    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,*)
    WRITE(6,*)'Initializing SCIPtool from '//TRIM(ini_file)
    WRITE(6,*)
    WRITE(6,'("Running Class1answers 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 InitToolClass1answers
!
!*************************************************************************
!
PROGRAM srf2smp

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

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, samHdr
CHARACTER (len=256)              :: line
CHARACTER (len=PATH_MAXLENGTH)   :: dosFile, depFile, dsiFile

INTEGER, EXTERNAL   :: iargc            ! to retrieve the command line
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             :: ios,irv          ! return value for SCIP calls
INTEGER             :: ij               ! index of lat,lon where answer is found
INTEGER             :: yr,mo            ! year and month
INTEGER             :: iClass, iChoice, iKind, iDOS, iDEP
INTEGER             :: i, j, iTime
INTEGER             :: cGrdI, dGrdI
INTEGER             :: lun_dos,lun_dep,lun_sam

LOGICAL             :: lexist           ! does a file exist?
LOGICAL             :: isMULT = .false. ! multi-component or tracer run?
LOGICAL             :: isLL = .false.   ! is LL

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

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

CHARACTER(128), EXTERNAL :: AddExtension

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


!--- Start up : Read Command line input

dsiFile  = ''
projName = ''

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

if (iargc() > 0) THEN
  call getarg(1,projName)
  if (iargc() > 1) THEN
    call getarg(2,dsiFile)
  end if
end if

!projName = "tva_980825_1yr_pr_nosam"
!fileName = "fc_egu1_07s3ab"

!dsiFile  = "tva_980825_sfc.sam"
!dsiFile  = "../fc_camxreceptors4.sam"

IF( LEN_TRIM(projName) == 0 )THEN
 WRITE(*,*)'Project Name ?'
 READ(*,*)projName
END IF

IF( LEN_TRIM(dsiFile) == 0 )THEN
  WRITE(*,*)'srf2smp location input file name ?'
  READ(*,*)dsiFile
END IF

ALLOCATE( nRec(1),STAT=astat )

nRec    = 0
lun_sam = 11

OPEN(lun_sam,FILE=TRIM(dsiFile),IOSTAT=ios)
IF( ios /= 0 )THEN
  nError = UK_ERROR
  WRITE(eMessage,*) 'Error opening sampler location file with ios = ',ios
  GO TO 9999
END IF

READ(lun_sam,'(A)',IOSTAT=ios)samHdr
IF( ios /= 0 )THEN
  nError = UK_ERROR
  WRITE(eMessage,*) 'Error reading sampler location file header'
  GO TO 9999
END IF

CALL CUPPER(samHdr)

isLL   = .FALSE.
nAreas = 1

IF( INDEX(samHdr,'LATLON') > 0 )isLL = .TRUE.

DO
  READ(lun_sam,*,IOSTAT=ios)x,y
  IF( ios < 0 )THEN
    EXIT
  ELSE IF( ios > 0 )THEN
    nError = UK_ERROR
    WRITE(eMessage,*) 'Error reading sampler locations'
    GO TO 9999
  END IF
  nRec(1) = nRec(1) +1
END DO

IF( ALLOCATED(xRec) )DEALLOCATE(xRec,yRec,STAT=astat)

maxRec = nRec(1)

ALLOCATE( xRec(nAreas,maxRec),yRec(nAreas,maxRec),STAT=astat )
IF( ios > 0 )THEN
  nError = UK_ERROR
  WRITE(eMessage,*) 'Error reading sampler locations'
  GO TO 9999
END IF

REWIND(lun_sam,IOSTAT=ios)
IF (ios /= 0) THEN
  nError = UK_ERROR
  WRITE(eMessage,*) 'Error rewinding sampler locations file'
  GO TO 9999
END IF

READ(lun_sam,*,IOSTAT=ios)samHdr
IF( ios /= 0 )THEN
  nError = UK_ERROR
  WRITE(eMessage,*) 'Error reading sampler location file header'
  GO TO 9999
END IF

DO i = 1,nAreas*maxRec
  READ( lun_sam,'(A)',IOSTAT=ios )line
  READ( line,*,IOSTAT=ios)xRec(1,i),yRec(1,i)
  IF( ios /= 0 )THEN
    READ(line,*,IOSTAT=ios)xRec(1,i),yRec(1,i)
    IF (ios /= 0) THEN
      nError = UK_ERROR
      WRITE(eMessage,*) 'Error reading sampler locations'
      WRITE(eAction,*)'Check line ',TRIM(line)
      GO TO 9999
    END IF
  END IF
END DO

CLOSE(unit=lun_sam,IOSTAT=ios)

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

!--- Start up : Initialize common

CALL InitClass1answers()

!--- Start up : Write Screen Header

WRITE(maxPath,'(I0)') PATH_MAXLENGTH

WRITE(6,'(/,A)')'Program : Class1answers 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 InitToolClass1answers()
IF( nError /= NO_ERROR )GOTO 9999


! -- Get Plot Times

CALL SplitName( ProjName,Project%name,Project%path ) !project name & path
Project%ID      = 0
Project%version = 0

!==== Action : Retrieve output times

WRITE(*,*)'Call SCIPNumPlotTimes'

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
WRITE(*,*)'After call to SCIPNumPlotTimes'

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
IF( isLL )THEN
  ALLOCATE( lon(maxRec),lat(maxRec),STAT=irv )
  DO j = 1,nAreas
    lon = xRec(j,:)
    lat = yRec(j,:)
    WRITE(*,*)
    irv = SCIPGetPrjCoord( callerID,maxRec,lon,lat,xRec(j,:),yRec(j,:) )
    IF( irv /= SCIPsuccess )THEN
      CALL toolError( 'Failed converting to project coordinates' )
      GOTO 9999
    END IF
  END DO
END IF

!==== 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
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'

do i = 1, nKind

   if (isMULT) then                        ! multi-species output
      Field%kind = i + 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 )

end do  ! do i = 1, nKind

!--- 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)

!--- 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
Field%class = iDEP

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

do i = 1, nKind

   if (isMULT) then                        ! multi-species output
      Field%kind = i + ClassChoiceArray(iDep,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
   j = INDEX(string1,'(')
   IF( j > 0 )string1 = string1(1:j-1)
   CALL CUPPER( string1 )

! These concentrations are needed for deposition

   if (TRIM(string1) == 'NO'   ) iNO    = i
   if (TRIM(string1) == 'NO2'  ) iNO2   = i
   if (TRIM(string1) == 'NO3'  ) iNO3   = i
   if (TRIM(string1) == 'N2O5' ) iN2O5  = i
   if (TRIM(string1) == 'HNO3' ) iHNO3  = i
   if (TRIM(string1) == 'HONO' ) iHONO  = i
   if (TRIM(string1) == 'PNA'  ) iPNA   = i
   if (TRIM(string1) == 'PAN'  ) iPAN   = i
   if (TRIM(string1) == 'PANX' ) iPANX  = i
   if (TRIM(string1) == 'NTR'  ) iNTR   = i
   if (TRIM(string1) == 'ANO3I') iANO3I = i
   if (TRIM(string1) == 'ANO3J') iANO3J = i
   if (TRIM(string1) == 'ANO3K') iANO3K = i
   if (TRIM(string1) == 'NH3'  ) iNH3   = i
   if (TRIM(string1) == 'ANH4I') iANH4I = i
   if (TRIM(string1) == 'ANH4J') iANH4J = i
   if (TRIM(string1) == 'ANH4K') iANH4K = i
   if (TRIM(string1) == 'SO2'  ) iSO2   = i
   if (TRIM(string1) == 'SULF' ) iSULF  = i
   if (TRIM(string1) == 'ASO4I') iASO4I = i
   if (TRIM(string1) == 'ASO4J') iASO4J = i
   if (TRIM(string1) == 'ASO4K') iASO4K = i

end do  ! do i = 1, nKind

!--- Find the units of this DEPOSITION 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(iDep,iChoice)%units)

call CLOWER( units )
call find_units_conv( units, "kg/ha", avg_time, 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 = 'class1answers'
   eMessage = 'Error allocating con'
   WRITE(eInform,'(A,I0,A,I0)')'Request =',MaxRec,' : Error=',irv
   GO TO 9999
END IF
con = 0.

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

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

CALL init_error()

irv = SAG_SetSpecialValue( .TRUE.,HP_SPV )

irv = SAG_InitList()

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


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

   yr = TimeOut(iTime)%time%year  ! this year
   mo = TimeOut(iTime)%time%month ! this month

!--- 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

   con = 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

   ddep(:,:,:) = srfFldVal(:,:,:,1)
   wdep(:,:,:) = srfFldVal(:,:,:,2)

!*****************************************************************
!   LOOP OVER CLASS I AREAS
!*****************************************************************
   WRITE(99,*)iTime

   DO i = 1,nKind

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

     WRITE(99,*)i,TRIM(string1),MAXVAL(con(1,:,i)),MAXVAL(ddep(1,:,i)),MAXVAL(wdep(1,:,i))

   END DO

   WRITE(99,*)

end do       ! do iTime = 1, nTimeOut

!--- DONE!

goto 9999

99 nError   = UK_ERROR
eRoutine = 'Class1answers'
eMessage = 'Error opening output files'

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(ddep)           ) DEALLOCATE(ddep,wdep)
IF( ALLOCATED(srfFldVal)      ) DEALLOCATE(srfFldVal)
IF( ALLOCATED(lat)            ) DEALLOCATE(lat,lon)
IF( ALLOCATED(xRec)           ) DEALLOCATE(xRec,yRec)
IF( ALLOCATED(nRec)           ) DEALLOCATE(nRec)


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

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()

STOP

END PROGRAM srf2smp
