!
#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_SOX
!
! !DESCRIPTION: data and methods for SOx emissions.
!
!
!    AR5 (August 2010): 
!        - SO2 emissions are used and split into SO2/SO4 
!        - in case of M7, additional partitioning is done into given modes

!    E.V.: 2.5 % of SO2 emissions are supposed to be SO4 (AEROCOM) to 
!          take into account the sulphate formation in the plumes, 
!          as sub-grid effect. For SO2 the remaining 97.5 % are 
!          calculated in the do_add routine. For SO4 they are calculated 
!          in the main emission_sox routine 
!
!          SO4 emissions:
!               industrial: (From Stier et al. 2004)
!                 100% accumulation mode number 
!                 median radius = 0.075 um sigma = 1.59
!               domestic - transport - biomass burning :
!                 50% Aitken mode number
!                 median radius = 0.03 um sigma = 1.59
!                 50% Accumulation mode number
!                 median radius = 0.075 um sigma = 1.59
!>>> TvN
!        The size distributions have been changed following
!        the AeroCom recommendations of Dentener et al. (2006),
!        adapted to the M7 modes (Stier et al., ACP, 2005).
!        For emissions from industry, power plants and shipping
!        a number median radius of 0.5 um was recommended,
!        at the boundary between accumulation and coarse mode.
!        Therefore, the emitted mass from these sources 
!        is equally divided between these modes,
!        with median radii of 0.075 and 0.75 um, respectively.
!        Emissions from road and off-road transport and the domestic sector
!        are assumed to be in the Aitken mode,
!        and from biomass burning in the accumulation mode.
!        Since only 2.5% of SOx emissions are emitted as particles,
!        these choices have only minor impacts.
!        See also comments in mo_aero.F90.
!<<< TvN
!\\
!\\
! !INTERFACE: 
!
MODULE EMISSION_SOX
  !
  ! !USES:
  !
  use GO,            only : gol, goErr, goPr
  use dims,          only : nregions, idate, okdebug
  use global_types,  only : emis_data, d3_data
  use emission_read, only : used_providers, has_emis

  use tm5_distgrid,  only : dgrid, get_distgrid, scatter
  use partools,      only : isRoot, par_broadcast

  implicit none
  private
  !
  ! !PUBLIC MEMBER FUNCTIONS:
  !
  public   :: Emission_SOx_Init    ! allocate dataset
  public   :: Emission_SOx_Done    ! deallocate dataset
  public   :: Emission_sox_declare ! read monthly input
  public   :: Emission_sox_apply   ! distribute & add emissions to tracer array
  !
  ! !PRIVATE DATA MEMBERS:
  !
  character(len=*), parameter  ::  mname = 'emission_sox'

  type( emis_data ), dimension(:,:), allocatable :: so2_emis_2d
  type( d3_data   ), dimension(:,:), allocatable :: so2_emis_3d

  integer              :: so2_3dsec,      so2_2dsec
  logical, allocatable :: has_data_3d(:), has_data_2d(:)

  !
  ! !REVISION HISTORY: 
  !          2005 - Elisabetta Vignati - changed for coupling with M7
  !    1 Oct 2010 - Achim Strunk  - AR5
  !    1 Dec 2011 - Narcisa Banda - added EDGAR 4
  !   27 Jun 2012 - P. Le Sager   - adapted for lon-lat MPI domain decomposition
  !
  ! !REMARKS:
  !
  !EOP
  !------------------------------------------------------------------------

CONTAINS

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    EMISSION_SOX_INIT
  !
  ! !DESCRIPTION: allocate memory
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE EMISSION_SOX_INIT( status )
    !
    ! !USES:
    !
    use dims,          only : lm
    use emission_read, only : providers_def, numb_providers, ed42_nsect_so2
    use emission_data, only : LAR5BMB
    use emission_read, only : n_ar5_ant_sec, n_ar5_shp_sec, n_ar5_air_sec, n_ar5_bmb_sec
    use emission_read, only : ar5_cat_ant,   ar5_cat_shp,   ar5_cat_air,   ar5_cat_bmb
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)        ::  status
    !
    ! !REVISION HISTORY: 
    !    1 Oct 2010 - Achim Strunk - AR5 format
    !   27 Jun 2012 - P. Le Sager  - adapted for lon-lat MPI domain decomposition
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

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

    status = 0
    if(.not. has_emis) return

    ! nb of sectors
    so2_2dsec = 0
    so2_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. 'AR5') then
             
             ! nb of available sectors in AR5 depends on category
             so2_2dsec = so2_2dsec + n_ar5_ant_sec*count('SO2'.eq.ar5_cat_ant) + &
                                     n_ar5_shp_sec*count('SO2'.eq.ar5_cat_shp)

             if (LAR5BMB) so2_2dsec = so2_2dsec + n_ar5_bmb_sec*count('SO2'.eq.ar5_cat_bmb)

             so2_3dsec = so2_3dsec + n_ar5_air_sec*count('SO2'.eq.ar5_cat_air)

