c**** READPTSA
c
      subroutine readptsa(idate,btim)
      use filunit
      use grid
      use camxcom
      use ptemiss
      use tracer
c
c----CAMx v7Beta6 190902
c
c-----------------------------------------------------------------------
c    Description:
c-----------------------------------------------------------------------
c
c   This routine reads one hour of emissions from each of the elevated
c   emissions files and calculates the emissions levels.  It then
c   places these emissions in the appropriate place in the point source
c   array used for tracer emissions.  The tracer to which the emissions
c   is assigned depends on the source group and the source region.
c
c     Copyright 1996 - 2018
c     Ramboll
c
c     Inputs:
c          idate  I  date of current hour (YYJJJ)
c          btim   R  time of current hour
c
c-----------------------------------------------------------------------
c    LOG:
c-----------------------------------------------------------------------
c
c     06/06/96   --gwilson--    Original development
c     12/10/96   --gwilson--    Added switch to skip emissions if PiG
c     02/01/97   --gwilson--    Put fuzz factor of 0.1 ppb to catch if
c                               emissions of leftover group is negative.
c     02/03/97   --gwilson--    Put code to ignore emissions on boundary.
c     01/20/99   --cemery---    Grid cell size from file should be meters
c                               for all cartesian projections (UTM, LCP,
c                               PSP)
c     11/06/01   --cemery--     Input dates are now Julian
c     07/19/02   --gwilson--    Added seperate source area map for each grids.
c     03/01/16   --gwilson--    Added partial source area map
c
c-----------------------------------------------------------------------
c    Include files:
c-----------------------------------------------------------------------
c
      include 'camx.prm'
      include 'flags.inc'
c
c-----------------------------------------------------------------------
c    Arguement declarations:
c-----------------------------------------------------------------------
c
      integer   idate
      real      btim
c
c-----------------------------------------------------------------------
c    Local parameters:
c-----------------------------------------------------------------------
c
c   FUZZ   R    fuzz factor used to determine if emissions are truly < 0
c
      real   FUZZ
c
      parameter( FUZZ = -0.0001 )
c
c-----------------------------------------------------------------------
c    Local variables:
c-----------------------------------------------------------------------
c
      character*200 fname
      integer imap, i, j, ii, jj, ipart, idxbase
      integer ndate, icel, jcel, num_ptsfiles
      real    xloctmp, yloctmp, ttime, frac, xloccrs, yloccrs
      logical is_netcdf_file
c
      integer, allocatable, dimension(:)   :: izcel
      real,    allocatable, dimension(:,:) :: lefcls
      real,    allocatable, dimension(:,:) :: emscls
c
c-----------------------------------------------------------------------
c    Entry point:
c-----------------------------------------------------------------------
c
c   ---- allocate the large local arrays ---
c
      allocate( izcel(MXPTSRC) )
      izcel = 0
      allocate( lefcls(MXTRCLS,MXPTSRC) )
      allocate( emscls(MXTRCLS,MXPTSRC) )
c
c   --- set the date and times ---
c
      ndate = idate
      ttime = btim/100.0
c
c   --- initialize emissions to zero ---
c
      do i=1,nptsrc
         do j=1,nsaspc
            sapnts(i,j) = 0.
         enddo
         do j=1,ntrcls
            emscls(j,i) = 0. 
         enddo
      enddo
c
c  --- loop over groups ---
c
      do igroup = 0,ngroup
         num_ptsfiles = num_iortpt(igroup)
         if( igroup .EQ. 0 ) num_ptsfiles = npoint_files
         do idxfile=1,num_ptsfiles
c
c   --- set the file name for elevated points emissions file ---
c
            is_netcdf_file = .FALSE. 
            if( igroup .EQ. 0 ) then
                write(fname,'(A,I3)') 'PTSOURCE -- UNIT ',iptem(idxfile)
                if( is_netcdf_iptem(idxfile) ) is_netcdf_file = .TRUE.
            else
                fname = tptfil(igroup,idxfile)
                if( is_netcdf_iortpt(igroup,idxfile) ) is_netcdf_file = .TRUE.
            endif
c
c  --- for each group, call routine to read the emissions for
c      the group ----
c
            if( is_netcdf_file ) then
               call ncf_rdptgrp(igroup,ntrcls,emscls,izcel)
            else
               call rdptgrp(ndate,ttime,igroup,ntrcls,emscls,izcel)
            endif
