#define TRACEBACK write (gol,'("in ",a," (",a,", line",i5,")")') rname, __FILE__, __LINE__; call goErr
#define IF_NOTOK_RETURN(action) if (status/=0) then; TRACEBACK; action; return; end if
#define IF_ERROR_RETURN(action) if (status> 0) then; TRACEBACK; action; return; end if
!
#include "tm5.inc"
!
!-----------------------------------------------------------------------------
!                    TM5                                                     !
!-----------------------------------------------------------------------------
!BOP
!
! !MODULE:      EMISSION_NMVOC
!
! !DESCRIPTION: hold data and methods for NMVOC emissions.
!
!   Non-methane VOC's :
!   nmhc(ncb5) = (/ipar,ieth,iole,iald2,imgly,ich2o/)
!\\
!\\
! !INTERFACE: 
!
MODULE EMISSION_NMVOC
  !
  ! !USES:
  !
  use GO,            only : gol, goPr, goErr
  use tm5_distgrid,  only : dgrid, get_distgrid, scatter
  use partools,      only : isRoot, par_broadcast
  use dims,          only : nregions, okdebug
  use global_types,  only : emis_data, d3_data
  use chem_param,    only : ncb5
  use emission_read, only : used_providers, has_emis

  IMPLICIT NONE
  PRIVATE
  !
  ! !PUBLIC MEMBER FUNCTIONS:
  !
  public   ::  Emission_NMVOC_Init     ! allocate memory
  public   ::  Emission_NMVOC_Done     ! deallocate memory
  public   ::  Emission_NMVOC_Declare  ! read input data
  public   ::  Emission_NMVOC_Apply    ! add emissions to tracer array
  !
  ! !PRIVATE DATA MEMBERS:
  !
  character(len=*), parameter  ::  mname = 'emission_nmvoc'

  type( emis_data ),dimension(:,:,:),allocatable :: hc_emis_2d
  type( d3_data   ),dimension(:,:,:),allocatable :: hc_emis_3d

  logical, allocatable :: has_data_3d(:,:), has_data_2d(:,:)
  integer :: hc_2dsec, hc_3dsec
  !
  ! !REVISION HISTORY: 
  !    1 Oct 2010 - Achim Strunk  - overhaul for AR5
  !    1 Dec 2011 - Narcisa Banda - added EDGAR 4
  !   29 Jun 2012 - Ph. Le Sager  - adapted for lon-lat MPI domain decomposition
  !
  ! !REMARKS:
  !
  !EOP
  !------------------------------------------------------------------------