!             so2_2dsec = so2_2dsec + providers_def(lprov)%nsect2d
!             so2_3dsec = so2_3dsec + count('SO2'.eq.ar5_cat_air)
             
          elseif (trim(providers_def(lprov)%name) .eq. 'ED42') then
             so2_2dsec = so2_2dsec + ed42_nsect_so2
             ! no 3d sectors in EDGAR 4.2
          else
             so2_2dsec = so2_2dsec + providers_def(lprov)%nsect2d
             so2_3dsec = so2_3dsec + providers_def(lprov)%nsect3d
          endif
       endif
    enddo
  
    allocate( so2_emis_2d( nregions, so2_2dsec ) )
    allocate( so2_emis_3d( nregions, so2_3dsec ) )
    allocate( has_data_2d(so2_2dsec)) ; has_data_2d=.false.
    allocate( has_data_3d(so2_3dsec)) ; 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 lsec=1,so2_2dsec
          allocate( so2_emis_2d(region,lsec)%surf(i1:i2,j1:j2) )
       end do

       do lsec=1,so2_3dsec
          allocate( so2_emis_3d(region,lsec)%d3(i1:i2,j1:j2,lmr) )
       end do

    enddo

    ! ok
    status = 0

  END SUBROUTINE EMISSION_SOX_INIT
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    EMISSION_SOX_DONE
  !
  ! !DESCRIPTION: Free memory.
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE EMISSION_SOX_DONE( status )
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)        ::  status
    !
    ! !REVISION HISTORY: 
    !    1 Oct 2010 - Achim Strunk - adapted for AR5
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

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

    status = 0
    if(.not. has_emis) return

    do region = 1, nregions
       do lsec=1,so2_2dsec
          deallocate( so2_emis_2d(region,lsec)%surf )
       end do
       do lsec=1,so2_3dsec
          deallocate( so2_emis_3d(region,lsec)%d3 )
       end do
    end do
    
    deallocate( so2_emis_2d, so2_emis_3d )
    deallocate( has_data_2d, has_data_3d )

    status = 0

  END SUBROUTINE EMISSION_SOX_DONE
  !EOC

  
  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    EMISSION_SOX_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. 
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE EMISSION_SOX_DECLARE( status )
    !
    ! !USES:
    !
    use binas,         only : pi
    use toolbox,       only : coarsen_emission
    use dims,          only : im, jm, lm, idate, sec_month, nlon360, nlat180, iglbsfc
    use chem_param,    only : xms, xmso2
    use emission_data, only : msg_emis, LAR5BMB
#ifdef with_m7
    use chem_param,    only : mode_ais, mode_acs, mode_cos
    use chem_param,    only : iso4ais, iso4acs, iso4cos
    use chem_param,    only : iais_n, iacs_n, icos_n
    use chem_param,    only : xmso4, xmnumb, sigma_lognormal
    use chem_param,    only : h2so4_factor
    use chem_param,    only : rad_so4_ait, rad_so4_acc, rad_so4_coa
    use chem_param,    only : frac_so4, so4_density
    use emission_data, only : emis_mass, emis_number
    use emission_data, only : emission_vdist_by_sector
#endif
    ! ---------------- AR5 - EDGAR 4 --------------------
    use emission_data, only : emis_input_year
    use emission_data, only : emis_input_dir_ed4,  emis_input_dir_mac
    use emission_data, only : emis_input_dir_gfed, emis_input_dir_retro
    use emission_read, only : emission_ar5_regrid_aircraft
    
    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_retro_ReadSector

    use emission_read, only : sectors_def, numb_sectors
    use emission_read, only : ar5_dim_3ddata
    use emission_read, only : ed42_so2_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
    !   27 Jun 2012 - P. Le Sager   - adapted for lon-lat MPI domain decomposition
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter :: rname = mname//'/emission_sox_declare'
    integer                     :: region, hasData
    integer, parameter          :: add_field=0
    integer, parameter          :: ayear=1, amonth=2
    integer                     :: imr, jmr, lmr, lsec, i1, i2, j1, j2
    type(d3_data)               :: work(nregions)
    type(emis_data)             :: wrk2D(nregions)