c
c  --- if doing source groups calculate the 
c      left-over and put emissions into proper place in the emissions
c      array ---
c
            if( ngroup .GT. 0 ) then
               do 10 i=1,nptsrc
c
c  --- find index into the master list ---
c
                  idxbase = idx_point_in_list(igroup,idxfile,i)
                  if( idxbase .LT. 0 ) goto 10
c
c  --- get the region for this cell from mapping array ----
c
                  if( .NOT.llatlon ) then
                     icel = FLOOR( (xlocpt(idxbase)/1000. - xorg) / delx ) + 1
                     jcel = FLOOR( (ylocpt(idxbase)/1000. - yorg) / dely ) + 1
                     xloccrs = xlocpt(idxbase)/1000. - xorg
                     yloccrs = ylocpt(idxbase)/1000. - yorg
                  else
                     icel = FLOOR( (xlocpt(idxbase) - xorg) / delx ) + 1
                     jcel = FLOOR( (ylocpt(idxbase) - yorg) / dely ) + 1
                     xloccrs = xlocpt(idxbase) - xorg
                     yloccrs = ylocpt(idxbase) - yorg
                  endif
                  if(icel .LE. 0 .OR. icel .GT. ncol(1)) goto 10
                  if(jcel .LE. 0 .OR. jcel .GT. nrow(1)) goto 10
c
c   --- find out if a nest contains this source  ---
c
                  igrd = 1
                  do ig = 2,ngrid
                    xloctmp = xloccrs - (inst1(ig)-1)*delx
                       yloctmp = yloccrs - (jnst1(ig)-1)*dely
                    if( xloctmp .GT. 0. .AND. yloctmp .GT. 0. ) then
                       ii = 2 + FLOOR( xloctmp/delx * FLOAT( meshold(ig) ))
                       jj = 2 + FLOOR( yloctmp/dely * FLOAT( meshold(ig) ))
                       if( ii .GT. 1 .AND. jj .GT. 1 .AND. ii .LT. ncol(ig)
     &                                    .AND. jj .LT. nrow(ig) ) then
                          igrd = ig
                          icel = ii
                          jcel = jj
                       endif
                    endif
                  enddo
                  imap = igrmap(0,1,igrd,icel,jcel)
                  frac = frcmap(0,1,igrd,icel,jcel)
                  if( izcel(idxbase) .LT. 0 ) then
                     if( .NOT. lptoverride ) then
                        write(iout,'(//,a)') 'ERROR in READPTSA:'
                        write(iout,'(/,1X,2A )') 'A source was read that has the point ',
     &                      'source override trigger turned on.'
                        write(iout,'(1X,2A)') 'If you intend to use point source ',
     &                     'override please set the namelist variable'
                        write(iout,'(1X,2A)') 'SA_PT_Override to TRUE in the CAMx_Control namelist.'
                        write(iout,'(1X,2A)') 'If you do not want point source override ',
     &                         'you need to set the kcell variable '
                        write(iout,'(1X,2A)') 'in your point source file to zero.'
                        call camxerr()
                     endif
                     imap = ABS( izcel(idxbase) )
                  endif
                  if( (imap .LE. 0 .AND. frac .GT. 0.) .OR. imap .GT. nregin ) then
                     write(iout,'(//,a)') 'ERROR in READPTSA:'
                     write(iout,'(/,1X,3A)') 'Invalid region found in',
     &                           ' point source override when reading ',
     &                                              'point source file.'
                     write(iout,'(1X,A,I4)') 'Region code      : ',imap
                     write(iout,'(1X,A,I4)') 'Number of regions: ',nregin
                     write(iout,'(10X,A,/,A)') 'Point source filename: ',
     &                                             fname(:istrln(fname))
                     write(iout,'(1X,2A)') 'Check the values in the ',
     &                                           'point source overide.'
                     call camxerr()
                  endif
c
c  --- calculate the leftover emissions to use in last source group ---
c
                  do icls=1,ntrcls
                     if( igroup .EQ. 0 ) then
                         lefcls(icls,idxbase) = emscls(icls,idxbase)
                     else
                         lefcls(icls,idxbase) = 
     &                         lefcls(icls,idxbase) - emscls(icls,idxbase)
c
c  --- put emissions into position in gridded tracer emissions array ---
c
                         ipos = iemcls(icls) - 1 + imap + (igroup-1)*nregin
                         sapnts(idxbase,ipos) = emscls(icls,idxbase) * frac
                     endif
                  enddo