CONTAINS

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    EMISSION_NMVOC_INIT
  !
  ! !DESCRIPTION: Allocate memory
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE EMISSION_NMVOC_INIT( status )
    !
    ! !USES:
    !
    use dims,           only : lm
    use emission_read,  only : providers_def, numb_providers 
    use emission_read , only : ed42_nsect_hc
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)     ::  status
    !
    ! !REVISION HISTORY: 
    !    1 Oct 2010 - Achim Strunk - adapted for AR5
    !   29 Jun 2012 - Ph. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter ::  rname = mname//'/Emission_NMVOC_Init'
    integer                     ::  region, lsec, icb5
    integer                     ::  lmr, lprov, i1, i2, j1, j2

    ! --- begin --------------------------------------
    
    status = 0
    if(.not. has_emis) return

    ! nb of sectors
    hc_2dsec = 0
    hc_3dsec = 0
    do lprov = 1, numb_providers
         if (count(used_providers.eq.providers_def(lprov)%name)/=0) then 
              if (trim(providers_def(lprov)%name) .eq. 'ED42') then
                 hc_2dsec = hc_2dsec + ed42_nsect_hc
                 ! no 3d sectors in EDGAR 4.2
              else
                 hc_2dsec = hc_2dsec + providers_def(lprov)%nsect2d
                 hc_3dsec = hc_3dsec + providers_def(lprov)%nsect3d
              endif
         endif
    enddo

    allocate( hc_emis_2d( nregions, hc_2dsec, ncb5 ) )
    allocate( hc_emis_3d( nregions, hc_3dsec, ncb5 ) )

    allocate( has_data_2d(hc_2dsec, ncb5 ) ) ; has_data_2d=.false.
    allocate( has_data_3d(hc_3dsec, ncb5 ) ) ; has_data_3d=.false.


    ! allocate information arrays (2d and 3d)
    do region=1,nregions

       CALL GET_DISTGRID( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )          
       lmr = lm(region)

       do icb5 = 1, ncb5

          do lsec=1,hc_2dsec
             allocate( hc_emis_2d(region,lsec,icb5)%surf(i1:i2, j1:j2) )
          end do

          do lsec=1,hc_3dsec
             allocate( hc_emis_3d(region,lsec,icb5)%d3(i1:i2, j1:j2, lmr) )
          end do
          
       end do
    enddo

    status = 0

  END SUBROUTINE EMISSION_NMVOC_INIT
  !EOC


  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    EMISSION_NMVOC_DONE
  !
  ! !DESCRIPTION: Free memory
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE EMISSION_NMVOC_DONE( status )
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)          ::  status
    !
    ! !REVISION HISTORY: 
    !    1 Oct 2010 - Achim Strunk - adapted to new structures
    !   29 Jun 2012 - Ph. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter ::  rname = mname//'/Emission_NMVOC_Done'
    integer                     ::  region, lsec, icb5

    ! --- begin --------------------------------------
    
    status = 0
    if(.not. has_emis) return

    do region = 1, nregions
       do icb5 = 1, ncb5
          do lsec=1,hc_2dsec
             deallocate( hc_emis_2d(region,lsec,icb5)%surf )
          end do
          do lsec=1,hc_3dsec
             deallocate( hc_emis_3d(region,lsec,icb5)%d3 )
          end do
       end do
    end do

    deallocate( hc_emis_2d )
    deallocate( hc_emis_3d )
    deallocate( has_data_2d, has_data_3d)

    status = 0

  END SUBROUTINE EMISSION_NMVOC_DONE
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    EMISSION_NMVOC_DECLARE
  !
  ! !DESCRIPTION: Opens, reads and evaluates input files (per month). 
  !               Provides emissions on 2d/3d-arrays which are then added
  !               to tracers in routine *apply. 
  !               Fields are communicated to all procs. 
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE EMISSION_NMVOC_DECLARE( status )
    !
    ! !USES:
    !
    use toolbox,       only : coarsen_emission
    use dims,          only : im, jm, lm, idate, sec_month, nlon360, nlat180, iglbsfc
    use chem_param,    only : xmc, ncb5, nmhc, names
    use emission_data, only : msg_emis, LAR5BMB, LMEGAN
    ! ---------------- AR5 - EDGAR 4 - ETC. --------------------
    use emission_data, only : emis_input_year_nmvoc, emis_input_year_nat
    use emission_data, only : emis_input_dir_mac
    use emission_data, only : emis_input_dir_megan
    use emission_data, only : emis_input_dir_retro
    use emission_data, only : emis_input_dir_gfed
    use emission_data, only : emis_input_dir_ed4

    use emission_read, only : emission_ar5_regrid_aircraft

    use emission_read, only : emission_cmip6_ReadSector
    use emission_read, only : emission_cmip6bmb_ReadSector
    use emission_read, only : emission_ar5_ReadSector
    use emission_read, only : emission_macc_ReadSector
    use emission_read, only : emission_ed4_ReadSector
    use emission_read, only : emission_gfed_ReadSector
    use emission_read, only : emission_megan_ReadSector
    use emission_read, only : emission_retro_ReadSector    

    use emission_read, only : sectors_def, numb_sectors
    use emission_read, only : ar5_dim_3ddata
    use emission_read, only : emis_cmip6_voc_name
    use emission_read, only : emis_cmip6_aircraft_tot2voc
    use emission_read, only : emis_cmip6_aircraft_tl_tot2voc
    use emission_read, only : emis_cmip6bmb_nvoc, emis_cmip6bmb_voc_name
    use emission_read, only : emis_cmip6bmb_voc2cbm5
    use emission_read, only : emis_ar5_nvoc, emis_ar5_voc_name
    use emission_read, only : emis_ar5_voc2cbm5_default
    use emission_read, only : emis_ar5_voc2cbm5_biomassb
    use emission_read, only : emis_ar5_voc2cbm5_biogenic
    use emission_read, only : emis_macc_nvoc, emis_macc_voc_name
    use emission_read, only : emis_megan_voc2cbm5_biogenic
    use emission_read, only : emis_macc_voc2cbm5_default
    use emission_read, only : emis_macc_voc2cbm5_biomassb
    use emission_read, only : emis_macc_voc2cbm5_biogenic
    use emission_read, only : emis_megan_voc2cbm5_biogenic
    use emission_read, only : emis_megan_nvoc, emis_megan_voc_name
    use emission_read, only : emis_gfed_nvoc, emis_gfed_voc_name, emis_voc2cbm5_gfed
    use emission_read, only : emis_retro_voc_name
    use emission_read, only : ed42_hc_sectors
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)        ::  status
    !
    ! !REVISION HISTORY: 
    !    1 Oct 2010 - Achim Strunk  - revamped for AR5
    !    1 Dec 2011 - Narcisa Banda - added EDGAR 4
    !   29 Jun 2012 - Ph. Le Sager  - adapted for lon-lat MPI domain decomposition
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter ::  rname = mname//'/emission_nmvoc_declare'

    ! --- local ---------------------------------------

    integer                    :: region
    logical                    :: hasData(ncb5)
    integer, parameter         :: add_field = 0
    integer, parameter         :: amonth = 2
    integer                    :: imr, jmr, lmr, i1, i2, j1, j2
    integer                    :: lsec, ivoc, icb5, ilev

    ! AR5
    real,dimension(:,:,:),  allocatable  :: field3d
    real,dimension(:,:,:,:),allocatable  :: field4d
    type(d3_data), dimension(nregions)   :: emis3d, work, work3d
    type(emis_data)                      :: wrk2D(nregions)
    integer                              :: seccount2d, seccount3d
    real,dimension(:,:),    allocatable  :: voc2cbm5

    ! --- begin ----------------------------------------
    
    status = 0
    if(.not. has_emis) return
    
    write(gol,'(" EMISS-INFO ------------- read NMVOC emissions -------------")'); call goPr

    ! reset arrays
    do region = 1, nregions
       do icb5 = 1, ncb5
          do lsec=1,hc_2dsec
             hc_emis_2d(region,lsec,icb5)%surf = 0.0
          end do
          do lsec=1,hc_3dsec
             hc_emis_3d(region,lsec,icb5)%d3   = 0.0
          end do
       end do

       CALL GET_DISTGRID( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )
       lmr = lm(region)
 
       allocate( work3d(region)%d3  (i1:i2,j1:j2, ar5_dim_3ddata) ) ; work3d(region)%d3   = 0.0
       allocate( emis3d(region)%d3  (i1:i2,j1:j2, lmr           ) ) ; emis3d(region)%d3   = 0.0
    end do

    ! global arrays for coarsening    
    do region = 1, nregions
       if (isRoot)then
          allocate(work(region)%d3(im(region),jm(region),ar5_dim_3ddata))
       else
          allocate(work(region)%d3(1,1,1))
       end if
    enddo

    do region = 1, nregions
       wrk2D(region)%surf => work(region)%d3(:,:,1)
    end do

    ! count 2d and 3d sectors
    seccount2d = 0
    seccount3d = 0

    ! always allocate here 3D/4D data set (for 2d sectors it will be filled in first layer only)
    if (isRoot) then
       allocate( field3d( nlon360, nlat180, ar5_dim_3ddata       ) ) ; field3d = 0.0
       allocate( field4d( nlon360, nlat180, ar5_dim_3ddata, ncb5 ) ) ; field4d = 0.0
    else
       allocate( field3d( 1, 1, 1    ) )
       allocate( field4d( 1, 1, 1, 1 ) )
    end if

    ! --------------------------------
    ! do a loop over available sectors
    ! --------------------------------
    sec : do lsec = 1, numb_sectors

       if (count(used_providers.eq.sectors_def(lsec)%prov).eq.0) cycle 

       if ((trim(sectors_def(lsec)%prov).eq.'ED42') .and. (count(ed42_hc_sectors.eq.sectors_def(lsec)%name) .eq. 0)) cycle 

       if( sectors_def(lsec)%f3d ) then 
          seccount3d = seccount3d + 1
       else 
          seccount2d = seccount2d + 1
       end if

       field3d = 0.0
       field4d = 0.0

       if (isRoot) then ! READ

          select case( trim(sectors_def(lsec)%prov) )

          case( 'CMIP6' )

             ! get appropriate array of splitting values from input to cbm5
             allocate( voc2cbm5(emis_ar5_nvoc,ncb5) )
             voc2cbm5 = reshape( emis_ar5_voc2cbm5_default,  (/emis_ar5_nvoc,ncb5/) )

             if (trim(sectors_def(lsec)%catname).ne.'aircraft') then 

               ! --------------------------------
               ! do a loop over available constituents (25)
               ! --------------------------------
               do ivoc = 1, emis_ar5_nvoc

                  ! skip anthropogenic emissions for isoprene and terpenes,
                  ! which are not provided:
                  if (trim(emis_cmip6_voc_name(ivoc)).eq.'VOC10-isoprene') cycle 
                  if (trim(emis_cmip6_voc_name(ivoc)).eq.'VOC11-terpenes') cycle

                  call emission_cmip6_ReadSector( trim(emis_cmip6_voc_name(ivoc)), &
                                emis_input_year_nmvoc, idate(2), lsec, field3d, status )
                  IF_NOTOK_RETURN(status=1)
                  do icb5 = 1, ncb5
                     field4d(:,:,:,icb5) = field4d(:,:,:,icb5) + field3d(:,:,:) * voc2cbm5(ivoc,icb5)
                  end do
               end do

             else

                ! For the aviation sector only the total NMVOC emissions are provided
                call emission_cmip6_ReadSector('NMVOC', emis_input_year_nmvoc, idate(2), lsec, field3d, status ) 

                ! Apply takeoff and landing VOC profile 
                ! to the lowest two layers (0.305 km and 0.915 km) 
                ! in the aircraft emissions file:
                do ivoc = 1, emis_ar5_nvoc
                   do icb5 = 1, ncb5
                      do ilev = 1, 2
                         field4d(:,:,ilev,icb5) = field4d(:,:,ilev,icb5) + field3d(:,:,ilev) * &
                                                  emis_cmip6_aircraft_tl_tot2voc(ivoc) * voc2cbm5(ivoc,icb5)
                      end do
                      do ilev = 3, ar5_dim_3ddata
                         field4d(:,:,ilev,icb5) = field4d(:,:,ilev,icb5) + field3d(:,:,ilev) * &
                                                  emis_cmip6_aircraft_tot2voc(ivoc) * voc2cbm5(ivoc,icb5) 
                      end do
                   end do
                end do

             endif

             deallocate( voc2cbm5 )

          case( 'CMIP6BMB' )

             ! get appropriate array of splitting values from input to cbm5
             allocate( voc2cbm5(emis_cmip6bmb_nvoc,ncb5) )
             voc2cbm5 = reshape( emis_cmip6bmb_voc2cbm5,  (/emis_cmip6bmb_nvoc,ncb5/) )

             ! --------------------------------
             ! do a loop over available constituents
             ! --------------------------------
             do ivoc = 1, emis_cmip6bmb_nvoc

                ! skip species for which biomass burning emissions are not provided
                if (trim(emis_cmip6bmb_voc_name(ivoc)).eq.'not provided') cycle
                ! skip isoprene or terpenes,
                ! as their contributions will be zero anyhow
                if (trim(emis_cmip6bmb_voc_name(ivoc)).eq.'C5H8') cycle
                if (trim(emis_cmip6bmb_voc_name(ivoc)).eq.'C10H16') cycle

                call emission_cmip6bmb_ReadSector( 'NMVOC-'//trim(emis_cmip6bmb_voc_name(ivoc)), &
                                      emis_input_year_nmvoc, idate(2), lsec, field3d, status )
                IF_NOTOK_RETURN(status=1)

                do icb5 = 1, ncb5
                   field4d(:,:,:,icb5) = field4d(:,:,:,icb5) + field3d(:,:,:) * voc2cbm5(ivoc,icb5)
                end do

             end do

             deallocate( voc2cbm5 )

          case( 'AR5' )

             ! screen out AR5 biomass burning if not wanted
             if ( .not. ((trim(sectors_def(lsec)%catname) .eq. 'biomassburning').and.(.not.LAR5BMB)) ) then

                ! get appropriate array of splitting values from input to cbm5
                allocate( voc2cbm5(emis_ar5_nvoc,ncb5) )

                ! Changing the tables here in order to get ALD and PAR etc from Biomassburning
                voc2cbm5 = reshape( emis_ar5_voc2cbm5_default,  (/emis_ar5_nvoc,ncb5/) )

                ! --------------------------------
                ! do a loop over available constituents (25)
                ! --------------------------------
                vocAR5: do ivoc = 1, emis_ar5_nvoc 
                   
                   ! skip missing categories
                   if (count(trim(emis_ar5_voc_name(ivoc)).eq.sectors_def(lsec)%species).eq.0) cycle

                   
                   call emission_ar5_ReadSector( trim(emis_ar5_voc_name(ivoc)), emis_input_year_nmvoc, idate(2), lsec, field3d, status )
                   IF_NOTOK_RETURN(status=1)

                   do icb5 = 1, ncb5
                      field4d(:,:,:,icb5) = field4d(:,:,:,icb5) + field3d(:,:,:) * voc2cbm5(ivoc,icb5)
                   end do

                end do vocAR5

                deallocate( voc2cbm5 )
             end if

          case( 'MACC' )

             ! screen out biomass burning (a/k/a emiss_bb), and 'soil', 'nat' and 'air' sectors (no NMVOC)
             ! skip 'bio' source if already provided by MEGAN
             if ( ( .not. (trim(sectors_def(lsec)%catname) .eq. 'biomassburning') ) .and. &
                  ( .not. (trim(sectors_def(lsec)%name)    .eq. 'emiss_soil')     ) .and. &
                  ( .not. (trim(sectors_def(lsec)%name)    .eq. 'emiss_nat')      ) .and. &
                  ( .not. (trim(sectors_def(lsec)%name)    .eq. 'emiss_air')      ) .and. &
                  ( .not. (LMEGAN .and. (trim(sectors_def(lsec)%name) .eq. 'emiss_bio'))) ) then
                
                ! get appropriate array of splitting values from input to cbm5
                allocate( voc2cbm5(emis_macc_nvoc,ncb5) )
                select case( trim(sectors_def(lsec)%catname) )
                case( 'anthropogenic', 'ships', 'aircraft' )
                   voc2cbm5 = reshape( emis_macc_voc2cbm5_default,  (/emis_macc_nvoc,ncb5/) )
                case( 'natural' )
                   voc2cbm5 = reshape( emis_macc_voc2cbm5_biogenic, (/emis_macc_nvoc,ncb5/) )
                case default
                   write(gol,'(80("-"))') ; call goPr
                   write(gol,'("ERROR: wrong category name `",a,"`!")') sectors_def(lsec)%catname ; call goErr
                   write(gol,'(80("-"))') ; call goPr
                end select

                ! --------------------------------
                ! do a loop over available constituents (14)
                ! --------------------------------
                do ivoc = 1, emis_macc_nvoc 

                   if (trim(sectors_def(lsec)%catname) .eq. 'natural') then
                      call emission_macc_ReadSector( emis_input_dir_mac, trim(emis_macc_voc_name(ivoc)), emis_input_year_nat, idate(2), &
                           '0.5x0.5_kg.nc', sectors_def(lsec)%name, 'kg / s', field3d, status )
                      IF_NOTOK_RETURN(status=1;deallocate(field3d,field4d))
                   else 
                      call emission_macc_ReadSector( emis_input_dir_mac, trim(emis_macc_voc_name(ivoc)), emis_input_year_nmvoc, idate(2), &
                           '0.5x0.5_kg.nc', sectors_def(lsec)%name, 'kg / s', field3d, status )
                      IF_NOTOK_RETURN(status=1;deallocate(field3d,field4d))
                   endif

                   do icb5 = 1, ncb5
                      field4d(:,:,:,icb5) = field4d(:,:,:,icb5) + field3d(:,:,:) * voc2cbm5(ivoc,icb5)
                   end do
                   
                end do
                
                deallocate( voc2cbm5 )
                
             endif

          case( 'ED41' )
             select case(trim(sectors_def(lsec)%name))
             case ('1A3b_c_e','1A3d_SHIP','1A3d1')

                ! AR5 NMVOC emissions and were used to split EDGAR NMVOC totals into separate (AR5) NMVOC species

                ! get appropriate array of splitting values from input to cbm5
                allocate( voc2cbm5(emis_ar5_nvoc,ncb5) )
                ! Changing the tables here in order to get ALD and PAR etc from Biomassburning
                voc2cbm5 = reshape( emis_ar5_voc2cbm5_default,  (/emis_ar5_nvoc,ncb5/) )

                ! --------------------------------
                ! do a loop over available constituents (25)
                ! --------------------------------
                vocED41: do ivoc = 1, emis_ar5_nvoc 
                   
                   ! screen out missing dataset
                   if ((ivoc==10).or.(ivoc==11)) cycle vocED41
                   if ( trim(sectors_def(lsec)%name) /= '1A3b_c_e') then
                      if ((ivoc==1).or.(ivoc>=18)) cycle vocED41
                   endif
                   if (((ivoc==1).or.(ivoc==24)).and.(emis_input_year_nmvoc<2005)) cycle vocED41 ! kludge (waiting for NB input)

                   
                   call emission_ed4_ReadSector( emis_input_dir_ed4, trim(emis_ar5_voc_name(ivoc)), trim(emis_ar5_voc_name(ivoc)),&
                                                 emis_input_year_nmvoc, idate(2), lsec,  trim(sectors_def(lsec)%prov), 'kg / s', field3d,&
                                                 status )
                   IF_NOTOK_RETURN(status=1;deallocate(field3d,field4d,voc2cbm5))

                   do icb5 = 1, ncb5
                      field4d(:,:,:,icb5) = field4d(:,:,:,icb5) + field3d(:,:,:) * voc2cbm5(ivoc,icb5)
                   end do

                end do vocED41
                deallocate( voc2cbm5 )

             end select

          case( 'ED42' )
             
             ! Biomass burning (GFED/RETRO/AR5BMB) and transport (ED41) are excluded through ED42_HC_SECTORS definition

             ! Same constituents as AR5 NMVOC are used
                   
             ! get appropriate array of splitting values from input to cbm5
             allocate( voc2cbm5(emis_ar5_nvoc,ncb5) )
             voc2cbm5 = reshape( emis_ar5_voc2cbm5_default,  (/emis_ar5_nvoc,ncb5/) )

             ! --------------------------------
             ! do a loop over available constituents (25)
             ! --------------------------------
             vocED42: do ivoc = 1, emis_ar5_nvoc 
                
                ! screen out missing dataset
                if ((ivoc==10).or.(ivoc==11)) cycle vocED42
                if ( trim(sectors_def(lsec)%name) == '3') then
                   if (count(ivoc.eq.(/2,3,4,5,7,8,9,12,13,16,21,22,24/)).eq.1) cycle vocED42
                endif

                if ( trim(sectors_def(lsec)%name) == '1A4') then
                   if (count(ivoc.eq.(/1,24/)).eq.1) cycle vocED42
                endif

                if ( trim(sectors_def(lsec)%name) == '4F') then
                   if (count(ivoc.eq.(/16,17,18,19,20,25/)).eq.1) cycle vocED42
                endif

                if ( ( trim(sectors_def(lsec)%name) == '1A1a')              .or. &
                     ( trim(sectors_def(lsec)%name) == '1A1b_c_1B_2C1_2C2') .or. &
                     ( trim(sectors_def(lsec)%name) == '1A2')               .or. &
                     ( trim(sectors_def(lsec)%name) == '2A_B_D_E_F_G')      .or. &
                     ( trim(sectors_def(lsec)%name) == '7A')                ) then
                   if (count(ivoc.eq.(/18,19,20/)).eq.1) cycle vocED42
                endif

                call emission_ed4_ReadSector( emis_input_dir_ed4, trim(emis_ar5_voc_name(ivoc)), &
                                              trim(emis_ar5_voc_name(ivoc)), emis_input_year_nmvoc, idate(2), lsec,  &
                                              trim(sectors_def(lsec)%prov), 'kg / s', field3d, status )
                IF_NOTOK_RETURN(status=1)
                do icb5 = 1, ncb5
                   field4d(:,:,:,icb5) = field4d(:,:,:,icb5) + field3d(:,:,:) * voc2cbm5(ivoc,icb5)
                end do

             end do vocED42
             deallocate( voc2cbm5 )

          case( 'GFEDv3' )

             ! We use the MACC voc split, since GFED3 contains the same species, 
             ! except for the last 3 species - Acetone, Acetaldehyde and MEK - which are not provided in GFED3

             ! get appropriate array of splitting values from input to cbm5
             allocate( voc2cbm5(emis_gfed_nvoc,ncb5) )
             voc2cbm5 = reshape( emis_voc2cbm5_gfed, (/emis_gfed_nvoc,ncb5/) )

             ! --------------------------------
             ! do a loop over available constituents (11)
             ! --------------------------------
             do ivoc = 1, emis_gfed_nvoc

                call emission_gfed_ReadSector(  emis_input_dir_gfed, trim(emis_gfed_voc_name(ivoc)), emis_input_year_nmvoc, idate(2), &
                                                sectors_def(lsec)%name, 'kg / s', field3d(:,:,1), status )
                IF_NOTOK_RETURN(status=1;deallocate(field3d,field4d))

                do icb5 = 1, ncb5
                   field4d(:,:,:,icb5) = field4d(:,:,:,icb5) + field3d(:,:,:) * voc2cbm5(ivoc,icb5)
                end do

             end do

             deallocate( voc2cbm5 )


          case( 'RETRO' )

             allocate( voc2cbm5(emis_ar5_nvoc,ncb5) )
             voc2cbm5 = reshape( emis_ar5_voc2cbm5_biomassb,  (/emis_ar5_nvoc,ncb5/) )

             ! --------------------------------
             ! do a loop over available constituents (25)
             ! --------------------------------
             do ivoc = 1, emis_ar5_nvoc

                if (trim(emis_retro_voc_name(ivoc)).eq.'not provided') cycle

                call emission_retro_ReadSector( emis_input_dir_retro, trim(emis_retro_voc_name(ivoc)), emis_input_year_nmvoc, idate(2), &
                     sectors_def(lsec)%name, 'kg / s', field3d(:,:,1), status )
                IF_NOTOK_RETURN(status=1;deallocate(field3d))

                do icb5 = 1, ncb5
                   field4d(:,:,:,icb5) = field4d(:,:,:,icb5) + field3d(:,:,:) * voc2cbm5(ivoc,icb5)
                end do

             end do

             deallocate( voc2cbm5 )
	     
	  case( 'MEGAN' ) 
	  
             allocate( voc2cbm5(emis_megan_nvoc,ncb5) )
             voc2cbm5 = reshape( emis_megan_voc2cbm5_biogenic,  (/emis_megan_nvoc,ncb5/) )

             ! --------------------------------
             ! do a loop over available constituents (25)
             ! --------------------------------
             do ivoc = 1, emis_megan_nvoc

                if (trim(emis_megan_voc_name(ivoc)).eq.'not provided') cycle

                call emission_megan_ReadSector( emis_input_dir_megan, trim(emis_megan_voc_name(ivoc)), emis_input_year_nat, idate(2), &
                     sectors_def(lsec)%name, 'kg / s', field3d(:,:,1), status )
                IF_NOTOK_RETURN(status=1;deallocate(field3d))

                do icb5 = 1, ncb5
                   field4d(:,:,:,icb5) = field4d(:,:,:,icb5) + field3d(:,:,:) * voc2cbm5(ivoc,icb5)
                end do

             end do

             deallocate( voc2cbm5 )

          case('DUMMY')

          case default
             write(gol,*) "Error in buidling list of providers USED_PROVIDERS"; call goErr
             status=1; TRACEBACK; return

          END SELECT

          ! nothing found?
          do icb5 = 1, ncb5

             if( sum(field4d(:,:,:,icb5)) < 100.*TINY(1.0) ) then
                if (okdebug) then
                   write(gol,'("EMISS-INFO - no NMVOC emissions found for ",a," ",a," for month ",i2 )') &
                        trim(sectors_def(lsec)%prov), trim(sectors_def(lsec)%name), idate(2) ; call goPr
                endif
                hasData(icb5)=.false.
             else
                if (okdebug) then
                   write(gol,'("EMISS-INFO - found NMVOC emissions for ",a," ",a," for month ",i2 )') &
                        trim(sectors_def(lsec)%prov), trim(sectors_def(lsec)%name), idate(2) ; call goPr
                endif
                field4d(:,:,:,icb5) = field4d(:,:,:,icb5) * sec_month ! from kg/s to kg/month
                hasData(icb5)=.true.
             end if
             
          enddo
          
       end if

       call Par_broadcast(hasData, status)
       IF_NOTOK_RETURN(status=1)
       