#ifdef with_m7                  
    real                        :: mass2numb_acs, mass2numb_ais, mass2numb_cos
    real                        :: numbscale, sfacacs, sfacais, sfaccos
    type(d3_data)               :: emis3d
#endif

    ! ---------------------------------------------------------------
    ! AR5
    real, dimension(:,:,:), allocatable  :: field3d,    field3d2
    integer                              :: seccount2d, seccount3d

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

    do region = 1, nregions
       do lsec=1,so2_2dsec
          so2_emis_2d(region,lsec)%surf = 0.0
       end do
       do lsec=1,so2_3dsec
          so2_emis_3d(region,lsec)%d3   = 0.0
       end do
    end do

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

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


    ! --------------------------------
    ! do a loop over available sectors
    ! --------------------------------

    ! count 2d and 3d sectors
    seccount2d = 0
    seccount3d = 0
    
    ! always allocate here 3d 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
    else
       allocate( field3d( 1, 1, 1 ) )
    end if


    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_so2_sectors.eq.sectors_def(lsec)%name) .eq. 0)) cycle

       if (associated(sectors_def(lsec)%species)) then  ! AR5 check
          if (count('SO2'.eq.sectors_def(lsec)%species).eq.0) cycle
          if ((trim(sectors_def(lsec)%catname) .eq. 'biomassburning').and.(.not.LAR5BMB)) cycle
       endif
       
       field3d = 0.0

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

       if (isRoot) then ! READ

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

          case( 'AR5' )
             
             ! Screen out agricultural and solvent sectors for SO2, 
             ! because they are zero in the RCPs
             ! and not present in the historical files.
             if (trim(sectors_def(lsec)%name) .ne. 'emiss_agr' .and. &
                 trim(sectors_def(lsec)%name) .ne. 'emiss_slv') then 
               call emission_ar5_ReadSector( 'SO2', emis_input_year, idate(2), lsec, field3d, status )
               IF_NOTOK_RETURN(status=1)
             endif

          case( 'MACC' )

             ! screen out sectors w/o SO2 (soil, bio, oc)
             if ( ( .not. (trim(sectors_def(lsec)%name) .eq. 'emiss_soil')) .and. &
                  ( .not. (trim(sectors_def(lsec)%name) .eq. 'emiss_bio') ) .and. &
                  ( .not. (trim(sectors_def(lsec)%name) .eq. 'emiss_air') ) .and. &
                  ( .not. (trim(sectors_def(lsec)%name) .eq. 'emiss_oc')  ) ) then 

                call emission_macc_ReadSector( emis_input_dir_mac, 'SO2', emis_input_year, idate(2), &
                                               '0.5x0.5_kg.nc', sectors_def(lsec)%name, 'kg / s', field3d, status )
                IF_NOTOK_RETURN(status=1)
               
             endif

          case( 'ED41' )

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

                call emission_ed4_ReadSector( emis_input_dir_ed4, 'SO2', 'so2', 2005, idate(2), &
                                              lsec, trim(sectors_def(lsec)%prov), 'kg / s', field3d, status )
                IF_NOTOK_RETURN(status=1)

                field3d = field3d * xmso2/xms
             end select

          case( 'ED42' )

             ! biomass burning (GFED/RETRO/AR5BMB) and transport (ED41) are excluded through ED42_SO2_SECTORS definition
             call emission_ed4_ReadSector( emis_input_dir_ed4, 'SO2', 'so2', emis_input_year, idate(2), &
                                           lsec, trim(sectors_def(lsec)%prov), 'kg / s', field3d, status )
             IF_NOTOK_RETURN(status=1)

             field3d = field3d * xmso2/xms

          case('GFEDv3')

             call emission_gfed_ReadSector( emis_input_dir_gfed, 'so2', emis_input_year, idate(2), &
                                            sectors_def(lsec)%name, 'kg / s', field3d(:,:,1), status )
             IF_NOTOK_RETURN(status=1)


          case('RETRO')

             call emission_retro_ReadSector( emis_input_dir_retro, 'SO2', emis_input_year, idate(2), &
                                             sectors_def(lsec)%name, 'kg / s', field3d(:,:,1), status )
             IF_NOTOK_RETURN(status=1)
	     
	  case('MEGAN')
             !
	     !   No biogenic emission of SO2 in the MEGAN inventory
	     !
          case('DUMMY')

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

          end select

          ! nothing found??? 
          if( sum(field3d) < 100.*TINY(1.0) ) then
             if (okdebug) then
                write(gol,'("EMISS-INFO - no SO2 emissions found for ",a," ",a," for month ",i2 )') &
                     trim(sectors_def(lsec)%prov), trim(sectors_def(lsec)%name), idate(2) ; call goPr
             endif
             hasData=0
          else
             if (okdebug) then
                write(gol,'("EMISS-INFO - found SO2 emissions for ",a," ",a," for month ",i2 )') &
                     trim(sectors_def(lsec)%prov), trim(sectors_def(lsec)%name), idate(2) ; call goPr
             endif
             field3d = field3d * sec_month ! from kg(SO2)/s to kg(SO2)/month
             hasData=1
          end if

       end if

       call Par_broadcast(hasData, status)
       IF_NOTOK_RETURN(status=1)

       if (hasData == 0) then
          cycle sec
       else
          if ( sectors_def(lsec)%f3d ) then
             has_data_3d(seccount3d)=.true.
          else
             has_data_2d(seccount2d)=.true.
          end if
       end if

       ! Distinguish 2d/3d sectors
       if( sectors_def(lsec)%f3d ) then 

          write(gol,'("EMISS-ERROR - Unexpected 3D data: Uncomment code below ")'); call goErr
          status=1; TRACEBACK; return

          ! ---------------------------
          ! 3d data (AIRCRAFT)
          ! ---------------------------