c
c  --- next point source ----
c
   10         continue
c
c  --- only one group, just load emissions into arrays from
c      position 0 in gridded array ----
c
            else
               do 20 i=1,nptsrc
c
c  --- find index into the master list ---
c
                  idxbase = idx_point_in_list(igroup,idxfile,i)
                  if( idxbase .LT. 0 ) goto 20
c
c  --- get the region for this cell from mapping array ----
c
                  if( .NOT.llatlon ) then
                     icel = FLOOR( (xlocpt(idxbase)/1000. - xorg) / delx ) + 1
                     jcel = FLOOR( (ylocpt(idxbase)/1000. - yorg) / dely ) + 1
                     xloctmp = xlocpt(idxbase)/1000. - xorg
                     yloctmp = ylocpt(idxbase)/1000. - yorg
                  else
                     icel = FLOOR( (xlocpt(idxbase) - xorg) / delx ) + 1
                     jcel = FLOOR( (ylocpt(idxbase) - yorg) / dely ) + 1
                     xloctmp = xlocpt(idxbase) - xorg
                     yloctmp = ylocpt(idxbase) - yorg
                  endif
                  if(icel .LE. 0 .OR. icel .GT. ncol(1)) goto 20
                  if(jcel .LE. 0 .OR. jcel .GT. nrow(1)) goto 20
c
c   --- find out if a nest contains this source  ---
c
                  igrd = 1
                  do ig = 2,ngrid
                    xloctmp = xloctmp - (inst1(ig)-1)*delx
                    yloctmp = yloctmp - (jnst1(ig)-1)*dely
                    ii = 2 + FLOOR( xloctmp/delx * FLOAT( meshold(ig) ) )
                    jj = 2 + FLOOR( yloctmp/dely * FLOAT( meshold(ig) ) )
                    if( ii .GT. 1 .AND. jj .GT. 1 .AND. 
     &                      ii .LT. ncol(ig) .AND. jj .LT. nrow(ig) ) then
                       igrd = ig
                       icel = ii
                       jcel = jj
                     endif
                  enddo
                  imap = igrmap(0,1,igrd,icel,jcel)
                  frac = frcmap(0,1,igrd,icel,jcel)
                  if( izcel(idxbase) .LT. 0 ) then
                     if( .NOT. lptoverride ) then
                        write(iout,'(//,a)') 'ERROR in READPTSA:'
                        write(iout,'(/,1X,2A )') 'A source was read that has the point ',
     &                   'source override trigger turned on.'
                        write(iout,'(1X,2A)') 'If you intend to use point source ',
     &                     'override please set the namelist variable'
                        write(iout,'(1X,2A)') 'SA_PT_Override to TRUE in the CAMx_Control namelist.'
                        write(iout,'(1X,2A)') 'If you do not want point source override ',
     &                         'you need to set the kcell variable '
                        write(iout,'(1X,2A)') 'in your point source file to zero.'
                        call camxerr()
                     endif
                     imap = ABS( izcel(idxbase) )
                  endif
                  if( (imap .LE. 0 .AND. frac .GT. 0) .OR. imap .GT. nregin ) then
                     write(iout,'(//,a)') 'ERROR in READPTSA:'
                     write(iout,'(/,1X,A,A,I4)') 'Invalid region found in',
     &                           ' point source override when reading ',
     &                                              'point source file.'
                     write(iout,'(1X,A,I4)') 'Region code      : ',imap
                     write(iout,'(1X,A,I4)') 'Number of regions: ',nregin
                     write(iout,'(10X,A,/,A)') 'Point source filename: ',
     &                                             fname(:istrln(fname))
                     write(iout,'(1X,2A)') 'Check the values in the ',
     &                                           'point source overide.'
                     call camxerr()
                  endif
c
c   --- put emissions in array at correct offset ---
c
                  do icls=1,ntrcls
                     ipos = iemcls(icls) - 1 + imap
                     sapnts(idxbase,ipos) = emscls(icls,idxbase) * frac
                  enddo
c
c  --- next point source ----
c
   20          continue
            endif
c
c  --- get emissions file from the next group ----
c
          enddo
      enddo
c
c   --- all groups read, load the leftover group ---
c
      if( ngroup .GT. 0 .AND. leftovr ) then
         do 30 i=1,nptsrc
c
c  --- find index into the master list ---
c
           idxbase = idx_point_in_list(igroup,idxfile,i)
           if( idxbase .LT. 0 ) goto 30
