ecearth3/sources/tm5mp/proj/cb05/emission_pom.F90

!
#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_POM
!
! !DESCRIPTION: data and methods for Particulate Organic Matter (POM) emissions. 
!
!  POM is a sum of fossil and biofuel, vegetation fire, and SOA emissions.
! 
!        Dimension properties:
!        1) from fossil fuel = 0.015 um and sigma = 1.59
!        2) from vegetation fires = 0.04 um sigma = 1.59
!        3) POM from SOA has the characteristics of secondary 
!           products from monoterpine oxidation, which probably 
!           condense on pre-existing particles. Therefore only mass
!           is calculated and added to the POM mass field if there 
!           are particles in the one or more soluble modes, if not 
!           they are formed in the Aitken soluble mode. 
!           The distribution of the total condensing SOA mass (M) in 
!           modes is
!        
!        Solubility properties (from Stier et al.ACP, 2005):
!        1) POM from fossil fuel emissions is considered 65% soluble
!        2) POM from vegetation fire emissions is considered 65% soluble
!        3) POM from SOA is considered 100% soluble (Kanakidou et al. ACP, 2004)
!
!>>> TvN
! The latest revision allows to apply
! different sizes for the emissions from
! different sectors, as well as for
! biofuel emissions.
! Also, different solubility fractions can
! now be assumed for emissions from biomass burning, 
! biofuel combustion and fossil fuel combustion.
! See settings in emission_data.F90.
!
! The soluble and insoluble SOA fractions are assumed to 
! condense onto particles in the Aitken modes, 
! as in the previous version by E. Vignati.
! However, this assumption is questionable
! in regions where biomass burning is important, 
! since emissions from vegetation fires
! are now assumed to occur in the accumulation modes.
! (Stier et al. assume condensation 
! in the soluble Aitken and accumulation modes, 
! but it is not clear how the condensation is/should be distributed 
! over the two soluble modes in that case.)
!
! The formation of new particles has been removed.
! Previously, a cut-off was applied based on emis_number
! (i.e. the number of particles emitted in a gridbox per month)
! from all sectors excluding the surrogate SOA emissions.
! This criterion is unrealistic 
! in the sense that it does not describe any physical process, 
! but was only used to select the cells without primary OA emissions.
! We have tested that the inclusion of such a cutoff for new particle formation
! has only very marginal impacts on the results, and can therefore be removed
! (in any case when starting from reasonable initial concentrations).
!
! Note that in the ECHAM version described by Zhang et al. 
! fossil fuel and biofuel emissions are considered 
! 100% and 65% soluble, respectively.
! Since these sources are not distinguished in the currently used 
! emission data sets, we keep the solubility of both sources
! at 65%, as was done by Stier et al.
!<<< TvN
!
!\\
!\\
! !INTERFACE: 
!
MODULE EMISSION_POM
  !
  ! !USES:
  !
  use GO,            only : gol, goErr, goPr
  use tm5_distgrid,  only : dgrid, get_distgrid, scatter, gather
  use dims,          only : nregions, okdebug
  use global_types,  only : emis_data, d3_data
  use emission_data, only : emis_input_dir_aerocom
  use emission_data, only : emis_input_year_oc
  use emission_read, only : used_providers_aer, has_aer_emis 

  implicit none
  private
  !
  ! !PUBLIC MEMBER FUNCTIONS:
  !
  public   :: emission_pom_init     ! allocate
  public   :: emission_pom_done     ! deallocate
  public   :: emission_pom_declare  ! read data
  !
  ! !PRIVATE DATA MEMBERS:
  !
  character(len=*), parameter  ::  mname = 'emission_pom'

  type(emis_data), dimension(:,:), allocatable :: pom_emis_2d, pom_bf_emis_2d
  type(d3_data),   dimension(:,:), allocatable :: pom_emis_3d

  integer                           :: pom_2dsec, pom_3dsec
  !
  ! !REVISION HISTORY: 
  !    ? ??? 2005 - Elisabetta Vignati - changed for coupling with M7
  !    1 Sep 2010 - Achim Strunk - introduced AR5 emissions
  !                              - reorganised the array structures and 
  !                                vertical distribution facility
  !                              - cleaning
  !   26 Jun 2012 - P. Le Sager  - adapted for lon-lat MPI domain decomposition
  !
  ! !REMARKS:
  !
  !EOP
  !------------------------------------------------------------------------