!           if (isRoot) then ! regrid
             
!              ! helper array for regridding
!              allocate( field3d2( nlon360,nlat180,lm(1) ) ) ; field3d2 = 0.0

!              ! aircraft data: regrid vertically to model layers
!              call emission_ar5_regrid_aircraft( iglbsfc, field3d, nlon360, nlat180, ar5_dim_3ddata, lm(1), field3d2, status )
!              IF_NOTOK_RETURN(status=1;deallocate(field3d,field3d2))
             
!              call msg_emis( amonth, trim(sectors_def(lsec)%prov)//' '//sectors_def(lsec)%name//' mass month', &
!                   & 'SO2', xmso2, sum(field3d2) )

!              ! distribute to so2_emis in regions
!              call Coarsen_Emission( 'SO2 '//trim(sectors_def(lsec)%name), nlon360,  nlat180, lm(1),            & 
!                                                                           field3d2, work,    add_field, status )
!              IF_NOTOK_RETURN(status=1;deallocate(field3d,field3d2))

!              deallocate( field3d2 ) 

!           end if

!           do region = 1, nregions
!              call scatter(dgrid(region), so2_emis_3d(region,seccount3d)%d3, work(region)%d3, 0, status)
!              IF_NOTOK_RETURN(status=1)
!           end do

       else ! ar5_sector is 2d

          ! ---------------------------
          ! 2d data (Anthropogenic, Ships, Biomassburning)
          ! ---------------------------
          if (isRoot) then ! print total & regrid

             call msg_emis( amonth, trim(sectors_def(lsec)%prov),sectors_def(lsec)%name, 'SO2', xmso2, sum(field3d(:,:,1)) )
             
             call coarsen_emission( 'SO2 '//sectors_def(lsec)%name,                            &
                                    nlon360, nlat180, field3d(:,:,1), wrk2D, add_field, status )
             IF_NOTOK_RETURN(status=1)

          end if

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

       end if ! 2D/3D

    end do sec ! sectors

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


#ifdef with_m7

    ! do SO4 partitioning from SO2 emissions

    ! --------------------------------
    ! do a loop over available sectors
    ! --------------------------------

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

    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_so2_sectors.eq.sectors_def(lsec)%name) .eq. 0)) cycle 
       
       if (associated(sectors_def(lsec)%species)) then  ! AR5 check
          if (count('SO2'.eq.sectors_def(lsec)%species).eq.0) cycle
          if ((trim(sectors_def(lsec)%catname) .eq. 'biomassburning').and.(.not.LAR5BMB)) cycle
       endif
       
       if( sectors_def(lsec)%f3d ) then 
          seccount3d = seccount3d + 1
          if (.not.has_data_3d(seccount3d)) cycle
       else 
          seccount2d = seccount2d + 1
          if (.not.has_data_2d(seccount2d)) cycle
       end if