c
c
c  --- get the region for this cell from mapping array ----
c
            if( .NOT.llatlon ) then
               icel = FLOOR( (xlocpt(idxbase)/1000. - xorg) / delx ) + 1
               jcel = FLOOR( (ylocpt(idxbase)/1000. - yorg) / dely ) + 1
               xloccrs = xlocpt(idxbase)/1000. - xorg
               yloccrs = ylocpt(idxbase)/1000. - yorg
            else
               icel = FLOOR( (xlocpt(idxbase) - xorg) / delx ) + 1
               jcel = FLOOR( (ylocpt(idxbase) - yorg) / dely ) + 1
               xloccrs = xlocpt(idxbase) - xorg
               yloccrs = ylocpt(idxbase) - yorg
            endif
            if(icel .LE. 0 .OR. icel .GT. ncol(1)) goto 30
            if(jcel .LE. 0 .OR. jcel .GT. nrow(1)) goto 30
c
c   --- find out if a nest contains this source  ---
c
            igrd = 1
            do ig = 2,ngrid
              xloctmp = xloccrs - (inst1(ig)-1)*delx
              yloctmp = yloccrs - (jnst1(ig)-1)*dely
              ii = 2 + FLOOR( xloctmp/delx * FLOAT( meshold(ig) ) )
              jj = 2 + FLOOR( yloctmp/dely * FLOAT( meshold(ig) ) )
              if( ii .GT. 1 .AND. jj .GT. 1 .AND.
     &                ii .LT. ncol(ig) .AND. jj .LT. nrow(ig) ) then
                 igrd = ig
                 icel = ii
                 jcel = jj
               endif
            enddo
            imap = igrmap(0,1,igrd,icel,jcel)
            frac = frcmap(0,1,igrd,icel,jcel)
            if( izcel(idxbase) .LT. 0 ) then
                if( .NOT. lptoverride ) then
                  write(iout,'(//,a)') 'ERROR in READPTSA:'
                  write(iout,'(/,1X,2A )') 'A source was read that has the point ',
     &                   'source override trigger turned on.'
                  write(iout,'(1X,2A)') 'If you intend to use point source ',
     &                     'override please set the namelist variable'
                  write(iout,'(1X,2A)') 'SA_PT_Override to TRUE in the CAMx_Control namelist.'
                  write(iout,'(1X,2A)') 'If you do not want point source override ',
     &                         'you need to set the kcell variable '
                  write(iout,'(1X,2A)') 'in your point source file to zero.'
                  call camxerr()
                endif
                imap = ABS( izcel(idxbase) )
            endif
            if( (imap .LE. 0 .AND. frac .GT. 0.) .OR. imap .GT. nregin ) then
               write(iout,'(//,a)') 'ERROR in READPTSA:'
               write(iout,'(/,1X,A,A,I4)') 'Invalid region found in',
     &          ' point source override when reading point source file.'
               write(iout,'(1X,A,I4)') 'Region code      : ',imap
               write(iout,'(1X,A,I4)') 'Number of regions: ',nregin
               write(iout,'(10X,A,/,A)') 'Point source filename: ',
     &                                             fname(:istrln(fname))
               write(iout,'(1X,2A)') 'Check the values in the point ',
     &                                                 'source overide.'
               call camxerr()
            endif
c
c   ---- make sure leftover group is not negative ---
c
            do icls=1,ntrcls
               if( lefcls(icls,idxbase) .LT. 0. ) then
                  if( lefcls(icls,idxbase) .GT. FUZZ ) then
                    lefcls(icls,idxbase) = 0.
                  else
                    write(iout,'(//,A)') 'ERROR in READPTSA:'
                    write(iout,'(/,4A,I5,A)')
     &                     'Negative ',clsnam(idxcls(icls)),
     &                       ' emissions calculated in leftover group ',
     &                                        'in point source: ',idxbase,' :'
                    write(iout,'(A,F20.6)')'   Value = ',lefcls(icls,idxbase)
                    call camxerr()
                  endif
               endif
               ipos = iemcls(icls) - 1 + imap + ngroup * nregin
               sapnts(idxbase,ipos) = lefcls(icls,idxbase) * frac
            enddo
   30    continue
      endif
c
c   ---- deallocate the large local arrays ---
c
      deallocate( izcel  )
      deallocate( lefcls )
      deallocate( emscls )
c
c  --- all files read ---
c
      return
      end
