C**** NCF_READAR
c
      subroutine ncf_readar(igrd,ifile,iunit,ncol,nrow,nlay,nlay_ems,
     &                                                aremis,numspcs)
      use chmstry
      use camxcom
      use filunit
      implicit none
c
c----CAMx v7Beta6 190902
c
c-----------------------------------------------------------------------
c    Description:
c-----------------------------------------------------------------------
c
c     NCF_READAR reads the NetCDF gridded emissions file for 
c     the given grid and loads the data into global arrays for emissions
c     injection.
c
c     Copyright 1996 - 2018
c     Ramboll
c
c     Input arguments:
c        igrd             grid index
c        ifile            index of file in list for this grid
c        ncol             number of columns
c        nrow             number of rows
c        nlay             number of layers in grid
c        nlay_ems         number of layers in emissions foles
c        iunit            file unit for area emissions file
c        numspcs          number of model species
c
c     Output arguments:
c        aremis              area emissions rate (mole/s)
c
c
c-----------------------------------------------------------------------
c    LOG:
c-----------------------------------------------------------------------
c
c     02/20/17   --gwilson--    Original development
c
c-----------------------------------------------------------------------
c    Include files:
c-----------------------------------------------------------------------
c
      include 'camx.prm'
      include 'flags.inc'
      include 'vbs.inc'
      include 'netcdf.inc'
c
c-----------------------------------------------------------------------
c    Argument declarations:
c-----------------------------------------------------------------------
c
      integer igrd
      integer ifile 
      integer iunit
      integer ncol
      integer nrow
      integer nlay
      integer nlay_ems
      integer numspcs
      real    aremis(ncol,nrow,nlay_ems,numspcs)
c
c-----------------------------------------------------------------------
c    External functions:
c-----------------------------------------------------------------------
c
      integer istrln
      integer ncf_get_tstep
c
c-----------------------------------------------------------------------
c    Local variables:
c-----------------------------------------------------------------------
c
      character*200 action
      character*10  this_var
      integer       this_tstep, ierr, nlays_in, ispc, i, j, k, l
      integer       this_varid, data_start(4), data_count(4)
      integer       this_time_tflag, this_time_etflag, this_dimid
      real          poa_gv_em(ncol,nrow)
      real          poa_dv_em(ncol,nrow)
      real          poa_mc_em(ncol,nrow)
      real          poa_op_em(ncol,nrow)
      real          poa_bb_em(ncol,nrow)
c
      real, allocatable, dimension(:,:,:) ::  argrid
c
c-----------------------------------------------------------------------
c    Entry point:
c-----------------------------------------------------------------------
c
      if( .NOT. is_netcdf_iarem(igrd,ifile) ) goto 9999
      write(action,'(2A,I3,A,I3)') 'Reading the NetCDF ',
     &                              'gridded emissions file. Grid: ',
     &                                           igrd,' File: ',ifile
c
c  ---- get the index for timestep containing this time ----
c
      this_tstep = ncf_get_tstep(iunit,action,date,time,
     &              this_time_tflag,this_time_etflag,le1day,.TRUE.)
c
c   ---- get the number of layers in this file ---
c
c  --- get the id for the layer dimension ---
c
      ierr = nf_inq_dimid(iunit, "LAY", this_dimid  )
      if( ierr .NE. NF_NOERR ) goto 7000
c
      ierr = nf_inq_dimlen(iunit,this_dimid,nlays_in)
      if( ierr .NE. NF_NOERR ) goto 7000
c
c  ---- allocate the temporary array ---
c
      allocate( argrid(ncol,nrow,nlays_in) )
c
c   --- if number of layers is not 1 make sure it matches grid ---
c
      if( nlays_in .NE. 1 .AND. nlays_in .NE. nlay ) goto 7002
c
c   --- set the indexes for what to read ---
c
      data_start(1) = 1
      data_count(1) = ncol
      data_start(2) = 1
      data_count(2) = nrow
      data_start(3) = 1
      data_count(3) = nlays_in
      data_start(4) = this_tstep
      data_count(4) = 1
c
c  ---- loop over the list of emissions spcies ---
c
      do ispc=1,nemspc
