ecearth3/sources/tm5mp/proj/cbm4/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 fuel emissions, vegetation fire emission, 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 emission radii have been modified
! following the AeroCom recommendations of Dentener et al. (ACP, 2006),
! adapted for the M7 modes by Stier et al. (JGR, 2005).
! The soluble fraction of the emissions from vegetation fires are now emitted
! in the accumulation mode (see mo_aero.F90).
! For the insoluble fraction, the mean emission radius is
! is assumed to be the same as for fossil-fuel emissions
! (as in Stier et al.)
!
! Furthermore, POM from SOA is considered 65% soluble,
! as recommended by AeroCom.
! The above mentioned paper by Kanakidou et al. 
! does not seem to support 100% solubility.
!
! 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
  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
  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_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) )
       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 )
       end do
       do lsec=1,pom_3dsec
          deallocate( pom_emis_3d(region,lsec)%d3 )
       end do
    end do

    deallocate( pom_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 chem_param,    only : rad_emi_ff, rad_emi_vg, rad_emi_vg_insol
    !use chem_param,    only : rad_soa
    use chem_param,    only : frac_pom_sol_ff, frac_pom_sol_vg, frac_soa_sol
    use MDF,           only : MDF_Open,  MDF_HDF4, 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 : oc2pom   
    use emission_data, only : oc2pom_ff, oc2pom_vg, oc2pom_soa
    use emission_data, only : msg_emis, emis_temp, emission_vdist_by_sector
    use emission_data, only : vd_class_name_len
    ! ---------------- 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_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
    !
    ! !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, numbscale
    integer                      :: imr, jmr, lmr, lsec, i1, i2, j1, j2

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

    integer                             :: fid, varid  ! hdf related

    ! ---------------------------------------------------------------
    ! AR5
    real, dimension(:,:),   allocatable :: field2d
    real, dimension(:,:,:), allocatable :: field3d,      field3d2
    integer                             :: mode_sol
    real                                :: rad_emi_sol, rad_emi_insol
    real                                :: mass2numb_sol, mass2numb_insol
    real                                :: fracpomsol
    real                                :: oc2pom
    type(d3_data),  dimension(nregions) :: emis3d,       work,       work3d
    type(emis_data),dimension(nregions) :: wrk2D,        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
    
    ! EV pom emitted in the insoluble aitken mode

    ! Reset emissions
    do region = 1, nregions
       do lsec=1,pom_2dsec
          pom_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
    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


    ! -----------------------------
    ! 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
    else
       allocate( field3d( 1, 1, 1 ) )
    end if

    
    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
       
       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_vg
       else
         oc2pom = oc2pom_ff
       endif

       if (isRoot) then ! READ

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

          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, idate(2), lsec, field3d, status )
               IF_NOTOK_RETURN(status=1;deallocate(field3d))
             endif

          case( 'MACC' )

             ! Screen out 'soil', 'nat', 'oc', bio', 'oc', and 'air' since they are not available for BC.
             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, 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, 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, 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 using oc2pom (set in emission_data)
             field3d = field3d * oc2pom

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

       end if

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

       if (hasData == 0) cycle sec

       ! special values for biomass burning emissions
       if( trim(sectors_def(lsec)%catname) == 'biomassburning' ) then 
          mode_sol = mode_acs
          rad_emi_sol = rad_emi_vg
          rad_emi_insol = rad_emi_vg_insol
          fracpomsol = frac_pom_sol_vg
       else
          mode_sol = mode_ais
          rad_emi_sol = rad_emi_ff
          rad_emi_insol = rad_emi_ff
          fracpomsol = frac_pom_sol_ff
       end if

       ! mass to number factors for the soluble and insoluble Aitken modes
       ! and soluble accumulation modes
       numbscale       = rad_emi_insol*EXP(1.5*(LOG(sigma_lognormal(mode_aii)))**2) 
       mass2numb_insol = 3./(4.*pi*(numbscale**3)*pom_density)
       numbscale      = rad_emi_sol*EXP(1.5*(LOG(sigma_lognormal(mode_sol)))**2)
       mass2numb_sol  = 3./(4.*pi*(numbscale**3)*pom_density)
       

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

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