!       if (.not.(any(hasData))) cycle sec ! early exit?
       
       if ( sectors_def(lsec)%f3d ) then
          has_data_3d(seccount3d,:)=hasData
       else
          has_data_2d(seccount2d,:)=hasData
       end if

       
       ! Loop over cb5 components and distinguish b/w 2d/3d sectors       
       do icb5 = 1, ncb5
          if( sectors_def(lsec)%f3d ) then 

             if (has_data_3d(seccount3d,icb5)) then

                ! ---------------------------------------
                ! 3d data (AIRCRAFT), available for CMIP6
                ! ---------------------------------------

                if (isRoot) then

                   ! write some numbers
                   call msg_emis( amonth, trim(sectors_def(lsec)%prov), sectors_def(lsec)%name, &
                                  trim(names(nmhc(icb5))), xmc, sum(field4d(:,:,:,icb5)) )

                   ! distribute to work arrays in regions
                   call Coarsen_Emission( trim(names(nmhc(icb5)))//trim(sectors_def(lsec)%name), &
                                          nlon360,  nlat180, ar5_dim_3ddata, field4d(:,:,:,icb5),  work,    add_field, status ) 
                   IF_NOTOK_RETURN(status=1)

                end if

                ! scatter, sum up on target array
                do region = 1, nregions

                   call scatter(dgrid(region), work3d(region)%d3, work(region)%d3, 0, status)
                   IF_NOTOK_RETURN( status=1 )
               
                   CALL GET_DISTGRID( dgrid(region), I_STRT=i1, J_STRT=j1)

                   ! aircraft data: regrid vertically to model layers
                   call emission_ar5_regrid_aircraft( region, i1, j1, work3d(region)%d3, emis3d(region)%d3, status )
                   IF_NOTOK_RETURN( status=1 )

                   hc_emis_3d(region,seccount3d,icb5)%d3 = hc_emis_3d(region,seccount3d,icb5)%d3 + emis3d(region)%d3

                end do

             endif

          else
             ! ---------------------------
             ! 2d data (Anthropogenic, Ships, Biomassburning)
             ! ---------------------------
             if (has_data_2d(seccount2d,icb5)) then 
             
                if (isRoot) then ! print total & regrid
                
                   call msg_emis( amonth, trim(sectors_def(lsec)%prov),sectors_def(lsec)%name, &
                                  trim(names(nmhc(icb5))), xmc, sum(field4d(:,:,1,icb5)) )

                   call coarsen_emission( trim(names(nmhc(icb5)))//sectors_def(lsec)%name, &
                                          nlon360, nlat180, field4d(:,:,1,icb5), wrk2D, add_field, status )
                   IF_NOTOK_RETURN(status=1)

                end if

                do region = 1, nregions
                   call scatter(dgrid(region), hc_emis_2d(region,seccount2d,icb5)%surf, work(region)%d3(:,:,1), 0, status)
                   IF_NOTOK_RETURN(status=1)
                end do

             endif
          endif ! 2D/3D
       enddo
       
    end do sec ! sectors


    deallocate( field3d, field4d )
    do region = 1, nregions
       if (associated(wrk2D(region)%surf)) nullify(wrk2D(region)%surf)
       deallocate( work(region)%d3 )
       deallocate( work3d(region)%d3 )
       deallocate( emis3d(region)%d3 ) 
    end do
    
    ! check sectors found
    if( seccount2d /= hc_2dsec ) then 
       write(gol,'(80("-"))') ; call goPr
       write(gol,'("ERROR: 2d sectors do not equal total number:",i4," /= ",i4," !")') seccount2d, hc_2dsec ; call goErr
       write(gol,'(80("-"))') ; call goPr
       status=1; return
    end if
    if( seccount3d /= hc_3dsec ) then 
       write(gol,'(80("-"))') ; call goPr
       write(gol,'("ERROR: 3d sectors do not equal total number:",i4," /= ",i4," !")') seccount3d, hc_3dsec ; call goErr
       write(gol,'(80("-"))') ; call goPr
       status=1; return
    end if
    
    status = 0

  end subroutine emission_nmvoc_declare
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:   EMISSION_NMVOC_APPLY
  !
  ! !DESCRIPTION: Take monthly emissions, and
  !               - split them vertically
  !               - apply time splitting factors 
  !               - add them up (add_3d) 
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE EMISSION_NMVOC_APPLY( region, status )
    !
    ! !USES:
    !
    use dims,          only : idate, itaur, nsrce, tref
    use dims,          only : im, jm, lm
    use chem_param,    only : xmcb5, xmc, ntracet, nmhc, names
    use datetime,      only : tau2date
    use emission_data, only : emission_vdist_by_sector
    use emission_data, only : do_add_3d, do_add_3d_cycle, bb_cycle
    use emission_data, only : emis_bb_trop_cycle
    use emission_read, only : sectors_def, numb_sectors
    use emission_read, only : ed42_hc_sectors
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in)         :: region
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)        :: status
    !
    ! !REVISION HISTORY: 
    !    1 Oct 2010 - Achim Strunk - rewritten for AR5
    !   29 Jun 2012 - Ph. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC
    
    character(len=*), parameter ::  rname = mname//'/emission_nmvoc_apply'

    ! --- local ---------------------------------------

    integer,dimension(6)       :: idater
    real                       :: dtime, fraction
    integer                    :: imr, jmr, lmr, lsec, icb5, i1, i2, j1, j2
    integer                    :: seccount2d, seccount3d
    type(d3_data)              :: emis3d

    ! --- begin -----------------------------------------

    status = 0
    if(.not. has_emis) return

    if( okdebug ) then
       write(gol,*) 'start of emission_nmvoc_apply'; call goPr
    end if

    call tau2date(itaur(region),idater)
    dtime=float(nsrce)/(2*tref(region))    !emissions are added in two steps...XYZECCEZYX.


    ! get a working structure for 3d emissions
    call get_distgrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )
    allocate( emis3d%d3(i1:i2,j1:j2,lm(region)) ) ; emis3d%d3 = 0.0


    ! count 2d and 3d sectors
    seccount2d = 0
    seccount3d = 0

    ! cycle over sectors
    do lsec = 1, numb_sectors

       if (count(used_providers.eq.sectors_def(lsec)%prov).eq.0) cycle 

       if ((trim(sectors_def(lsec)%prov).eq.'ED42') .and. (count(ed42_hc_sectors.eq.sectors_def(lsec)%name) .eq. 0)) cycle

       if( sectors_def(lsec)%f3d ) then ! count
          seccount3d = seccount3d + 1
       else 
          seccount2d = seccount2d + 1
       end if
       
       fraction = 1.0           ! default: no additional splitting

       ! ----------------------------------------------------------------------------------------
       ! distinguish here between sectors and whether they should have additional splitting 
       ! if( sectors_def(lsec)%catname == 'biomassburning' ) fraction = fraction * bb_frac etc...
       ! ----------------------------------------------------------------------------------------
       cb05: do icb5=1,ncb5

          if( nmhc(icb5).gt. ntracet ) cycle cb05   ! only transported species, skip short lived species

          ! distinguish between 2d/3d sectors
          if( sectors_def(lsec)%f3d ) then 
             if (.not.has_data_3d(seccount3d,icb5)) cycle cb05
             emis3d%d3 = hc_emis_3d(region,seccount3d,icb5)%d3
          else
             if (.not.has_data_2d(seccount2d,icb5)) cycle cb05

             ! vertically distribute according to sector
             emis3d%d3 = 0.0
             call emission_vdist_by_sector( sectors_def(lsec)%vdisttype, trim(names(nmhc(icb5))), region, &
                                            hc_emis_2d(region,seccount2d,icb5), emis3d, status )
             IF_NOTOK_RETURN(status=1)
          endif

          ! add dataset according to sector and category
          if( emis_bb_trop_cycle .and. trim(sectors_def(lsec)%catname) == "biomassburning" ) then 
             call do_add_3d_cycle( region, nmhc(icb5), i1, j1, emis3d%d3, bb_cycle(region)%scalef, &
                                   xmcb5(icb5), xmcb5(icb5), status, fraction )
             IF_NOTOK_RETURN(status=1)
          else
             call do_add_3d( region, nmhc(icb5), i1, j1, emis3d%d3, xmcb5(icb5), xmcb5(icb5), status, fraction )
             IF_NOTOK_RETURN(status=1)
          endif
          
       enddo cb05
          
    enddo ! sectors

    deallocate( emis3d%d3 )

    if(okdebug) then
       write(gol,*) 'end of emission_nmvoc_apply'; call goPr
    end if

    status = 0

  END SUBROUTINE EMISSION_NMVOC_APPLY
  !EOC


END MODULE EMISSION_NMVOC