!>>> TvN
       select case( trim(sectors_def(lsec)%vdisttype) )
       !case( 'industry' )
       case ( 'industry', 'combenergy' )
          sfacais = 0.0 * xmso4 / xmso2 * frac_so4
          !sfacacs = 1.0 * xmso4 / xmso2 * frac_so4
          sfacacs = 0.5 * xmso4 / xmso2 * frac_so4
          sfaccos = 0.5 * xmso4 / xmso2 * frac_so4
       case ( 'volcanic' )
         sfacais = 0.5 * xmso4 / xmso2 * frac_so4
         sfacacs = 0.5 * xmso4 / xmso2 * frac_so4
         sfaccos = 0.0 * xmso4 / xmso2 * frac_so4
       case default
         select case ( trim(sectors_def(lsec)%catname) )
         case ('ships')
           ! as industry and power plants
           sfacais = 0.0 * xmso4 / xmso2 * frac_so4
           sfacacs = 0.5 * xmso4 / xmso2 * frac_so4
           sfaccos = 0.5 * xmso4 / xmso2 * frac_so4
         case ( 'biomassburning')
           sfacais = 0.0 * xmso4 / xmso2 * frac_so4
           sfacacs = 1.0 * xmso4 / xmso2 * frac_so4
           sfaccos = 0.0 * xmso4 / xmso2 * frac_so4
         case default
           ! e.g. traffic, domestic
           sfacais = 1.0 * xmso4 / xmso2 * frac_so4
           sfacacs = 0.0 * xmso4 / xmso2 * frac_so4
           sfaccos = 0.0 * xmso4 / xmso2 * frac_so4
         end select
       end select
!<<< TvN

       ! mass to number factors for these two modes
       !>>> TvN
       ! The mass2numb conversion factors have been
       ! slightly enhanced to account for the small contribution
       ! from the hydrogen ions to the mass, 
       ! since so4_density is the density of sulfuric acid.
       numbscale     = rad_so4_ait*EXP(1.5*(LOG(sigma_lognormal(mode_ais)))**2) 
       mass2numb_ais = h2so4_factor*3./(4.*pi*numbscale**3*so4_density)
       numbscale     = rad_so4_acc*EXP(1.5*(LOG(sigma_lognormal(mode_acs)))**2) 
       mass2numb_acs = h2so4_factor*3./(4.*pi*numbscale**3*so4_density)
       numbscale     = rad_so4_coa*EXP(1.5*(LOG(sigma_lognormal(mode_cos)))**2)
       mass2numb_cos = h2so4_factor*3./(4.*pi*numbscale**3*so4_density)
       !<<< TvN


       do region = 1, nregions

          ! allocate helper array
          call get_distgrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )          
          lmr = lm(region)
          allocate( field3d(i1:i2,j1:j2,lmr) )

          ! distinguish 2d/3d sectors
          if( sectors_def(lsec)%f3d ) then 

             field3d = so2_emis_3d(region,seccount3d)%d3

          else 

             ! helper struct for vertical distribution
             allocate( emis3d%d3(i1:i2,j1:j2,lm(region)) ) ; emis3d%d3 = 0.0
             
             ! do vertical distribution
             call emission_vdist_by_sector( sectors_def(lsec)%vdisttype, 'SO2', region, so2_emis_2d(region,seccount2d), &
                                            emis3d, status )
             IF_NOTOK_RETURN(status=1;deallocate(field3d,emis3d%d3))

             field3d = emis3d%d3

             deallocate(emis3d%d3)

          end if

          ! finally add it to the emis target arrays according to the splitting above
          ! AIS
          emis_mass  (region,mode_ais)%d4(:,:,:,1) = &
               emis_mass  (region,mode_ais)%d4(:,:,:,1) + field3d * sfacais
          ! scale mass to get number
          emis_number(region,mode_ais)%d4(:,:,:,1) = &
               emis_number(region,mode_ais)%d4(:,:,:,1) + field3d * sfacais * mass2numb_ais
          ! ACS
          emis_mass  (region,mode_acs)%d4(:,:,:,1) = &
               emis_mass  (region,mode_acs)%d4(:,:,:,1) + field3d * sfacacs
          ! scale mass to get number
          emis_number(region,mode_acs)%d4(:,:,:,1) = &
               emis_number(region,mode_acs)%d4(:,:,:,1) + field3d * sfacacs * mass2numb_acs