!DBG          
!           ! up to now only aircraft (it should be this)
!           if( trim(sectors_def(lsec)%name) /= 'emiss_air' ) then 
!              write(gol,'(80("-"))') ; call goPr
!              write(gol,'("ERROR: 3d sector `",a,"` not available !")') sectors_def(lsec)%name ; call goErr
!              write(gol,'(80("-"))') ; call goPr
!              deallocate( field3d )
!              status=1; return
!           end if
          
          if (isRoot) then

             ! write some numbers
             call msg_emis( amonth, trim(sectors_def(lsec)%prov), sectors_def(lsec)%name, 'POM', oc2pom*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.-fracpomsol) * pom_emis_3d(region,seccount3d)%d3(:,:,:)
             emis_mass  (region,mode_ais)%d4(:,:,:,3) = &
               emis_mass  (region,mode_ais)%d4(:,:,:,3) + &
               (   fracpomsol) * pom_emis_3d(region,seccount3d)%d3(:,:,:)
             
             emis_number(region,mode_aii)%d4(:,:,:,2) = &
               emis_number(region,mode_aii)%d4(:,:,:,2) + &
               (1.-fracpomsol) * pom_emis_3d(region,seccount3d)%d3(:,:,:) * mass2numb_insol
             emis_number(region,mode_ais)%d4(:,:,:,3) = &
               emis_number(region,mode_ais)%d4(:,:,:,3) + &
               (   fracpomsol) * pom_emis_3d(region,seccount3d)%d3(:,:,:) * mass2numb_sol
          end do

       else ! 2d sector

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

             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)

          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)
             
             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)

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

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

             ! soluble emissions from vegetation fires to accumulation mode,
             ! other soluble emissions to Aitken mode
             if( trim(sectors_def(lsec)%catname) == 'biomassburning' ) then
               emis_mass  (region,mode_acs)%d4(:,:,:,3) = &
               emis_mass  (region,mode_acs)%d4(:,:,:,3) + &
                 (   fracpomsol) * emis3d(region)%d3(:,:,:)

               emis_number(region,mode_acs)%d4(:,:,:,3) = &
               emis_number(region,mode_acs)%d4(:,:,:,3) + &
                 (   fracpomsol) * emis3d(region)%d3(:,:,:) * mass2numb_sol
             else
               emis_mass  (region,mode_ais)%d4(:,:,:,3) = &
               emis_mass  (region,mode_ais)%d4(:,:,:,3) + &
                 (   fracpomsol) * emis3d(region)%d3(:,:,:) 

               emis_number(region,mode_ais)%d4(:,:,:,3) = &
               emis_number(region,mode_ais)%d4(:,:,:,3) + &
                 (   fracpomsol) * emis3d(region)%d3(:,:,:) * mass2numb_sol
             endif

          end do

       end if

    end do sec ! sectors

    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 /= 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


    ! ---------------------------------------------------
    ! 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

       call MDF_Open( trim(emis_input_dir_aerocom)//'/SOA.hdf', MDF_HDF4, 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(:,:,:) 
       emis_mass(region,mode_ais)%d4(:,:,:,3) = &
       emis_mass(region,mode_ais)%d4(:,:,:,3) + &
            frac_soa_sol * emis3d(region)%d3(:,:,:)
       emis_mass(region,mode_aii)%d4(:,:,:,2) = &
       emis_mass(region,mode_aii)%d4(:,:,:,2) + &
            (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  )
       deallocate(  work3d(region)%d3  )
    end do
    
    status = 0
    
  END SUBROUTINE EMISSION_POM_DECLARE
  !EOC

END MODULE EMISSION_POM