ecearth3/sources/tm5mp/proj/cb05/emission_bc.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_BC
!
! !DESCRIPTION: Data and methods fo Black Carbon (BC) emissions.
!
!>>> TvN
! The latest revision allows to apply
! different sizes for the emissions from
! different sectors, as well as for
! biofuel emissions.
! The assumption that all BC 
! is emitted in the insoluble Aitken mode
! has been relaxed.
! See settings in emission_data.F90.
!<<< TvN
!
!\\
!\\
! !INTERFACE: 
!
MODULE EMISSION_BC
  !
  ! !USES:
  !
  use GO,            only : gol, goErr, goPr
  use tm5_distgrid,  only : dgrid, get_distgrid, scatter, gather
  use dims,          only : nregions, okdebug
  use global_data,   only : emis_data, d3_data
  use emission_data, only : emis_input_year_bc
  use emission_read, only : used_providers_aer, has_aer_emis

  IMPLICIT NONE
  PRIVATE
  !
  ! !PUBLIC MEMBER FUNCTIONS:
  !
  public :: emission_bc_init, emission_BC_done, emission_bc_declare
  !
  ! !PRIVATE DATA MEMBERS:
  !
  character(len=*), parameter  ::  mname = 'emission_bc'

  type(emis_data), dimension(:,:), allocatable :: bc_emis_2d, bc_bf_emis_2d
  type(d3_data),   dimension(:,:), allocatable :: bc_emis_3d

  integer                           :: bc_2dsec, bc_3dsec
  !
  ! !REVISION HISTORY: 
  !    ? ??? 2005 - Elisabetta Vignati - changed for coupling with M7
  !    1 Sep 2010 - Achim Strunk - revised for AR5 emissions
  !   22 Jun 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
  !
  ! !REMARKS:
  !
  !EOP
  !------------------------------------------------------------------------