c
c  --- load data into temporary arrray ---
c
          this_var = emspcname(ispc)
          ierr = nf_inq_varid(iunit, this_var, this_varid)
          if( ierr .NE. NF_NOERR) cycle
          ierr = nf_get_vara_real(iunit,this_varid,data_start,
     &                                               data_count,argrid)
          if( ierr .NE. NF_NOERR) goto 7001
c
c-----Put the data into the global array ----
c
          do k=1,nlays_in
             do j=1,nrow
                do i=1,ncol
                   if (argrid(i,j,k).lt.0.) then
                     write(iout,'(//,a)') 'ERROR in NCF_READAR:'
                     write(iout,'(A)') action(:istrln(action))
                     write(iout,'(a)') 'Negative emissions read.'
                     write(iout,'(a,3i3,a,i3)') 'Cell (i,j,k): ',i,j,k,
     &                                                      ' Species: ',ispc
                     call camxerr()
                   endif
                   aremis(i,j,k,ispc) = aremis(i,j,k,ispc) + 
     &                                      argrid(i,j,k)/(60.*dtems) 
                enddo
             enddo
          enddo
      enddo
c
cgwilson          if(lvbs .and. LVBSPREPROC) then
cgwilson            if( lemmap(kpap_c(0)) .EQ. ll ) then ! POA_OP in VBS emiss
cgwilson              poa_op_em(:,:) = argrid(:ncol,:nrow)
cgwilson              goto 50
cgwilson            endif
cgwilson            if( lemmap(kpap_c(1)) .EQ. ll ) then ! POA_GV in VBS emiss
cgwilson              poa_gv_em(:,:) = argrid(:ncol,:nrow)
cgwilson              goto 50
cgwilson            endif
cgwilson            if( lemmap(kpap_c(2)) .EQ. ll ) then ! POA_DV in VBS emiss
cgwilson              poa_dv_em(:,:) = argrid(:ncol,:nrow)
cgwilson              goto 50
cgwilson            endif
cgwilson            if( lemmap(kpcp_c(0)) .EQ. ll ) then ! POA_MC in VBS emiss
cgwilson              poa_mc_em(:,:) = argrid(:ncol,:nrow)
cgwilson              goto 50
cgwilson            endif
cgwilson            if( lemmap(kpfp_c(0)) .EQ. ll ) then ! POA_BB in VBS emiss
cgwilson              poa_bb_em(:,:) = argrid(:ncol,:nrow)
cgwilson              goto 50
cgwilson            endif
cgwilson        endif

c
cgwilson 50   continue
cgwilson      iF (LVBS .and. LVBSPREPROC) then
cgwilson        do l = 0, NVOLBIN
cgwilson          aremis(:,:,kpap_c(l)) = poa_op_em(:,:) * poa_op_ef(l)
cgwilson     &                          + poa_gv_em(:,:) * poa_gv_ef(l)
cgwilson     &                          + poa_dv_em(:,:) * poa_dv_ef(l)
cgwilson          aremis(:,:,kpcp_c(l)) = poa_mc_em(:,:) * poa_mc_ef(l)
cgwilson          aremis(:,:,kpfp_c(l)) = poa_bb_em(:,:) * poa_bb_ef(l)
cgwilson        enddo
cgwilson      endif
c
c  ---- deallocate the temporary array ---
c
      deallocate( argrid )
c
      goto 9999
c
c-----------------------------------------------------------------------
c    Error messages:
c-----------------------------------------------------------------------
c
 7000 continue
      write(iout,'(//,a)') 'ERROR in NCF_READAR:'
      write(iout,'(A)') action(:istrln(action))
      write(iout,'(A)') 'Cannot find the dimension id for layer dimension.'
      call camxerr()
c
 7001 continue
      write(iout,'(//,a)') 'ERROR in NCF_READAR:'
      write(iout,'(A)') action(:istrln(action))
      write(iout,'(2A)') 'Cannot read data for variable: ',
     &                                      this_var(:istrln(this_var))
      call camxerr()
c
 7002 continue
      write(iout,'(//,a)') 'ERROR in NCF_READAR:'
      write(iout,'(A)') action(:istrln(action))
      write(iout,'(A)') 'Number of layers does not match: '
      write(iout,'(A,I4)') 'User supplied: ',nlay
      write(iout,'(A,I4)') 'Value in file: ',nlays_in
      call camxerr()
c
 9999 continue
      return
      end