!>>> TvN
          ! COS
          emis_mass  (region,mode_cos)%d4(:,:,:,1) = &
               emis_mass  (region,mode_cos)%d4(:,:,:,1) + field3d * sfaccos
          ! scale mass to get number
          emis_number(region,mode_cos)%d4(:,:,:,1) = &
               emis_number(region,mode_cos)%d4(:,:,:,1) + field3d * sfaccos * mass2numb_cos
!<<< TvN

          deallocate( field3d )

       end do ! regions

    end do ! sectors

#endif /* ifdef M7 */

  END SUBROUTINE EMISSION_SOX_DECLARE
  !EOC

  
  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    EMISSION_SOX_APPLY
  !
  ! !DESCRIPTION: Take monthly emissions, and
  !               - split them vertically
  !               - apply time splitting factors 
  !               - add them up (add_3d) 
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE EMISSION_SOX_APPLY( region, status )
    !
    ! !USES:
    !
    use dims,          only : itaur, nsrce, tref
    use dims,          only : im, jm, lm
    use datetime,      only : tau2date
    use emission_data, only : emission_vdist_by_sector,   LAR5BMB
    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 : ed42_so2_sectors
    use chem_param,    only : inh3, xmnh3, xmn, frac_so4
    use emission_read, only : sectors_def, numb_sectors
    use chem_param,    only : xmso2, xms, iso2, iso4
    use chem_param,    only : xmh2so4
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in)        :: region
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)       :: status
    !
    ! !REVISION HISTORY: 
    !    1 Oct 2010 - Achim Strunk - AR5 version
    !   27 Jun 2012 - Ph. Le Sager - adapted for lon-lat MPI domain decomposition 
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter ::  rname = mname//'/emission_sox_apply'
    integer, dimension(6)       :: idater
    real                        :: dtime, fraction, sfrac
    integer                     :: imr, jmr, lmr, lsec, 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_sox_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_so2_sectors.eq.sectors_def(lsec)%name) .eq. 0)) cycle 
       
       if (associated(sectors_def(lsec)%species)) then  ! AR5 check
          if (count('SO2'.eq.sectors_def(lsec)%species).eq.0) cycle
          if ((trim(sectors_def(lsec)%catname) .eq. 'biomassburning').and.(.not.LAR5BMB)) cycle
       endif


       ! default: no additional splitting
       fraction = 1.0

       ! ----------------------------------------------------------------------------------------
       ! distinguish here between sectors and whether they should have additional splitting, e.g.,
       ! if( sectors_def(lsec)%catname == 'biomassburning' ) fraction = fraction * bb_frac etc...
       ! ----------------------------------------------------------------------------------------

       ! distinguish between 2d/3d sectors
       if( sectors_def(lsec)%f3d ) then 

          seccount3d = seccount3d + 1
          if (.not.has_data_3d(seccount3d)) cycle
          emis3d%d3 = so2_emis_3d(region,seccount3d)%d3

       else

          seccount2d = seccount2d + 1
          if (.not.has_data_2d(seccount2d)) cycle
          emis3d%d3 = 0.0 
          ! vertically distribute according to sector
          call emission_vdist_by_sector( sectors_def(lsec)%vdisttype, 'SO2', region, so2_emis_2d(region,seccount2d), emis3d, status)
          IF_NOTOK_RETURN(status=1;deallocate(emis3d%d3))
       end if


       ! SO2: sfrac = fraction * (1.-frac_so4)
       sfrac = fraction * (1.0 - frac_so4)
       ! 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, iso2, i1, j1, emis3d%d3, bb_cycle(region)%scalef, xmso2, xmso2, status, sfrac )
          IF_NOTOK_RETURN(status=1)
       else
          call do_add_3d( region, iso2, i1, j1, emis3d%d3, xmso2, xmso2, status, sfrac )
          IF_NOTOK_RETURN(status=1)
       end if
       

#ifndef with_m7
       ! SO4: sfrac = fraction * frac_so4
       sfrac = fraction * frac_so4
       ! 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, iso4, i1, j1, emis3d%d3, bb_cycle(region)%scalef, xmh2so4, xmso2, status, sfrac )
          IF_NOTOK_RETURN(status=1)
       else
          call do_add_3d( region, iso4, i1, j1, emis3d%d3, xmh2so4, xmso2, status, sfrac )
          IF_NOTOK_RETURN(status=1)
       end if
#endif 


    end do ! sectors_def

    deallocate( emis3d%d3 )

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

    ! OK 
    status = 0

  END  SUBROUTINE EMISSION_SOX_APPLY
  !EOC

  
END MODULE EMISSION_SOX