CONTAINS

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    EMISSION_BC_INIT
  !
  ! !DESCRIPTION: Allocate space needed to handle the emissions.
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE EMISSION_BC_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 -
    !   25 Jun 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

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

    status = 0
    if(.not. has_aer_emis) return

    ! nb of sectors
    bc_2dsec = 0
    bc_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
             bc_2dsec = bc_2dsec + n_ar5_ant_sec*count('BC'.eq.ar5_cat_ant) + &
                                   n_ar5_shp_sec*count('BC'.eq.ar5_cat_shp)

             if (LAR5BMB) bc_2dsec = bc_2dsec + n_ar5_bmb_sec*count('BC'.eq.ar5_cat_bmb)

             bc_3dsec = bc_3dsec + count('BC'.eq.ar5_cat_air)
          else
             bc_2dsec = bc_2dsec + providers_def(lprov)%nsect2d
             bc_3dsec = bc_3dsec + providers_def(lprov)%nsect3d  
          endif
       endif
    enddo

    allocate( bc_emis_2d( nregions, bc_2dsec ) )
    allocate( bc_bf_emis_2d( nregions, bc_2dsec ) )
    allocate( bc_emis_3d( nregions, bc_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,bc_2dsec
          allocate( bc_emis_2d(region,lsec)%surf(i1:i2,j1:j2) )
          allocate( bc_bf_emis_2d(region,lsec)%surf(i1:i2,j1:j2) )
       end do

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

    enddo

    ! ok
    status = 0

  END SUBROUTINE EMISSION_BC_INIT
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    EMISSION_BC_DONE
  !
  ! !DESCRIPTION: Free space after handling of BC emissions.
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE EMISSION_BC_DONE( status )
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)          ::  status
    !
    ! !REVISION HISTORY: 
    !    1 Oct 2010 - Achim Strunk -
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

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

    status = 0
    if(.not. has_aer_emis) return

    do region = 1, nregions
       do lsec=1,bc_2dsec
          deallocate( bc_emis_2d(region,lsec)%surf )
          deallocate( bc_bf_emis_2d(region,lsec)%surf ) 
       end do
       do lsec=1,bc_3dsec
          deallocate( bc_emis_3d(region,lsec)%d3 )
       end do
    end do

    deallocate( bc_emis_2d )
    deallocate( bc_bf_emis_2d )
    deallocate( bc_emis_3d )

    ! ok
    status = 0

  END SUBROUTINE EMISSION_BC_DONE
  !EOC


  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    EMISSION_BC_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_BC_DECLARE( status )
    !
    ! !USES:
    !
    use binas,         only : pi
    use partools,      only : isRoot, par_broadcast
    use toolbox,       only : coarsen_emission
    use dims,          only : im, jm, lm, idate, sec_month, nlon360, nlat180
    use chem_param,    only : xmbc, sigma_lognormal, carbon_density
    use chem_param,    only : mode_aii, mode_ais, mode_acs
    use emission_data, only : emis_mass, emis_number, msg_emis, LAR5BMB
    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 : frac_bc_sol_ff, frac_bc_sol_bf, frac_bc_sol_bb
    use emission_data, only : emission_vdist_by_sector
    ! ----------------  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 : sectors_def, numb_sectors
    use emission_read, only : ar5_dim_3ddata
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)        ::  status
    !
    ! !REVISION HISTORY: 
    !    1 Oct 2010 - Achim Strunk - adapted for AR5
    !   22 Jun 2012 - P. Le Sager  - adapted for lon-lat MPI domain decomposition
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

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

    real,dimension(:,:,:),  allocatable  :: field3d
    real,dimension(:,:),    allocatable  :: field2d_bf
    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 
    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
    integer                              :: seccount2d, seccount3d

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

    status = 0
    if(.not. has_aer_emis) return
    
    write(gol,'(" EMISS-INFO ------------- read BC emissions -------------")'); call goPr
    
    ! Reset emissions
    do region = 1, nregions
       do lsec=1,bc_2dsec
          bc_emis_2d(region,lsec)%surf = 0.0
          bc_bf_emis_2d(region,lsec)%surf = 0.0 
       end do
       do lsec=1,bc_3dsec
          bc_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


    ! --------------------------------
    ! 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)*carbon_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)*carbon_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)*carbon_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 ! skip unused providers
       
       if (associated(sectors_def(lsec)%species)) then  ! AR5 check
          if (count('BC'.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 (isRoot) then ! READ

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

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

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

          case( 'AR5' )

             ! Screen out agricultural and solvent sectors for BC, 
             ! 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( 'BC', emis_input_year_bc, 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 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, 'BC', emis_input_year_bc, 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, 'bc', emis_input_year_bc, idate(2), &
                                            'GFED', 'kg / s', field3d(:,:,1), status )
             IF_NOTOK_RETURN(status=1)

          case('RETRO')
             call emission_retro_ReadSector( emis_input_dir_retro, 'BC', emis_input_year_bc, 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 BC"; 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 BC 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 BC emissions for ",a," ",a," for month ",i2 )') &
                     trim(sectors_def(lsec)%prov), trim(sectors_def(lsec)%name), idate(2) ; call goPr
             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)
          ! ---------------------------

          if (isRoot) then

             ! write some numbers
             call msg_emis( amonth, trim(sectors_def(lsec)%prov), sectors_def(lsec)%name, 'BC', xmbc, sum(field3d) )

             ! distribute to work arrays in regions
             call Coarsen_Emission( 'BC '//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 )

             bc_emis_3d(region,seccount3d)%d3 = bc_emis_3d(region,seccount3d)%d3 + emis3d(region)%d3
          end do

          ! add to emis target arrays (scale mass to get number)
          do region = 1, nregions
             emis_mass  (region,mode_aii)%d4(:,:,:,1) = &
             emis_mass  (region,mode_aii)%d4(:,:,:,1) + &
               (1.-frac_bc_sol_ff) * bc_emis_3d(region,seccount3d)%d3(:,:,:)
             emis_mass  (region,mode_ais)%d4(:,:,:,2) = &
             emis_mass  (region,mode_ais)%d4(:,:,:,2) + &
               (   frac_bc_sol_ff) * bc_emis_3d(region,seccount3d)%d3(:,:,:)

             emis_number(region,mode_aii)%d4(:,:,:,1) = &
             emis_number(region,mode_aii)%d4(:,:,:,1) + &
               (1.-frac_bc_sol_ff) * bc_emis_3d(region,seccount3d)%d3(:,:,:) * mass2numb_air_insol
             emis_number(region,mode_ais)%d4(:,:,:,2) = &
             emis_number(region,mode_ais)%d4(:,:,:,2) + &
               (   frac_bc_sol_ff) * bc_emis_3d(region,seccount3d)%d3(:,:,:) * mass2numb_air_sol 
          end do

       else ! 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, 'BC', xmbc, sum(field3d(:,:,1)) )

             call coarsen_emission( 'BC '//sectors_def(lsec)%name, nlon360, nlat180,   field3d(:,:,1), &
                                                                   wrk2D,   add_field, status )
             IF_NOTOK_RETURN(status=1)

             call coarsen_emission( 'BC 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), bc_emis_2d(region,seccount2d)%surf, work(region)%d3(:,:,1), 0, status)
             IF_NOTOK_RETURN(status=1)

             call scatter(dgrid(region), bc_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, 'BC', region, bc_emis_2d(region,seccount2d), &
                                            emis3d(region), status )
             IF_NOTOK_RETURN(status=1)

             if( trim(sectors_def(lsec)%catname) == 'biomassburning' ) then
               
                emis_mass  (region,mode_aii)%d4(:,:,:,1) = &
                emis_mass  (region,mode_aii)%d4(:,:,:,1) + emis3d(region)%d3(:,:,:) * &
                  (1.-frac_bc_sol_bb)

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

                emis_mass  (region,mode_acs)%d4(:,:,:,2) = &
                emis_mass  (region,mode_acs)%d4(:,:,:,2) + emis3d(region)%d3(:,:,:) * &
                  frac_bc_sol_bb

                emis_number(region,mode_acs)%d4(:,:,:,2) = &
                emis_number(region,mode_acs)%d4(:,:,:,2) + emis3d(region)%d3(:,:,:) * &
                  frac_bc_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 ( bc_emis_2d(region,seccount2d)%surf(:,:) > tiny(0.0) )
                     frac_bf(region)%d3(:,:,1) = bc_bf_emis_2d(region,seccount2d)%surf(:,:) / &
                                                 bc_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

               endif

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

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

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

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

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

                   emis_number(region,mode_acs)%d4(:,:,l,2) = &
                   emis_number(region,mode_acs)%d4(:,:,l,2) + emis3d(region)%d3(:,:,l) * &
                     (     frac_bf(region)%d3(:,:,1)  * frac_bc_sol_bf * mass2numb_bf_sol )

               enddo

             endif

          end do

       end if

    end do sec


    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( emis3d(region)%d3 )
       deallocate( frac_bf(region)%d3 )
    end do

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

    ! ok
    status = 0

  END SUBROUTINE EMISSION_BC_DECLARE
  !EOC

END MODULE EMISSION_BC