CONTAINS

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    EMISSION_POM_INIT
  !
  ! !DESCRIPTION: Allocate space needed to handle the emissions. 
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE EMISSION_POM_INIT( status )
    !
    ! !USES:
    !
    use dims,          only : lm
    use emission_read, only : providers_def, numb_providers
    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 - v0
    !   26 Jun 2012 - P. Le Sager  - adapted for lon-lat MPI domain decomposition 
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

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

    status = 0
    if(.not. has_aer_emis) return

    ! nb of sectors
    pom_2dsec = 0
    pom_3dsec = 0
    do lprov = 1, numb_providers
       if (count(used_providers_aer.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
             pom_2dsec = pom_2dsec + n_ar5_ant_sec*count('OC'.eq.ar5_cat_ant) + &
                                     n_ar5_shp_sec*count('OC'.eq.ar5_cat_shp)
         
             if (LAR5BMB) pom_2dsec = pom_2dsec + n_ar5_bmb_sec*count('OC'.eq.ar5_cat_bmb)

             pom_3dsec = pom_3dsec + count('OC'.eq.ar5_cat_air)
          else
             pom_2dsec = pom_2dsec + providers_def(lprov)%nsect2d
             pom_3dsec = pom_3dsec + providers_def(lprov)%nsect3d
          endif
       endif
    enddo

    allocate( pom_emis_2d( nregions, pom_2dsec ) )
    allocate( pom_bf_emis_2d( nregions, pom_2dsec ) )
    allocate( pom_emis_3d( nregions, pom_3dsec ) )

    ! 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,pom_2dsec
          allocate( pom_emis_2d(region,lsec)%surf(i1:i2,j1:j2) )
          allocate( pom_bf_emis_2d(region,lsec)%surf(i1:i2,j1:j2) )
       end do
       
       do lsec=1,pom_3dsec
          allocate( pom_emis_3d(region,lsec)%d3(i1:i2,j1:j2,lmr) )
       end do

    enddo

    ! ok
    status = 0

  END SUBROUTINE EMISSION_POM_INIT
  !EOC

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

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

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

    status = 0
    if(.not. has_aer_emis) return

    do region = 1, nregions
       do lsec=1,pom_2dsec
          deallocate( pom_emis_2d(region,lsec)%surf )
          deallocate( pom_bf_emis_2d(region,lsec)%surf )
       end do
       do lsec=1,pom_3dsec
          deallocate( pom_emis_3d(region,lsec)%d3 )
       end do
    end do

    deallocate( pom_emis_2d )
    deallocate( pom_bf_emis_2d )
    deallocate( pom_emis_3d )

    ! ok
    status = 0

  END SUBROUTINE EMISSION_POM_DONE
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    EMISSION_POM_DECLARE
  !
  ! !DESCRIPTION: Opens, reads and evaluates input files (per month). 
  !               Provides emissions on 2d/3d-arrays which are then given 
  !               to emis_number and emis_mass, which are used in 
  !               sedimentation. (no *_apply!)
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE EMISSION_POM_DECLARE( status )
    !
    ! !USES:
    !
    use binas,         only : pi
    use partools,      only : isRoot, par_broadcast
    use toolbox,       only : coarsen_emission, distribute_emis2D
    use dims,          only : im, jm, lm, idate, sec_month, nlon360, nlat180
    use chem_param,    only : xmc, sigma_lognormal, pom_density
    use chem_param,    only : mode_aii, mode_ais, mode_acs
    use MDF,           only : MDF_Open,  MDF_NETCDF, MDF_Get_Var
    use MDF,           only : MDF_Close, MDF_READ, MDF_Inq_VarID
    use emission_data, only : emis_mass, emis_number
    use emission_data, only : rad_emi_ff_sol,  rad_emi_ff_insol
    use emission_data, only : rad_emi_ene_sol, rad_emi_ene_insol
    use emission_data, only : rad_emi_ind_sol, rad_emi_ind_insol
    use emission_data, only : rad_emi_tra_sol, rad_emi_tra_insol
    use emission_data, only : rad_emi_shp_sol, rad_emi_shp_insol
    use emission_data, only : rad_emi_air_sol, rad_emi_air_insol
    use emission_data, only : rad_emi_bf_sol,  rad_emi_bf_insol
    use emission_data, only : rad_emi_bb_sol,  rad_emi_bb_insol
    !use emission_data, only : rad_soa
    use emission_data, only : frac_pom_sol_ff, frac_pom_sol_bf, frac_pom_sol_bb, frac_soa_sol
    !use emission_data, only : oc2pom   
    use emission_data, only : oc2pom_ff, oc2pom_bf, oc2pom_bb, oc2pom_soa
    use emission_data, only : msg_emis, emis_temp, emission_vdist_by_sector
    use emission_data, only : vd_class_name_len
    ! ---------------- CMIP6 - AR5 - GFED - RETRO - MACC --------------------
    use emission_data, only : LAR5BMB
    use emission_data, only : emis_input_dir_mac
    use emission_data, only : emis_input_dir_retro
    use emission_data, only : emis_input_dir_gfed

    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_gfed_ReadSector
    use emission_read, only : emission_retro_ReadSector
 
    use emission_read, only : sector_name_len
    use emission_read, only : sectors_def, numb_sectors
    use emission_read, only : ar5_dim_3ddata
    !RM
    use mo_aero,       only : nsoa
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)        ::  status
    !
    ! !REVISION HISTORY: 
    !    1 Oct 2010 - Achim Strunk - revamped for AR5
    !   26 Jun 2012 - P. Le Sager  - adapted for lon-lat MPI domain decomposition
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter  :: rname = mname//'/emission_pom_declare'
    integer                      :: region, hasData
    integer, parameter           :: add_field=0
    integer, parameter           :: amonth=2
    !real                         :: rad_aver_mass
    real                         :: numbscale_exp
    integer                      :: imr, jmr, lmr, lsec, i1, i2, j1, j2, l

    real(kind=4), dimension(:,:,:  ), allocatable   :: hdfr3

    integer                             :: fid, varid  ! hdf related

    ! ---------------------------------------------------------------

    real, dimension(:,:),   allocatable :: field2d
    real, dimension(:,:,:), allocatable :: field3d, field3d2
    real, dimension(:,:),   allocatable :: field2d_bf
    integer                             :: mode_sol
    real                                :: mass2numb_fact
    real                                :: mass2numb_ff_sol,    mass2numb_ff_insol
    real                                :: mass2numb_ene_sol,   mass2numb_ene_insol
    real                                :: mass2numb_ind_sol,   mass2numb_ind_insol
    real                                :: mass2numb_tra_sol,   mass2numb_tra_insol
    real                                :: mass2numb_shp_sol,   mass2numb_shp_insol
    real                                :: mass2numb_air_sol,   mass2numb_air_insol
    real                                :: mass2numb_bf_sol,    mass2numb_bf_insol
    real                                :: mass2numb_bb_sol,    mass2numb_bb_insol
    real                                :: mass2numb_nonbf_sol, mass2numb_nonbf_insol
    real                                :: oc2pom
    type(d3_data),  dimension(nregions) :: emis3d, work, work_bf, work3d
    type(d3_data),  dimension(nregions) :: frac_bf
    type(emis_data),dimension(nregions) :: wrk2D, wrk2D_bf, emis_glb
    integer                             :: seccount2d,   seccount3d
    character(len=sector_name_len)      :: tmpsector
    character(len=vd_class_name_len)    :: tmpvsplit

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

    status = 0
    if(.not. has_aer_emis) return
    
    write(gol,'(" EMISS-INFO ------------- read POM emissions -------------")'); call goPr
    
    ! Reset emissions
    do region = 1, nregions
       do lsec=1,pom_2dsec
          pom_emis_2d(region,lsec)%surf = 0.0
          pom_bf_emis_2d(region,lsec)%surf = 0.0
       end do
       do lsec=1,pom_3dsec
          pom_emis_3d(region,lsec)%d3   = 0.0
       end do
    end do

    ! Allocate work arrays
    do region = 1, nregions

       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
       allocate( frac_bf(region)%d3 (i1:i2,j1:j2, 1             ) )
       
    end do

    ! Global arrays for coarsening    
    do region = 1, nregions

       if (isRoot)then
          allocate(work(region)%d3(im(region),jm(region),ar5_dim_3ddata))
          allocate(work_bf(region)%d3(im(region),jm(region),1))
       else
          allocate(work(region)%d3(1,1,1))
          allocate(work_bf(region)%d3(1,1,1)) 
       end if
    enddo

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


    ! -----------------------------
    ! get emissions (ATTENTION: THIS IS Organic Carbon! )
    write(gol,'(1x,80("-"))') ; call goPr
    write(gol,*) '    WARNING: OC emissions are used instead of POM !' ; call goPr
    !write(gol,*) '             masses are transformed using constant factor oc2pom' ; call goPr
    write(gol,*) '             masses are transformed using constant factors' ; call goPr
    write(gol,*) '             for vegetation fires and other sources' ; call goPr
    write(gol,'(1x,80("-"))') ; call goPr


    ! --------------------------------
    ! 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
       allocate( field2d_bf( nlon360, nlat180 ) ) ; field2d_bf = 0.0
    else
       allocate( field3d( 1, 1, 1 ) )
       allocate( field2d_bf( 1, 1 ) )
    end if

    ! mass to number conversion factors for the relevant modes
    numbscale_exp  = EXP(1.5*(LOG(sigma_lognormal(mode_aii)))**2)
    mass2numb_fact = 3./(4.*pi*(numbscale_exp**3)*pom_density)
    mass2numb_ff_insol  = mass2numb_fact/(rad_emi_ff_insol**3)
    mass2numb_ene_insol = mass2numb_fact/(rad_emi_ene_insol**3)
    mass2numb_ind_insol = mass2numb_fact/(rad_emi_ind_insol**3)
    mass2numb_tra_insol = mass2numb_fact/(rad_emi_tra_insol**3)
    mass2numb_shp_insol = mass2numb_fact/(rad_emi_shp_insol**3)
    mass2numb_air_insol = mass2numb_fact/(rad_emi_air_insol**3)
    mass2numb_bf_insol  = mass2numb_fact/(rad_emi_bf_insol**3)
    mass2numb_bb_insol  = mass2numb_fact/(rad_emi_bb_insol**3)

    numbscale_exp  = EXP(1.5*(LOG(sigma_lognormal(mode_ais)))**2)
    mass2numb_fact = 3./(4.*pi*(numbscale_exp**3)*pom_density)
    mass2numb_ff_sol  = mass2numb_fact/(rad_emi_ff_sol**3)
    mass2numb_ene_sol = mass2numb_fact/(rad_emi_ene_sol**3)
    mass2numb_ind_sol = mass2numb_fact/(rad_emi_ind_sol**3)
    mass2numb_tra_sol = mass2numb_fact/(rad_emi_tra_sol**3)
    mass2numb_shp_sol = mass2numb_fact/(rad_emi_shp_sol**3)
    mass2numb_air_sol = mass2numb_fact/(rad_emi_air_sol**3)

    numbscale_exp  = EXP(1.5*(LOG(sigma_lognormal(mode_acs)))**2)
    mass2numb_fact = 3./(4.*pi*(numbscale_exp**3)*pom_density)
    mass2numb_bf_sol = mass2numb_fact/(rad_emi_bf_sol**3)
    mass2numb_bb_sol = mass2numb_fact/(rad_emi_bb_sol**3) 
    
    sec : do lsec = 1, numb_sectors

       if (count(used_providers_aer.eq.sectors_def(lsec)%prov).eq.0) cycle 
       
       if (associated(sectors_def(lsec)%species)) then
          if (count('OC'.eq.sectors_def(lsec)%species).eq.0) cycle
          if ((trim(sectors_def(lsec)%catname) .eq. 'biomassburning').and.(.not.LAR5BMB)) cycle
       endif

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

       if( trim(sectors_def(lsec)%catname) == 'biomassburning' ) then
         oc2pom = oc2pom_bb
       else
         oc2pom = oc2pom_ff
       endif

       if (isRoot) then ! READ

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

          case( 'CMIP6' )
            call emission_cmip6_ReadSector( 'OC', emis_input_year_oc, idate(2), lsec, field3d, status, field2d_bf )
            IF_NOTOK_RETURN(status=1;deallocate(field3d,field2d_bf))

          case( 'CMIP6BMB' )
            call emission_cmip6bmb_ReadSector( 'OC', emis_input_year_oc, idate(2), lsec, field3d, status )
            IF_NOTOK_RETURN(status=1;deallocate(field3d))

          case( 'AR5' )

             ! Screen out agricultural and solvent sectors for OC, 
             ! 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( 'OC', emis_input_year_oc, idate(2), lsec, field3d, status )
               IF_NOTOK_RETURN(status=1;deallocate(field3d))
             endif

          case( 'MACC' )

             ! Screen out 'soil', 'nat', 'oc', bio', and 'air' since they are not available for OC.
             if ( ( .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_oc') ) .and. &
                  ( .not. (trim(sectors_def(lsec)%name) .eq. 'emiss_bio' ) ) .and. &
                  ( .not. (trim(sectors_def(lsec)%name) .eq. 'emiss_air') ) ) then

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

             end if

          case('GFEDv3')
             call emission_gfed_ReadSector( emis_input_dir_gfed, 'oc', emis_input_year_oc, idate(2), &
                                            'GFED', 'kg / s', field3d(:,:,1), status )
             IF_NOTOK_RETURN(status=1)

          case('RETRO')
             call emission_retro_ReadSector( emis_input_dir_retro, 'OC', emis_input_year_oc, idate(2), &
                                             sectors_def(lsec)%name, 'kg / s', field3d(:,:,1), status )
             IF_NOTOK_RETURN(status=1)

          case default
             write(gol,*) "Error in list of providers for POM"; 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 POM 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 POM emissions for ",a," ",a," for month ",i2 )') &
                     trim(sectors_def(lsec)%prov), trim(sectors_def(lsec)%name), idate(2) ; call goPr
             endif

             ! convert from OC to POM
             if( trim(sectors_def(lsec)%catname) == 'biomassburning' ) then
                field3d = field3d * oc2pom_bb
             else
                field3d = field3d * oc2pom_ff
                ! correct biofuel contribution in surface level
                field3d(:,:,1) = field3d(:,:,1) + field2d_bf(:,:) * (oc2pom_bf - oc2pom_ff)
                field2d_bf = field2d_bf * oc2pom_bf
             endif

             ! scale from kg/s to kg/month
             field3d = field3d * sec_month ! kg / month 
             field2d_bf = field2d_bf * sec_month
             hasData=1
          end if

       end if

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

       if (hasData == 0) cycle sec

       ! distinguish 2d/3d sectors
       if( sectors_def(lsec)%f3d ) 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, 'POM', oc2pom_ff*xmc, sum(field3d) )

             ! distribute to work arrays in regions
             call Coarsen_Emission( 'POM '//trim(sectors_def(lsec)%name), nlon360,  nlat180, ar5_dim_3ddata,   &
                                     field3d,  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 )

             pom_emis_3d(region,seccount3d)%d3 = pom_emis_3d(region,seccount3d)%d3 + emis3d(region)%d3

          end do

          ! add to emis target arrays
          do region = 1, nregions
             emis_mass  (region,mode_aii)%d4(:,:,:,2) = &
             emis_mass  (region,mode_aii)%d4(:,:,:,2) + &
               (1.-frac_pom_sol_ff) * pom_emis_3d(region,seccount3d)%d3(:,:,:)
             emis_mass  (region,mode_ais)%d4(:,:,:,3) = &
             emis_mass  (region,mode_ais)%d4(:,:,:,3) + &
               (   frac_pom_sol_ff) * pom_emis_3d(region,seccount3d)%d3(:,:,:)
             
             emis_number(region,mode_aii)%d4(:,:,:,2) = &
             emis_number(region,mode_aii)%d4(:,:,:,2) + &
               (1.-frac_pom_sol_ff) * pom_emis_3d(region,seccount3d)%d3(:,:,:) * mass2numb_air_insol
             emis_number(region,mode_ais)%d4(:,:,:,3) = &
             emis_number(region,mode_ais)%d4(:,:,:,3) + &
               (   frac_pom_sol_ff) * pom_emis_3d(region,seccount3d)%d3(:,:,:) * mass2numb_air_sol
          end do

       else ! 2d sector

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

             ! in the following print statement
             ! a possible difference between oc2pom_bf and oc2pom_ff 
             ! is not accounted for,
             ! but this has no further consequences
             call msg_emis( amonth, trim(sectors_def(lsec)%prov), sectors_def(lsec)%name, 'POM', oc2pom*xmc, &
                            sum(field3d(:,:,1)) )
             IF_NOTOK_RETURN(status=1)
             
             call coarsen_emission( 'POM '//sectors_def(lsec)%name, nlon360, nlat180, field3d(:,:,1), wrk2D, add_field, status)
             IF_NOTOK_RETURN(status=1)

             call coarsen_emission( 'POM solid biofuel '//sectors_def(lsec)%name, nlon360, nlat180, field2d_bf(:,:), wrk2D_bf, add_field, status)
             IF_NOTOK_RETURN(status=1)

          end if
             
          ! get temporary array for vertical distribution
          do region = 1, nregions

             call scatter(dgrid(region), pom_emis_2d(region,seccount2d)%surf, work(region)%d3(:,:,1), 0, status)
             IF_NOTOK_RETURN(status=1)

             call scatter(dgrid(region), pom_bf_emis_2d(region,seccount2d)%surf, work_bf(region)%d3(:,:,1), 0, status)
             IF_NOTOK_RETURN(status=1)
             
             emis3d(region)%d3 = 0.0

             ! do vertical distribution
             call emission_vdist_by_sector( sectors_def(lsec)%vdisttype, 'POM', region, pom_emis_2d(region,seccount2d), &
             emis3d(region), status )
             IF_NOTOK_RETURN(status=1)

             if( trim(sectors_def(lsec)%catname) == 'biomassburning' ) then

                ! add to emis target arrays
                emis_mass  (region,mode_aii)%d4(:,:,:,2) = &
                emis_mass  (region,mode_aii)%d4(:,:,:,2) + emis3d(region)%d3(:,:,:) * &
                  (1.-frac_pom_sol_bb)

                emis_number(region,mode_aii)%d4(:,:,:,2) = &
                emis_number(region,mode_aii)%d4(:,:,:,2) + emis3d(region)%d3(:,:,:) * &
                  (1.-frac_pom_sol_bb) * mass2numb_bb_insol 

                emis_mass  (region,mode_acs)%d4(:,:,:,3) = &
                emis_mass  (region,mode_acs)%d4(:,:,:,3) + emis3d(region)%d3(:,:,:) * &
                  frac_pom_sol_bb 

                emis_number(region,mode_acs)%d4(:,:,:,3) = &
                emis_number(region,mode_acs)%d4(:,:,:,3) + emis3d(region)%d3(:,:,:) * &
                  frac_pom_sol_bb * mass2numb_bb_sol
               
             else 

                ! currenly only for CMIP6 sector names
                select case( trim(sectors_def(lsec)%name) )
                   case ('ENE')
                      mass2numb_nonbf_sol   = mass2numb_ene_sol
                      mass2numb_nonbf_insol = mass2numb_ene_insol
                   case ('IND')
                      mass2numb_nonbf_sol   = mass2numb_ind_sol
                      mass2numb_nonbf_insol = mass2numb_ind_insol
                   case ('TRA')
                      mass2numb_nonbf_sol   = mass2numb_tra_sol
                      mass2numb_nonbf_insol = mass2numb_tra_insol
                   case ('SHP')
                      mass2numb_nonbf_sol   = mass2numb_shp_sol
                      mass2numb_nonbf_insol = mass2numb_shp_insol 
                   case default 
                      mass2numb_nonbf_sol   = mass2numb_ff_sol
                      mass2numb_nonbf_insol = mass2numb_ff_insol
                end select

                if ( trim(sectors_def(lsec)%prov) == 'CMIP6' ) then

                   ! calculate mass fraction related to solid biofuel
                   where ( pom_emis_2d(region,seccount2d)%surf(:,:) > tiny(0.0) )
                      frac_bf(region)%d3(:,:,1) = pom_bf_emis_2d(region,seccount2d)%surf(:,:) / &
                                                  pom_emis_2d(region,seccount2d)%surf(:,:)
                   elsewhere
                      frac_bf(region)%d3(:,:,1) = 0.0
                   endwhere

                   ! for safety, prevent fractions larger than unity.
                   where (frac_bf(region)%d3(:,:,1) > 1.0 )
                      frac_bf(region)%d3(:,:,1) = 1.0
                   endwhere
                else
                   frac_bf(region)%d3(:,:,1) = 0.0
                endif

                do l = 1, lmr
                   emis_mass  (region,mode_aii)%d4(:,:,l,2) = &
                   emis_mass  (region,mode_aii)%d4(:,:,l,2) + emis3d(region)%d3(:,:,l) * &
                     ( (1.-frac_bf(region)%d3(:,:,1)) * (1.-frac_pom_sol_ff) + &
                           frac_bf(region)%d3(:,:,1)  * (1.-frac_pom_sol_bf) ) 

                   emis_number(region,mode_aii)%d4(:,:,l,2) = &
                   emis_number(region,mode_aii)%d4(:,:,l,2) + emis3d(region)%d3(:,:,l) * &
                     ( (1.-frac_bf(region)%d3(:,:,1)) * (1.-frac_pom_sol_ff) * mass2numb_nonbf_insol + &
                           frac_bf(region)%d3(:,:,1)  * (1.-frac_pom_sol_bf) * mass2numb_bf_insol )

                   emis_mass  (region,mode_ais)%d4(:,:,l,3) = &
                   emis_mass  (region,mode_ais)%d4(:,:,l,3) + emis3d(region)%d3(:,:,l) * &
                     ( (1.-frac_bf(region)%d3(:,:,1)) * frac_pom_sol_ff )

                   emis_number(region,mode_ais)%d4(:,:,l,3) = &
                   emis_number(region,mode_ais)%d4(:,:,l,3) + emis3d(region)%d3(:,:,l) * &
                     ( (1.-frac_bf(region)%d3(:,:,1)) * frac_pom_sol_ff * mass2numb_nonbf_sol )

                   emis_mass  (region,mode_acs)%d4(:,:,l,3) = &
                   emis_mass  (region,mode_acs)%d4(:,:,l,3) + emis3d(region)%d3(:,:,l) * &
                     (     frac_bf(region)%d3(:,:,1)  * frac_pom_sol_bf )

                   emis_number(region,mode_acs)%d4(:,:,l,3) = &
                   emis_number(region,mode_acs)%d4(:,:,l,3) + emis3d(region)%d3(:,:,l) * &
                     (     frac_bf(region)%d3(:,:,1)  * frac_pom_sol_bf * mass2numb_bf_sol )
                enddo

             endif
               
          end do

       end if

    end do sec ! sectors

    deallocate( field3d )
    deallocate( field2d_bf )
    do region = 1, nregions
       if (associated(wrk2D(region)%surf)) nullify(wrk2D(region)%surf)
       if (associated(wrk2D_bf(region)%surf)) nullify(wrk2D_bf(region)%surf)
       deallocate( work(region)%d3 )
       deallocate( work_bf(region)%d3 )       
       deallocate( work3d(region)%d3 ) 
       deallocate( frac_bf(region)%d3 )
    end do

    ! check sectors found
    if( seccount2d /= pom_2dsec ) then 
       write(gol,'(80("-"))') ; call goPr
       write(gol,'("ERROR: 2d sectors do not equal total number:",i4," /= ",i4," !")') seccount2d, pom_2dsec ; call goErr
       write(gol,'(80("-"))') ; call goPr
       status=1; return
    end if
    if( seccount3d /= pom_3dsec ) then 
       write(gol,'(80("-"))') ; call goPr
       write(gol,'("ERROR: 3d sectors do not equal total number:",i4," /= ",i4," !")') seccount3d, pom_3dsec ; call goErr
       write(gol,'(80("-"))') ; call goPr
       status=1; return
    end if


    ! Old AEROCOM SOA emissions available only when nsoa == 0

    ! Otherwise the production of SOA is calculated from isoprene and monoterpene 
    ! through chemical reactions in ebischeme.f90.

    if (nsoa .EQ. 0) then

       ! ---------------------------------------------------
       ! SOA condenses onto existing particles in the Aitken soluble mode
       ! ---------------------------------------------------
       
       !    select case( emis_input_provider )
       
       !    case( 'AR5','MACC' )
       write(gol,'(80("-"))') ; call goPr
       write(gol,*) '    WARNING: SOA emissions are employed in boxes ' ; call goPr
       write(gol,*) '             with too few particles in mode_ais! ' ; call goPr
       write(gol,*) '             Masses of SOA are neglected !!! ' ; call goPr
       write(gol,'(80("-"))') ; call goPr
       !    end select
       
       tmpsector = 'anthropogenic'
       ! assume SOA emissions to be vertically split due to this class:
       tmpvsplit = 'nearsurface'
       
       ! work arrays
       if (isRoot) then
          do region = 1, nregions
             allocate(emis_glb(region)%surf(im(region), jm(region)))
          end do
       else
          do region = 1, nregions
             allocate(emis_glb(region)%surf(1,1))
          end do
       end if
       
       ! Read
       if (isRoot) then
          
          allocate( hdfr3(nlon360,nlat180,1) ) ; hdfr3   = 0.0
          allocate( field2d(nlon360,nlat180) ) ; field2d = 0.0

          ! FIXME - this is untested, SOA.nc4 needs to be created - but not needed in EC-Earth
          call MDF_Open( trim(emis_input_dir_aerocom)//'/SOA.nc4', MDF_NETCDF, MDF_READ, fid, status)
          IF_NOTOK_RETURN(status=1)
          
          call MDF_Inq_VarID( fid, 'FIELD', varid, status)
          IF_NOTOK_RETURN(status=1)
          
          call MDF_Get_Var( fid, varid, hdfr3, status, start=(/1,1,idate(2)/), count=(/nlon360,nlat180,1/))
          IF_NOTOK_RETURN(status=1)
          
          CALL MDF_Close( fid, status )
          IF_NOTOK_RETURN(status=1)       
          
          ! paste to 2d field
          field2d = hdfr3(:,:,1)
          
          !>>> TvN
          ! SOA emissions are provided in Tg POM, and are based on 
          ! the assumption that 15% of the emitted terpene mass is converted to SOA
          ! (Kanakidou et al., ACP, 2005; Dentener et al., ACP, 2006).
          ! However, it seems that such an assumption
          ! is inconsistent with a POM to OC ratio of about 2.0 for freshly formed SOA
          ! (Aiken et al., Environ. Sci. Technol., 2008).
          ! The VOCs with the largest potential for SOA formation 
          ! are mono-terpenes (C10H16),  which account for 40-80%
          ! of the overall terpene emissions, and seqsquiterpenes (C15H24),
          ! which have a 100% SOA yield.
          ! The full weight to carbon weight ratios for these components
          ! are 1.133 and 1.16, respectively.
          ! Therefore, it seems reasonable to increase the SOA emissions
          ! as provided in the AeroCom input file 
          ! by a factor of about ocpom_soa/1.15.
          ! Increasing the SOA emissions from AeroCom,
          ! which amount to 19.2 Tg/yr, is also consistent with AR5,
          ! where 20 Tg/yr is assumed to be a lower bound.
          field2d = field2d * oc2pom_soa/1.15
          !<<< TvN
          
          deallocate(hdfr3)
          
          ! prompt some numbers
          call msg_emis( amonth, 'AEROCOM', 'SOA-'//trim(tmpsector), 'POM', oc2pom_soa*xmc, sum(field2d) )
          IF_NOTOK_RETURN(status=1;deallocate(field2d))
          
          ! initialise/reset emis_temp(regions)
          do region = 1, nregions
             emis_temp(region)%surf = 0.0
          end do
          
          ! regrid from field2d to emis_glb
          call coarsen_emission( 'SOA '//trim(tmpsector), nlon360, nlat180, field2d, emis_glb, add_field, status)
          IF_NOTOK_RETURN(status=1;deallocate(field2d))
          
          deallocate( field2d )

       end if
       
       !>>> TvN
       ! EV total particle number emitted in a gridbox
       !rad_aver_mass = rad_soa*EXP(1.5*(LOG(sigma_lognormal(mode_ais)))**2) 
       !mass_to_numb  = 3./(4.*pi*(rad_aver_mass**3)*pom_density)
       
       ! Scatter, vertical distribution
       ! ------------------------------       
       do region = 1, nregions
          
          call scatter(dgrid(region), emis_temp(region)%surf, emis_glb(region)%surf, 0, status)
          IF_NOTOK_RETURN(status=1)
          
          emis3d(region)%d3 = 0.0
          
          call emission_vdist_by_sector( tmpvsplit, 'SOA', region, emis_temp(region), emis3d(region), status )
          IF_NOTOK_RETURN(status=1)
          
          ! --------------
          ! Comment by EV:
          ! Let the SOA condense on the existing aerosols. Not emit the aerosol numbers.
          ! Only if there are insufficient aerosols (<1.0), we emit the aerosol numbers.
          
          ! add masses
          !emis_mass(region,mode_ais)%d4(:,:,:,3) = &
          !     emis_mass(region,mode_ais)%d4(:,:,:,3) + emis3d(region)%d3(:,:,:) 
          
          ! TB use the new compounds available to track the emissions from SOA.hdf
          ! emis_mass(region,mode_ais)%D4(:,:,:,4) put mass as soa-compound
          ! emis_mass(region,mode_aii)%D4(:,:,:,3) put mass as soa-compound 
          emis_mass(region,mode_ais)%d4(:,:,:,4) = &
               emis_mass(region,mode_ais)%d4(:,:,:,4) + &
               frac_soa_sol * emis3d(region)%d3(:,:,:)
          emis_mass(region,mode_aii)%d4(:,:,:,3) = &
               emis_mass(region,mode_aii)%d4(:,:,:,3) + &
               (1.-frac_soa_sol) * emis3d(region)%d3(:,:,:)


       ! add numbers where appropriate
!       where(emis_number(region,mode_ais)%d4(:,:,:,3) < 1.)
          !emis_number(region,mode_ais)%d4(:,:,:,3) = &
          !     emis_number(region,mode_ais)%d4(:,:,:,3) + emis3d(region)%d3(:,:,:) * mass_to_numb
!          emis_number(region,mode_ais)%d4(:,:,:,3) = &
!          emis_number(region,mode_ais)%d4(:,:,:,3) + &
!               frac_soa_sol * emis3d(region)%d3(:,:,:) * mass_to_numb
!       endwhere
       
!       where(emis_number(region,mode_aii)%d4(:,:,:,2) < 1.)
!          emis_number(region,mode_aii)%d4(:,:,:,2) = &
!          emis_number(region,mode_aii)%d4(:,:,:,2) + &
!               (1.-frac_soa_sol) * emis3d(region)%d3(:,:,:) * mass_to_numb
!       endwhere
!<<< TvN

       enddo
   
       ! Done
       do region = 1, nregions
          deallocate(emis_glb(region)%surf)
          deallocate(  emis3d(region)%d3  )
       end do
    end if
    status = 0
    
  END SUBROUTINE EMISSION_POM_DECLARE
  !EOC

END MODULE EMISSION_POM